Hemodynamic and perfusion analyses are currently performed by injecting an intravenous bolus with gadolinium (Gd), and the lung is imaged using 4D magnetic resonance imaging (MRI) with the signal variation observed at the input (blood vessels) and at the output (lung tissue). However, 2 problems remain to be solved: (1) contrast calibration and (2) blood reperfusion.
Quantification of the Gd tracer concentration is an essential step and needs to be determined at each point to allow hemodynamic analysis. A simplification procedure assumes that the change in signal intensity or relative signal intensity because of the tracer is directly proportional to the tracer concentration.^{1,2} However, the effect of nonlinearity is particularly large in areas of high concentration,^{3} such as the input position in perfusion analysis, and any input component may be overestimated. The proposed is a calibration procedure using phantoms, which works even for high concentrations of the contrast agent, or regions where the contrast characteristic reverses.^{4} Some issues of the proposed method were found. When a highconcentration contrast agent was injected into the patient, the image quality was reduced and some artifacts occurred. The phantoms deteriorate over time, and it was necessary to use new phantoms for each imaging procedure. These issues motivated us to propose a nonlinear calibration method based on the minimum variation of the mass of the contrast agent, which does not require phantoms.
The pulmonary arteries leave the heart and carry deoxygenated blood to the lungs. After oxygenation of the blood, the blood leaves the lungs and returns to the heart. Most of the oxygenated blood goes to the body, and a small part enters the lungs to keep them alive. This phenomenon will be called lung reperfusion. Lung reperfusion is generally considered to be approximately 1%^{5} of the whole blood flow (BF) and is neglected by conventional methods.^{1,6} However, lung reperfusion increases and cannot be neglected in cases of lung cancer.^{7} The multiple input model proposed here considers lung reperfusion and its applications for analyzing lung cancer and performing functional and nutritional vasculature analysis separately. Although 2input perfusion has been proposed in the past,^{8} they used the maximum slope method, which assumes rapid intravenous injection of a contrast agent, and depending on the patient's condition, this is not true. The multipleinput analysis method proposed here has characteristics that do not depend on the input waveform.
MATERIALS AND METHODS
Calibration of Gd Tracer Concentration and Signal Intensity in MRI
In this work, lower levels of Gd tracer concentration are used, and the proposed tracer calibration procedure does not require any phantom. Figure 1 presents the route of the blood with Gd. The Gd concentrations in the pulmonary artery, left atrium, and aorta are represented as ${C}_{art}\left(t\right)$, ${C}_{la}\left(t\right)$, and ${C}_{aor}\left(t\right)$. The BF term is the blood flow in that vessel, and the M term is the total mass of the contrast tracer in the first perfusion of that vessel.
Initially, contrast blood enters the right heart through the vena cava and pumps to the lungs through the pulmonary artery with BF. Blood is oxygenated in the lungs and goes to the left atrium through the pulmonary vein. This oxygenated blood has 3 routes: aorta, coronary artery, and bronchial artery. Most goes to the entire body through the aorta, but some of it goes back to the lungs through the bronchial artery. The remainder goes to the heart through the coronary artery.
Now, blood from the bronchial artery becomes deoxygenated after passing through the lungs, and there are 2 possible routes: the bronchial and pulmonary veins. Blood that goes through the pulmonary veins to the left heart has a small portion of deoxygenated blood (lung shunt) combined with oxygenated blood. Based on this discussion, ${C}_{art}\left(t\right)$ is the first input of blood to the lungs and ${C}_{aor}\left(t\right)$ is the second input.
Note that there are 3 circulatory pathways. The first is a circulatory path going through the pulmonary artery and has $BF$. This circulatory path will be called circulatory path of pulmonary artery (CPPA). The second is a circulatory path with BF equal to $l\times BF$ and represents the pulmonary shunt. This circulatory path will be called circulatory path of pulmonary shunt (CPPS). The third is a circulatory path going through the coronary artery and has BF equal to $n\times BF$. This circulatory path will be called CPCA.
Considering exclusively CPPS, after infinite iterations, the amount of contrast mass that passes through the left atrium can be determined by(1)$\begin{array}{l}M+\frac{l}{1+l}M+{\left(\frac{l}{1+l}\right)}^{2}M+\cdots =\\ {\displaystyle \sum}_{n=0}^{\infty}M{\left(\frac{l}{1+l}\right)}^{n}=\left(1+l\right)M.\end{array}$
Because the BF in the left atrium is $\left(1+l\right)BF$, the amount of contrast going through the left atrium is as follows:(2)$\begin{array}{l}\left(1+l\right)M=\left(1+l\right)BF{\int}_{0}^{\infty}{C}_{la}\left(t\right)dt\\ M=BF{\int}_{0}^{\infty}{C}_{la}\left(t\right)dt\end{array}$
The aorta undergoes a similar procedure (considering exclusively CPPS), resulting in(3)$\begin{array}{l}\left(1mn\right)M=\left(1mn\right)BF{\int}_{0}^{\infty}{C}_{aor}\left(t\right)dt\\ M=BF{\int}_{0}^{\infty}{C}_{aor}\left(t\right)dt\end{array}$
A similar procedure can be applied to the bronchial artery. From the above, the integration of C(t) in time is constant, independent of the part, and can be written as(4)${\int}_{0}^{\infty}C\left(t\right)dt=\frac{M}{BF}$
Eq. (4) is also valid for the pulmonary artery as all contrast passes through it. This result considers an ideal world, where no errors occur. In the real world, it is necessary to deal with errors and the variation of the integration of C(t) in time is considered minimum.
In this research, the contrast concentration used is low and preserves a nonlinear and monotonic relationship with the intensity of the MRI signal. Therefore, we propose a calibration using(5)$C\left(t\right)=\alpha \left\{\mathit{exp}\left[\beta \left(\frac{SI\left(t\right)}{SI\left(0\right)}1\right)\right]1\right\}$
The values of α and β are determined by minimizing the tracer concentration in the following compartments: the pulmonary artery, the aorta, and the left atrium. Minimization is solved by the LevenbergMarquardt method. The function to be minimized is(6)$\underset{\alpha ,\beta}{\mathrm{arg}\hspace{0.17em}\mathrm{min}}={\left(\int {C}_{aor}\int {C}_{art}\right)}^{2}+{\left(\int {C}_{la}\int {C}_{art}\right)}^{2}$
Hemodynamic Analysis
Theory
The indicator dilution method was established as a theory for hemodynamic analysis,^{9,10} and each voxel can be used to analyze a section and determine regional blood volume (rBV; mL_{Blood}/mL_{Voxel}), regional mean transit time (rMTT; s) and regional blood flow (rBF; mL_{Blood}/mL_{Voxel}/s).^{11}(7)$rBV=\frac{\underset{0}{\overset{\infty}{{\displaystyle \int}}}{C}_{tis}\left(t\right)dt}{\underset{0}{\overset{\infty}{{\displaystyle \int}}}{C}_{art}\left(t\right)dt}$(8)$rMTT=\underset{0}{\overset{\infty}{{\displaystyle \int}}}h\left(t\right)tdt$(9)$rBF=\frac{rBV}{rMTT}$where ${C}_{tis}\left(t\right)$ is the average contrast concentration that occupies the tissue inside a voxel. $h\left(t\right)$ is the probability distribution function of the tracer transit time, and using the outflow concentration ${C}_{out}\left(t\right)$, it is given by(10)$h\left(t\right)=\frac{{C}_{out}\left(t\right)}{\underset{0}{\overset{\infty}{{\displaystyle \int}}}{C}_{out}\left(t\right)dt}.$
Because outflow cannot be observed in voxels, the residue function of the contrast inside the voxel R(t) is used instead.(11)$R\left(t\right)=rBF\hspace{0.17em}\left[1\underset{0}{\overset{\infty}{{\displaystyle \int}}}h\left(\tau \right)d\tau \right]$${C}_{tis}\left(t\right)$ can be evaluated from the convolution of the accumulated amount of the tracer inside the voxel, defined by(12)${C}_{tis}\left(t\right)\left[mmol/{ml}_{Blood}\right]=\underset{0}{\overset{t}{{\displaystyle \int}}}{C}_{in}\left(\tau \right)R\left(t\tau \right)d\tau .$
Because the impulse input is impractical, the r(t) estimation must be performed using an inverse filter. The inverse filter is sensitive to noise, so γvariate function was fitted to the input and output waveforms.^{1,6} However, this ignores CPPS.(13)$\gamma \left(t\right)=a{\left(tb\right)}^{c}\hspace{0.17em}\mathit{exp}\left[\frac{\left(tb\right)}{d}\right]$
Multiple Input Method
We used a γvariate function Eq. (13) to model the impulse response, instead of fitting γvariate functions to the input and output waveforms. Substituting Eq. (4) and (8) into expression Eq. (11), it is possible to obtain a new expression to R(t) by(14)$R\left(t\right)=rBF\left\{\hspace{0.17em}1\underset{0}{\overset{t}{{\displaystyle \int}}}\left[\frac{\gamma \left(\tau \right)}{\underset{0}{\overset{\infty}{{\displaystyle \int}}}\gamma \left(t\right)dt}\right]d\tau \right\}$
The γvariate function parameters that minimize the difference(15)$\underset{\left\{a,b,c,d\right\}}{\mathrm{arg}\hspace{0.17em}\mathrm{min}}{\displaystyle \sum}_{t=0}^{n}{\left[{C}_{tisobs}\left(t\right){C}_{tiscal}\left(t\right)\right]}^{2}$between the observed ${C}_{tis\_obs}\left(t\right)$ and the calculated ${C}_{tis\_cal}\left(t\right)$ values are obtained by using the LevenbergMarquardt estimation.^{12}${C}_{tis\_cal}\left(t\right)$ can be obtained from ${C}_{in}\left(t\right)$ and R(t) from Eq. (12) and (14). This is an inverse problem, and the initial values influence the final solution. The formulation of the inverse problem can be greatly modified to incorporate meaningful variables from which good initial values can be obtained.
Instead of using the γvariate function from Eq. (13), it is possible to use the incomplete Γ function and regularized upper incomplete Γ function (Qfunction) defined by^{13}(16)$\Gamma \left(a,x\right)={\int}_{x}^{\infty}{t}^{a1}{e}^{t}dt$(17)$\mathrm{Q}\left(\mathrm{a},\mathrm{x}\right)=\frac{\Gamma \left(a,x\right)}{\Gamma \left(a\right)}$and after some algebraic manipulation, r(t) can be represented by(18)$R\left(t\right)=\{\begin{array}{c}rBF\times Q\left(\frac{{\sigma}^{2}}{{d}^{2}},\frac{trMTT}{d}+{\sigma}^{2}\right),\hspace{0.17em}\hspace{0.17em}t\ge rMTT\\ rBF,\hspace{0.17em}otherwise\end{array}$where(19)${\sigma}^{2}={\int}_{0}^{\infty}{\left(trMTT\right)}^{2}h\left(t\right)dt={d}^{2}\left(c+1\right)$
Expression Eq. (18) is the new singleinput method with meaningful variables. When considering multiple inputs, we assign a subscript index to each input ${C}_{in,i}\left(t\right)$ and${R}_{i}\left(t\right)$. ${C}_{tis\_cal}$ expressed by Eq. (12), can be evaluated by accumulating ${C}_{in,i}\left(t\right)$ and ${R}_{i}\left(t\right)$, as(20)${C}_{tis\_cal}\left(t\right)={\displaystyle \sum}_{i}{\int}_{0}^{t}{C}_{in,i}\left(\tau \right){R}_{i}\left(t\tau \right)d\tau $
Expressions Eq. (18) and (20) are used to minimize(21)$\underset{\left\{{rBF}_{i},{rMTT}_{i},{\sigma}_{i},{d}_{i}\right\}}{\mathrm{arg}\hspace{0.17em}\mathrm{min}}{\displaystyle \sum}_{t=0}^{n}{\left[{C}_{tis\_obs}\left(t\right){C}_{tis\_cal}\left(t\right)\right]}^{2}$
The initial values for ${rBF}_{i},$${rMTT}_{i}$ and ${\sigma}_{i}$ can be obtained, and the convergence to an unintended local solution can be reduced. rBV can be determined with Eq. (9) where rBF and rMTT are already known.
The proposed method is not influenced by reperfusion or noise generated by aliasing. Therefore, it is possible to directly use the input and output waveforms without performing any curve fitting.
Experiment
The Spoiled Gradient Recalled Acquisition in Steady State (SPGR) method was used with a highspeed time series MRI of the chest (12 [slice], 12.0 [mm/slice], 1.62 [s/frame], 14 [frame], 256 × 256 [pixels], 1.68 [mm/pixels], coronal). Table 1 presents the device information and imaging parameters. The contrast agent Gd was adjusted to 0.05 mmol/kg according to the weight of the patient. The capture began concurrently with the subject holding his breath, and the contrast agent was injected immediately afterward. The analyses were performed exclusively on frames in the breathhold state. The results shown here are related to 6 patients (see Table 2 for details).
Table 1 
Equipment and Parameters, the Only Difference Among the Patients' Parameters Used is the TR
MRI Type 
Achieva 1.5T ASeries (Philips) 
Type of Gd contrast 
Magnevist 
Bolus administration 
5 mL/seconds flush with saline 
TE 
0.796 ms 
TR 
2349 ms (A) 2677 ms (BF) 
FA 
25° 
Number of phase encoding steps 
132 
Protocol name 
Perfusion SENSE 
Gd relaxivity 
2500 mL/sec/mmol 
CT type 
Aquilion ONE (TOSHIBA) 
ECG triggering 
Not used 
Repetition time was different for patient A.
FA indicates flip angle; Gd, gadolinium; MRI, magnetic resonance imaging; TE, echo time; TR, repetition time.
Table 2 
Patient Characteristics
Patient 
Sex 
Disease 
A 
Male 
Lung cancer 
B 
Male 
Lung cancer 
C 
Male 
Lung cancer 
D 
Male 
Lung cancer 
E 
Male 
Lung cancer 
F 
Female 
Lung cancer 
Two calibration and BF analysis experiments were performed. In both experiments, the pulmonary artery, aorta, and left atrium were all specified as circles with a radius of 3 pixels. Then, the average signal intensity for each region was determined.
Calibration Experiment
The proposed calibration is compared with the conventional linear approximation and phantom calibration. In the phantom calibration,^{4} 11 contrast agent phantoms were placed near the subject's shoulder. Each phantom has a different value of Gd concentration, and its position and dimension in the image are known. The signal intensity associated with each phantom is determined by manually drawing a circle with a radius of 3 pixels. The associated pixel intensity average and standard deviation are determined for each phantom.
To validate the proposed calibration model represented by Eq. (5), we also compared it with four dimensional computerized tomography. The phantom study showed that the calibration of the iodine contrast concentration and the intensity of the CT signal present a high proportional relationship.^{14} The imaging conditions were (320 [slice], 0.5 [mm/slice], 2.05 [s/frame], 16 [frame], 512 × 512 [pixels], and 0.625 [mm/pixels]) of a healthy female subject.
Our BF Analysis
Multiinput BF analysis was performed to confirm the results of ${C}_{tis}\left(t\right)$ waveform optimization and to confirm the images analyzed. The lung region is determined by the user manually. A percentile filter of the top 5% of signal intensity was used to remove the heart and some thick vessels. The initial values for ${rBF}_{i}$, ${rMTT}_{i}$ and ${\sigma}_{i}$ are determined, and the LevenbergMarquardt method is used to solve Eq. (21). A color mask shows the BF analysis results.
To evaluate the amount of CPPS perfusion, the difference between the left and right cardiac output is calculated and compared. It is obtained using a Fast GE Cine. The diastolic and systolic areas of the right and left ventricles are segmented, and the difference between systolic output per stroke is calculated.^{15} The imaging parameters used are followed, and the epicardial and endocardial contours of the basal and apical slices were outlined manually.
Another approach to evaluate CPPS is to integrate the rBV from each input contribution: pulmonary artery and aorta. Both the conventional linear approximation and the proposed calibration were performed, and $\int {rBV}_{art}$ and $\int {rBV}_{aor}$ were compared.
RESULTS
Hemodynamic analysis Experiment
Figure 2 shows a comparison between the 3 calibration procedures: the conventional linear approximation, the calibration procedure defined by Eq. (5), and the phantom calibration. For pixels with higher intensity, the contrast concentration is underestimated by the conventional linear approximation compared with the one proposed here. Considering the phantom and linear calibrations, the linear model deviates from the phantom calibration at higher Gd concentrations up to 1.5 × 10^{2}.
Figure 3 shows CT values in the pulmonary artery, left atrium, and aorta. The values of $\int {C}_{art}\left(t\right)dt$ are calculated for both calibration procedures, and they are considered reference values (the value for both calibration procedures is 1.00). If well calibrated, such as CT, these values were shown to be equal within 2% error.
Figure 4 shows the tracer concentrations in the same 3 regions as in Figure 3 taken by MRI. The value of $\int {C}_{la}\left(t\right)dt$ is 16% higher, and $\int {C}_{aor}\left(t\right)dt$ is 15% higher than the reference value when using the conventional linear approximation. With the proposed calibration, $\int {C}_{la}\left(t\right)dt$ is 2% and $\int {C}_{aor}\left(t\right)dt$ is 0% higher than the reference. The proposed model showed a smaller difference in the integral values of the pulmonary artery, left atrium, and aorta than the linear model.
Table 3 presents the results and averages for the 6 patients. In the linear approximation, the left atrium and aorta were 13% and 10% larger than the pulmonary artery, but these were suppressed to 5% and 1% by the proposed method.
Table 3 
Results for the Calibration Obtained for the 6 Patients Considering
Eq. (5) and the Conventional Linear Calibration
Patients 
(5)
$C\left(t\right)=\alpha \left\{\mathit{exp}\left[\beta \left(\frac{SI\left(t\right)}{SI\left(0\right)}1\right)\right]1\right\}$

Linear 
Art 
L Atr 
Aor 
Art 
L Atr 
Aor 
A 
1.00 
1.02 
1.00 
1.00 
1.16 
1.15 
B 
1.00 
1.09 
1.01 
1.00 
1.16 
1.11 
C 
1.00 
1.16 
0.94 
1.00 
1.21 
1.14 
D 
1.00 
1.01 
0.99 
1.00 
1.06 
1.03 
E 
1.00 
1.00 
1.00 
1.00 
1.07 
1.04 
F 
1.00 
1.04 
0.98 
1.00 
1.10 
1.12 
Average 
1.00 
1.05 
0.99 
1.00 
1.13 
1.10 
Multiple Input BF Analysis Results
An example of a patient A, with the contribution of each input to the pulmonary perfusion, is shown in Figure 5. The contribution from each input can be analyzed separately because of the time difference. The peak is reached around 6.5 seconds, and the correspondence between ${C}_{art}\left(t\right)\u229b{R}_{art}\left(t\right)$ and the observed data is good considering the time interval without the influence of the second input, the aorta contribution. Near 13 seconds, the contribution from the pulmonary artery reaches approximately 0 mmol/mL, and the correspondence between ${\hspace{0.17em}C}_{aor}\left(t\right)\u229b{R}_{aor}\left(t\right)$ and the observed data has good correspondence. By summing the contributions from both inputs, a good correspondence with the observed data is obtained. There is only 1 possible outlier near the time instant of 19 seconds.
Figure 6 shows the results of the local analysis. In particular, (Figure 6F) and (Figure 6G) show that the blood contribution of the bronchial artery is greater and supports the cancer cells receiving arterial blood.
Table 4 presents the results for left and right cardiac output. The stroke volume for 1 cardiac output and the cardiac output volume per second are estimated according to the estimated heart rate. The difference between the cardiac output is 4.3 mL/s and represents CPPS and 4.3 mL/s represents 6.2% of the output from the right ventricle, which is 69.2 mL/s. CPPS has 6.2% of the whole BV. CPPS is the l parameter used in expression Eq. (1).
Table 4 
Cardiac Output Using Fast GE Cine
 Left Ventricle 
Right Ventricle 
Heart rate 
81 bpm 
Stroke volume 
54.4 mL 
51.2 mL 
Cardiac output 
73.5 mL/s 
69.2 mL/s 
Difference leftright 
4.3 mL/s 
CPPS perfusion 
6.2% 
The same subject with lung cancer is considered.
Table 5 presents the results for $\int {rBV}_{art}$ and $\int {rBV}_{aor}$ considering both calibration procedures: the conventional linear approximation and the one proposed here. For both calibration methods, the CPPS was greater than 20%.
Table 5 
Blood Volume Contribution of Each Input and Comparison Between the 2 Calibration Procedures
Calibration 
$\int {rBV}_{art}$

$\int {rBV}_{aor}$

CPPS Perfusion 
Linear 
318.0 mL 
74.7 mL 
23.5% 
(5)
$C\left(t\right)=\alpha \left\{\mathit{exp}\left[\beta \left(\frac{SI\left(t\right)}{SI\left(0\right)}1\right)\right]1\right\}$

303.3 mL 
62.4 mL 
20.5% 
DISCUSSION
Calibration Model Consideration
As presented in Table 3, the proposed calibration still produces errors in $\int {C}_{la}\left(t\right)dt$ and $\int {C}_{aor}\left(t\right)dt$. This might originate from a different source of errors: imaging errors (such as patient movement during acquisition and segmentation problems), the inverse characteristic of contrast agents, and others.
The first is the imaging error. Considering Figure 3, even in CT, $\int {rBV}_{la}$ has a value 1% higher than the reference and $\int {rBV}_{aor}$ has a value 2% higher than the reference. The sampling rate was cited as the cause of this: 2.05 seconds/frame for CT and 1.62 seconds/frame for MRI. This low sampling rate may have caused errors in the integrated values.
The second is the inverse characteristic of contrast agents. When the concentration of the MRI contrast agent is too high, the contrast characteristic is reversed and the intensity of the MR signal reduces, as shown in Figure 2. In this study, the reverse characteristic was not considered. Because the contrast agent concentration reached near 0.03 mmol/mL in the pulmonary artery, the effect is not serious. However, a calibration procedure that includes this feature is desired and can be used in future experiments.
Difference Between Left and Right Cardiac Output
In this section, we discuss the error between the CPPS obtained from the difference between the left and right cardiac output of patient A (6.2%) and the CPPS obtained by the proposed multipleinput analysis (20.5%) The discrepancy between both CPPS is 14.3%.
The first is the influence of thick blood vessels near the lung boundary. Figure 7 shows the results of a local analysis ${rBV}_{aor}$ in a thick slice. It is possible to observe an overestimation of BV values in the region near the spinal cord. One possibility is the inclusion of the aorta in the region of interest (ROI). Figure 8 also shows this overestimation for a specific voxel. The concentration associated with the second input increases greatly, and one possibility is that the parenchymal tissue, such as the intercostal muscle, has reached the thickness of the slice of 12 mm.
In this additional experiment, a threshold value was set to exclude thick blood vessels from the ROI. The threshold cannot be too reduced because it is necessary to detect lung cancer, and consequently, some large vessels are included in the ROI. Trying to overcome the influence of some thick blood vessels near the lung boundary, a morphology contraction operation was performed and ROI was reduced by 3 pixels. Table 6 presents the results with a 3pixel contraction; the ratio between $\int {rBV}_{art}$ and $\int {rBV}_{aor}$ improved by approximately 7%, from 20.5% to 13.5%. However, because the lung boundary BV was not considered for pulmonary artery input, the $\int {rBV}_{art}$ value decreased by 109.9 mL. Blood oxygenates at the lung boundary through the pulmonary alveolus, and a higher contrast concentration is expected. Considering this fact, a better proposal is required to define ROI, in which the contribution of each input of contrast concentration is considered. A higher concentration of the first input represents blood oxygenation near the lung boundary, and a higher concentration of the second input represents oxygenated blood going through thick vessels. Figure 8 combines the internal part of the lung with a high peak originating from the pulmonary artery and the external part of the lung with a high peak of oxygenated blood originating from the aorta. The cause of this combination may be the proximity between the lung boundary and a thick vessel and/or the occurrence of internal involuntary movement.
Table 6 
Contribution of Blood Volume From Each Input Without Boundary Contraction (0 pixels) and With Boundary Contraction (3 pixels)
Contraction 
$\int {rBV}_{art}$

$\int {rBV}_{aor}$

CPPS Perfusion 
0 pixels 
303.3 mL 
62.4 mL 
20.5% 
3 pixels 
193.4 mL 
26.1 mL 
13.5% 
The second is that we consider the influence of the CPCA (coronary artery) as an overestimation of the perfusion of the bronchial artery. Figure 9 shows the influence of 3 circulatory paths (CPPA, CPPS, and CPCA) at 3 different locations (pulmonary artery, lungs, and aorta). The contrast is injected and goes through the pulmonary artery and then through the lungs and aorta. This contribution comes from CPPA and is represented by a blue curve. When it reaches the aorta, the BF splits into 3 different circulatory paths. Two of them are represented here: CPPS and CPCA. CPPS enters the lungs and then returns to the aorta. It is represented by a red curve. CPCA enters the heart and then enters the pulmonary artery, lungs, and aorta. It is represented by a black curve. In CPCA, it takes approximately 8 ± 2 seconds to go from the coronary artery to the coronary vein.^{16} However, the subject can hold his/her breath for a limited time and starts breathing before coronary circulation is observed in the lung. Therefore, it is inferred that the acquisition of the contribution from CPCA was interrupted, causing an overestimation of blood perfusion in the pulmonary artery. Considering Eq. (5), the perfusion from CPCA was included in ${C}_{in}$, but it was not included in ${C}_{tis}$; this situation caused a relative overestimation of the pulmonary artery.
Figure 4 shows a small increase in perfusion in the pulmonary artery near 19 seconds, approximately 8 seconds after the peak of aorta perfusion (the peak occurred around 11 seconds). This is approximately the time that the blood must travel from the coronary artery to the coronary vein.^{16} The same happens in Figure 3, where a second peak appears near 26 seconds that is 10 seconds after the peak of aorta perfusion (the peak occurred around 16 seconds). The time difference indicates that the additional perfusion appearing in Figures 3 and 4 can be derived from CPCA.
An estimation of the contribution of the coronary artery can be performed by displacing the aorta graph by 10 seconds; this way, the peak of the aorta graph will be placed near 26 seconds. The translated graph is scaled to coincide with the peak near 26 seconds on the pulmonary artery curve. The resulting graph is shown in yellow in Figure 10. This is a virtual coronary contribution to the pulmonary artery. The $\int {C}_{cor}$ (t)dt in the pulmonary artery has 13% of BV in the pulmonary artery. The value of 13% is an overestimate that could include the perfusion of some organs near the lung, such as the liver. CPCA is known to have approximately 5% of the right cardiac output.
The last is the reperfusion error. Eq. (1) considers an infinite number of iterations. However, only 2 iterations are considered in the experiments. It is necessary to evaluate the error present in the calibration performed considering only 2 iterations. The results presented in Table 4 show that l = 6.2%. Considering just 2 iterations after truncating Eq. (1), we have$\frac{1+2l}{1+l}=1.058.$
Changing l = 6.2% in Eq. (1), after infinite iterations, we get 1.062, the error is approximately 0.4% and can be ignored.
CONCLUSIONS AND FUTURE WORKS
A new calibration method is proposed based on the circulatory system, considering the minimum variation of the integral of the contrast concentration in time. Because no phantom is required in the proposed calibration method, it is cheaper. This is a great advantage.
The proposed singleinput method is not influenced by reperfusion or noise generated by aliasing. The twoinput method allowed the separation of functional and nutritional vasculature analysis. Six subjects considered having lung cancer, and a large perfusion with oxygenated blood from the aorta was observed. Comparison of rBV integrals for each input contribution with the Fast GE Cine method showed a discrepancy of approximately 14%. This discrepancy might happen because of the inability to perform perfect lung segmentation. Some thick blood vessels are near the lung boundary and influence local analysis.
Future research is intended to improve the determination of ROI, removing higher oxygenated blood peaks near the lung boundary. Current BF analysis is limited by the subject's ability to hold his breath. For example, Eq. (5) should use not only 1 SI(0) but also the average of multiple images before contrast injection. Future research also intends to improve the registration procedure to allow the same analysis to occur without this limitation.
ACKNOWLEDGMENTS
We are indebted to all patients and especially to the residents who worked in our unit during the study period for their dedication and collaboration in providing care to the patients who participated in this research.
References
1. Ohno Y, Hatabu H, Higashino T, et al. Dynamic perfusion MRI versus perfusion scintigraphy: prediction of postoperative lung function in patients with lung cancer. AJR Am J Roentgenol. 2004;182:73–78.
2. Gunnar B, Semmler W, Port R, et al. Pharmacokinetic parameters in CNS GdDTPA enhanced MR imaging. J Comput Assist Tomogr. 1991;15:621–628.
3. Jeong HK, Lee KH, Kim MH, et al. Signal intensity of contrast enhancement according to TE in 3.0T MRI T1 imaging. Appl Sci. 2018;8(7):1138.
4. Saka T, Ichikawa M, Kagei S, et al. Perfusion analysis for lung MR images considering nonmonotonic response of Gdcontrast agent. IFAC Proc Vol. 2014;47(3):3587–3592.
5. Walker CM, RosadodeChristenson ML, MartínezJiménez S, et al. Bronchial arteries: anatomy, function, hypertrophy, and anomalies. Radiographics. 2015;35(1):32–49.
6. Thompson HK, Starmer CF, Whalen RE, Mcintosh HD. Indicator transit time considered as a gamma variate. Circ Res. 1964;14:502–515.
7. Buckley DL, Roberts C, Parker GJM, et al. Prostate cancer: evaluation of vascular characteristics with dynamic contrast enhanced T1weighted MR imaging—initial experience. Radiology. 2004;233:709–715.
8. Fang J, Li H, Liang M, et al. Pulmonary dual hemodynamic changes in severe copd patients: a quantitative study using lowdose ct lung perfusion scan. Iran J Radiol. 2019;16(1):1–7.
9. Meier P, Zierler KL. On the theory of the indicatordilution method for measurement of blood flow and volume. J Appl Physiol. 1954;6:731–744.
10. Doriot PA, Dorsaz PA, Dorsaz L, Rutishauser WJ. Is the indicator dilution theory really the adequate base of many blood flow measurement techniques? Med Phys. 1997;24:1889–1898.
11. Wirestam R, Knutsson L, Ostergaard L, et al. Assessment of regional cerebral blood flow by dynamic susceptibility contrast MRI using different deconvolution techniques. Magn Reson Med. 2000;43:691–700.
12. More JJ. The LevenbergMarquardt algorithm: implementation and theory. Numer Anal. 1978;630:105–116.
13. Saka T, Gotoh T, Kagei S, Iwasawa T. Blood flow contribution analysis for pulmonary artery and aorta using contrast enhanced images. Med Imaging Technol. 2016;34:245–258.
14. Chandarana H, Megibow AJ, Cohen BA, et al. Iodine quantification with dualenergy CT: phantom study and preliminary experience with renal masses. Am J Roentgenology. 2011;196:693–700.
15. Rominger MB, Bachmann GF, Pabst W, Rau WS. Right ventricular volumes and ejection fraction with fast cine MR imaging in breathhold technique: applicability, normal values from 52 volunteers, and evaluation of 325 adult cardiac patients. J Magn Reson Imaging. 1999;10:908–918.
16. Ishikawa K, Haneda T, Ikeda S, et al. Coronary circulation time in man. Tohoku J Exp Med. 1974;114:253–261.