|Home | About | Journals | Submit | Contact Us | Français|
The purpose of this research was to develop a novel numerical procedure to deconvolute arterial input function from contrast concentration vs. time curves and to obtain the impulse response functions from dynamic contrast enhanced MRI data. Numerical simulations were performed to study variations of contrast concentration vs. time curves and the corresponding impulse response functions. The simulated contrast media concentration curves were generated by varying the parameters of an empirical mathematical model within reasonable ranges based on a previous experimental study. The arterial input function was calculated from plots of contrast media concentration vs. time in muscle under assumption that they are well approximated by the two-compartment model. A general simple mathematical model of the impulse response function was developed and the physiological meaning of the model parameters was determined by comparing them with the widely accepted ‘two compartment model’. The results demonstrate that the deconvolution procedure developed in this research is a simple, robust, and useful technique. In addition, ‘impulse response analysis’ leads to the derivation of novel parameters relating to tumor vascular architecture, and these new parameters may have clinical utility.
Dynamic contrast enhanced MRI (DCEMRI) detects early cancer based on rapid uptake and washout of contrast media. The clinical importance of DCEMRI for detection of cancers is increasing – particularly for breast cancer screening and diagnosis (1,2). Regular MRI scans are recommended for people who are at high risk for breast cancer (3). Early detection of malignant breast tumors when they are small and relatively easy to treat may significantly decrease morbidity and mortality, and in addition significantly reduce the financial cost of treatment. However, methods currently available to analyze the DCEMRI data have poor specificity and yield many false positives (4). In addition, improvements are needed in sensitivity of DCEMRI to early pre-invasive cancers (5).
DCEMRI data are frequently analyzed clinically based on the change in signal-intensity as a function of time following contrast injection. The change in signal intensity is evaluated using qualitative parameters, for example using the Kuhl classification system (6), or using semi-quantitative parameters related to contrast media uptake and washout. More quantitative analysis is based on the change in contrast media concentration as a function of time (7) to reduce sensitivity to instrumental/pulse sequence characteristics. The kinetics of contrast media uptake and washout can be evaluated using ‘model-free’ approaches such as the ‘area under the curve’ (8), or the empirical mathematical model (EMM) (9) developed in this laboratory. Alternatively, physiological models can be used to extract parameters more directly related to perfusion/capillary permeability. The two-compartment model (TCM) is widely accepted (10) and is used to fit DCEMRI data and extract physiological parameters that are markers for disease. This model assumes that tissue can be modeled as two well-mixed compartments – a vascular and an extra-vascular extra-cellular compartment. However, tumors are generally heterogeneous, and are not well approximated by two well mixed compartments. As a result, the TCM often does not provide good fits to experimental data from tumors (9). Inaccurate fitting of the change in contrast media concentration as a function of time is a primary factor limiting diagnostic accuracy of DCEMRI.
A general problem with conventional analysis of DCEMRI data is that variation in the arterial input function (AIF) is not properly accounted for. Qualitative approaches such as the Kuhl classification method (6), and semi-quantitative analysis of signal intensity changes (11) do not take the AIF into account. This leads to errors because variations in contrast media injection, cardiac output, and major arterial supply to the tumor containing tissue, cause the AIF to vary significantly among patients (12). To account for these variations, the AIF for each patient is sometimes determined from the change in signal intensity as a function of time in a major artery (13). However, these measurements are subject to significant systematic errors due to the very high concentration of contrast agent in the artery and partial volume effects. It is rarely possible to find a large artery near the tumor, so that the AIF’s obtained from direct measurements in arteries are usually global rather than local. Other methods such as ‘Reference Tissue’ (14) approaches correct for the regional AIF. However, they do not produce a truly local AIF that reflect the delivery of the contrast media bolus directly to the tumor.
Ideally one would like to study contrast media uptake and washout from the tumor itself, without systematic errors or bias caused by the use of inappropriate physiologic models and incorrect AIFs. To find the hemodynamic characteristics of the tumor, it is necessary to separate the effects of the AIF from the kinetics of uptake and washout of contrast media. This can be done if the tissue is treated as a system that gives a linear, time invariant, and purely ‘causal’ response to the ‘local’ AIF (15). The local AIF can be deconvolved from the contrast concentration curves, so that the impulse response function (IRF) can be determined (16) as a fundamental characteristic of tissue.
Theoretically the deconvolution is easy to accomplish utilizing the properties of Laplace or Fourier transform. However, in practice the deconvolution is an ill-posed problem and noise in the data is magnified during processing of real DCEMRI data. The most commonly used deconvolution technique is singular value decomposition (SVD) (17). However, the conventional singular value decomposition introduces unwanted oscillations across a substantial portion of the IRF. Some regularization techniques have to be used as alternative singular value decomposition method to reducing the oscillations (18). These techniques are quite complicated and subject to error, and may not be appropriate for routine clinical DCEMRI data analysis. Therefore, a simple, robust deconvolution method that does not amplify noise is necessary to improve diagnosis of cancers.
Here, we describe a new deconvolution technique for determining the IRF. We test the method using numerical simulations of realistic DCEMRI data. Simulated plots of contrast media concentration as a function of time (referred to in the following as C(t)) for both tumor and muscle were generated using an EMM that has been demonstrated to accurately fit a wide range of experimental DCEMRI data (19,20). In additional to numerical simulations, the deconvolution technique was demonstrated directly with sample rodent tumor DCEMRI data acquired at 4.7 T. The TCM was also used to fit the contrast media concentration vs. time plots and the results of this procedure were compared with the results obtained with the deconvolution method. The physiological meaning of IRF parameters was determined by comparing them with Ktrans and ve derived from the TCM.
Under assumption that muscle is well approximated by the TCM, the AIF was derived from the simulated muscle C(t) using the ‘reference tissue method’ (14). To obtain the IRF, we used two steps, namely ‘prediction’, and ‘correction’ to deconvolute the AIF from C(t). In the prediction step, a numerical IRF was calculated using a recursive formula. Due to the rapid variations in the tissue C(t) and the AIF at very early times, there were oscillations at the early times (within the first one or two sampling periods) in the resulting numerical IRF. Therefore, in correction step, an equation with a limited number of parameters was used to fit the numerical IRF to eliminate these oscillations. These parameters were varied over a limited range so that the convolution between the AIF and the final IRF matched the original C(t).
The analysis described above produced IRF’s for tumors that generally included a rapid enhancement phase. This is not physically realistic, since in principle, we are determining the purely causal response to a delta function input of contrast agent. This ‘enhancement phase’ is most likely due to the fact that the AIF derived from a reference tissue reflects blood flow to a reference tissue near the tumor, i.e. a regional AIF. The contrast bolus is likely to experience additional delay and dispersion on the way to specific voxels or regions within the tumor, and this delay and dispersion produces an unrealistic IRF. Thus, we modified the AIF derived from the reference tissue by introducing additional delay and dispersion to produce a truly local AIF, and to obtain a monotonically decaying IRF.
For a causal system, the convolution integral gives the relationship between the input function of a linear, time-invariant system with the impulse response and output response functions. The contrast concentration curve vs. time (C(t)) can be considered as the convolution between AIF (Ca(t)) and the IRF (Cδ(t)), i.e., mathematically:
Here we assume that the AIF is zero for time t < 0 and the system has no initial contrast agent (zero initial conditions). The lower limit of zero is due to Cδ(t) = 0 for t < 0, and the upper limit of t is due to the system being causal, i.e. Ca(t) = 0 for t < 0. To obtain the Cδ(t) from above equation, a deconvolution has to be performed.
We use two steps namely ‘prediction’ and ‘correction’ to obtain the IRF from the Eq. . In the prediction step, we divide the whole time period into N very small time intervals in Eq. . Then at time tj simply applying Riemann sums using right-hand endpoints yields:
where we assume that Cδ(t) does not change much in the small time interval Δti = [ti−1, ti] and take the value at ti. This is a good approximation when C(t) and Ca(t) were calculated with high temporal resolution. Then a predicted IRF is obtained by solving Eq. , i.e.,
Obviously the Eq.  requires that Ca(0) ≠ 0. When Ca(0) = 0, we modify the Eq.  as follows. In Eq. , we take Ca(tj-τ) as the left endpoint instead of right. In this way, the Ca(0) in Eq.  was replaced by Ca(1) if Ca(0) = 0.
As previously mentioned, due to rapid variations in the Ca(t) and C(t) at very early times, there are oscillations at the early times in the numerical IRF. Please notice that the oscillations produced by the above calculation only appear during a short period of time (generally less than 10 seconds, which is less than one or two sampling periods for experimental MRI data). In contrast, the oscillations produced by ‘singular value decomposition’ occur throughout the entire sampling period (21). To eliminate these oscillations and obtain an accurate IRF, the following equation was used to fit the numerical IRF:
where K (min−1) is the transfer constant of the tissue, κi (i=1, 2) are the decay constants, λ (min−1) is uptake rate, ρ is related the slope of uptake, and ε is a scaling factor. Equation 4 can be related to the TCM with K = Ktrans, and κ1 = kep = Ktrans/ve, when λ 1 and ε ≈ 0. Please notice that the IRF includes an uptake phase here. In principle the response of tissue to a delta function input should be a monotonic decay. However, the results demonstrate that in tumors the IRF contains an uptake phase when a regional or global AIF is used rather than a truly local AIF (see below). In the correction step, we use the direct search technique to vary the parameters of the fitted IRF, so that the result of Ca(t) convoluted with Cδ(t) (AIFIRF) accurately fits the original C(t) with highest possible goodness-fit (R2) values. We refer to ‘AIFIRF’ as a ‘convolution fit’ in the remainder of this paper.
Kovar et al. developed a method to estimate the AIF from contrast media dynamics in a reference tissue with (approximately) known physiology (14). Under the assumption that the reference tissue, such as muscle, was well approximated by the two-compartment model, the AIF could be obtained from the muscle C(t) using the literature values of Ktrans and ve following:
where Ktrans(min−1) is the volume transfer constant between blood plasma and extravascular extracellular space (EES), and ve is the volume of EES per unit volume of tissue. For muscle, Ktrans = 0.11 min−1 and ve = 0.2 were used in the simulations (14). In practice Ktrans and ve can be varied within a reasonable range around the literature values to achieve an optimal fit to experimental data, but in the case of the present simulation study the literature values were used without modification.
Previously studies showed that the following empirical mathematical model (EMM) accurately fits the contrast media concentration vs. time curves generated from the DCEMRI data for a variety of different tissues including tumors (9):
where A is the upper limit of the tracer concentration, α (min−1) is the rate of uptake, β (min−1) is the overall rate of washout, γ (min−1) is the initial rate of washout, and q is related to the slope of uptake. The EMM was used to produce a range of simulated contrast media concentration vs. time curves. For numerical simulations, the temporal resolution of 0.6 s was used in the calculations of C(t) for muscle and tumor.
To select realistic EMM parameters in the simulations, the ranges of parameters were determined based on previous analysis of data from DCEMRI studies of rodent tumors (9). In that previous study, T1-weighted gradient echo images (TR/TE = 15/6 ms) with temporal resolution of ~5 s were acquired for 2 minutes before and ~60 min after I.V. contrast injection on a SIGNA 1.5T scanner. A total of 24 Copenhagen rats with transplanted prostate tumors (AT2.1 and AT3.1) in the hind limbs were studied. The dose normalized contrast agent concentration vs. time curves were calculated for region of interests (ROIs) of normal muscle, tumor rim, and tumor center. Those curves were accurately fitted with the EMM. The pooled data provided EMM parameters (average ± standard deviations) for normal muscle: A = 0.68 ± 0.11, q = 1.33 ± 0.33, α = 3.53 ± 1.35, β = 0.027 ± 0.009, γ = 0.17 ± 0.08; for tumor rim: A = 1.96 ± 0.53, q = 0.80 ± 0.32, α = 2.09 ± 1.39, β = 0.014 ± 0.008, γ = 0.04 ± 0.03; and for tumor center: A = 1.59 ± 0.96, q = 1.09 ± 0.39, α = 2.23 ± 1.90, β = 0.012 ± 0.009, γ = 0.06 ± 0.06. Therefore, the EMM parameters used in simulations (Table 1) for muscle and tumor were selected close to the average, and the average plus and minus one standard deviation. For each EMM parameter we tested values equal to the mean value plus and minus one standard deviation, while the other EMM parameters were kept at their mean value. This meant that a total of 10 different combinations of parameters were tested. No simulated noise was added to the EMM fits; the EMM fits came from experimental data and therefore the effect of noise is already incorporated.
Real DCEMRI data acquired at 4.7 T was used to test the deconvolution technique. T1-weighted gradient echo images (TR/TE = 40/3.5 ms, FOV = 4.0 cm, array size = 1282, slice thickness = 1 mm, number of slice = 3, NEX =1) of rodent transplanted prostate tumor were acquired with temporal resolution of ~5 s before and after 0.2 mmol/kg Gd-DTPA (Omniscan, GE Healthcare, Piscataway, NJ) injection for a total of ~10 min. After contrast media injection, the T1 in each image voxel was estimated by comparing signal intensity S(t) in a region of interest (ROI), to the control signal intensity S(0) in a reference tissue (for instance, muscle) with known T1. Since TR T1 for T1-weighted gradient echo imaging, the signal intensity could be approximated as a linear function of T1 and contrast media concentration (22), i.e.,
where R1 (4.3 mM−1s−1) is the longitudinal relaxivity of the contrast agent Gd-DTPA at 4.7 Tesla (23), Sref(0) is the signal at the ROI containing the reference tissue – muscle for this case, and T1(ref) = 1280 ms (24). To apply the deconvolution technique, the EMM was used to fit the C(t) obtained in Eq.  first, then high temporal resolution (0.6 s) C(t) was calculated from the EMM parameters.
The average values of EMM parameters shown in Table 1 were used to produce typical concentration versus time plots for muscle and tumor. Figure 1 (top) shows the contrast concentration curves for muscle (open circles) and tumor (open triangles); (middle) the corresponding AIF derived from the muscle curve; and (bottom) the IRFs obtained from deconvolution for muscle (gray line) and tumor (black line). The final results of convolution between AIF and IRFs (AIFIRF) are also shown in the top row for muscle (gray line) and tumor (black line). The AIFIRF fits the simulated data very well (R2 > 0.99). This demonstrates that the derived IRFs are accurate, under the assumption that the AIF calculated from the reference tissue is accurate. The IRF model parameters are given in Table 2 for the muscle and tumor. Thus, the parameter λ(=20.5) was about five times larger and ε (=0.003) was about 30 times smaller in muscle compared to tumor. In other words, the uptake component of the IRF for muscle was almost instantaneous and the amplitude of the second exponential decay was negligible. As expected, the IRF for muscle could be modeled as a single exponential decay, which was just the same as the two-compartment model with K = 0.11 = Ktrans and K/κ1 = 0.20 = ve, so that the IRF parameters for muscle provide no new information. However, in the tumor, a significant initial uptake component of the IRF is evident, and the washout phase is bi-exponential. For tumor, the IRF parameter K was almost the same as the corresponding TCM parameter Ktrans. However, K/κ1 = 0.63, which is much less than the corresponding TCM parameter ve= 0.81. As a result, the TCM often required unrealistically larger ve for tumor (see last two column in Table 2) in order to fit the C(t) accurately. If ve was forced into the range between 0.10 and 0.35, the fits to the simulated data were poor.
Figures 2 – 6 show the results of changing one of the EMM parameters - A, q, α, β, or γ –while keeping the other parameters at their mean values. In each case, AIFIRF was compared with the TCM fit. The detailed results of these simulations are given in Table 2.
Figure 2 shows the results for the parameter A set to 0.5 (mM) and 3.5 (mM). The contrast concentration vs. time plots were scaled when parameter A was changed and the shapes of IRFs were similar for these two cases but the value of K differed by a factor of seven. AIFIRF fitted the simulated data very well (R2 > 0.99), but the TCM had trouble fitting the curve for A = 3.5 (mM) with R2 = −1.56. Figure 3 shows the results for q equal to 0.5 and 1.5, resulting in small changes in the slope of the uptake portion of C(t). This resulted in a large change in the IRF parameter λ = 25.5 min−1 and 2.8 min−1 for q = 0.5 and 1.5, respectively. Therefore, the IRF for q = 0.5 could be modeled simply with a double exponential decay using four parameters, while the IRF for q = 1.5 required an uptake component. Again, the fit obtained with AIFIRF (R2 > 0.96) was better than the TCM fit (R2 > 0.66).
Figure 4 shows the results for the parameter α set to 0.25 min−1 and 3.75 min−1. For slow contrast uptake (α = 0.25 min−1), the IRF had a strong uptake component and double exponential decay. On the other hand, for fast contrast uptake, the IRF was accurately modeled simply as double exponential decay due to the large λ (= 52.1 min−1) value. The TCM fit was more accurate for slow contrast uptake than for fast uptake (R2 = 0.93 vs. 0.78), but overall the AIFIRF fits (R2 > 0.98) were better than the TCM fits. Figure 5 shows the results when the parameter β changed to −0.005 and 0.025 min−1. The IRFs obtained for these two values of β were similar at early time points, but very different at later time points (where curve shape was dictated primarily by the second, slower component of the exponential decay) because of large differences in the values of ε and κ2. The TCM fitted the curve better for the β = 0.025 (R2 = 0.99) than β = −0.005 min−1 (R2 = 0.95), but even for the large β value, fits obtained for the TCM were not as good as those obtained with convolution fits (R2 > 0.99). For the parameter γ, there were a similar trends, see Fig. 6. This is because both parameters control the contrast media washout.
In above results, we changed EMM parameters one at a time from their average values. We also changed all five EMM parameters simultaneously (Table 2) to generate two more examples: (i) slow uptake and slow washout; and (ii) fast uptake and fast washout. The results of deconvolution between concentration curves and the AIF are shown in Fig. 7. The IRF for case (i) requires a rapid uptake component and a double exponential decay; and for case (ii) only a double exponential decay is required. Once again, the convolution fits (R2 > 0.96) of the contrast concentration curves were better than the TCM fits (R2 > 0.45).
Figure 8 (top dash line) shows the IRF deconvoluted from a tumor contrast concentration curve using a regional AIF (Fig. 8 bottom dash line). The IRF contains a rapid enhancement phase with significant amplitude. To minimize the rapid enhancement phase, the regional AIF was modified with increased delay and dispersion. An AIF with increased delay of approximately 7 s and width-at-half-height increased by a factor of two (shown in Fig. 8 bottom, gray line), produced an IRF with a minimal uptake phase, but with bi-exponential decay (shown in Fig. 8 top gray line).
Figure 9 illustrates the application of this approach to raw data from a non-metastatic rodent prostate tumor. The Figure shows dose normalized Gd-DTPA C(t) curves (open circles) for ROIs in muscle (30 pixels) and tumor rim (35 and 33 pixels for ROI1 and ROI2, respectively) fitted with the TCM (red line), and convolution (blue line). The convolution (AIFIRF) fits the data more accurately than the TCM. The calculated goodness-fit (R2) values for the three ROIs were 0.98, 0.99, and 0.99 for convolution fit, and 0.95, 0.98 and 0.45 for the TCM. The corresponding IRFs (black line – calculated directly from Eq.  and dash line - after correction) from muscle and two different tumor rim ROIs are shown in Fig. 10. As in previous examples, the tumor IRFs include a large uptake component, and the two-compartment model does not provide good fits to the raw data from tumor ROI2.
The results presented here demonstrate a simple and effective algorithm for deconvoluting the global AIF from DCEMRI data and evaluating the resulting IRF. The deconvolution technique was tested on simulated contrast concentration vs. time curves. The simulations were based on EMM fits to a large number of experimental datasets, and thus accurately represent the range of data that would be seen in typical experiments. The results demonstrate that the deconvolution technique produces excellent fits to DCEMRI data without artifacts and is straight forward and easy to use. This is a new approach to analysis of DCEMRI data that could lead to improvements in diagnostic accuracy and assessment of response to therapy.
Because muscle is well approximated by the simple two compartment model, the IRF’s for muscle were accurately fit with a single exponential decay equivalent to the two compartment model, with K = Ktrans, and κ1 = kep = Ktrans/ve, when λ 1 and ε ≈ 0. No initial uptake phase was required for muscle impulse response functions.
However, the IRF required for tumors was much more complicated than for normal tissue; this reflects the complex physiology and micro-anatomy of tumors. When a ‘regional’ AIF derived from a nearby muscle reference tissue was used, the resulting tumor IRF included an initial exponential enhancement phase. This is physically unrealistic, since the response of the tumor to a delta function bolus should decay monotonically. This suggests that the ‘regional’ AIF derived from reference tissue must be modified by adding delay and dispersion terms to produce a local AIF for each tumor region. The results illustrated in Fig. 8 demonstrate that the ‘local’ AIF gives a physically realistic IRF that is a (approximately) bi-exponential decay. The additional delay and dispersion terms that modify the ‘regional’ AIF to produce a ‘local’ AIF may be related to tumor interstitial pressure, and blockage and intermittent flow in tumor microvessels. As a result, these are novel parameters that may have diagnostic utility and may be sensitive to changes in tumors during therapy. However, the potential clinical application of these new parameters is speculative at present, and further work is needed to determine whether they can be useful in routine clinical practice.
In this study, the deconvolution to obtain an IRF assumed that an accurate AIF was used. The AIF was derived from reference tissue – muscle, and did not include a second pass. Although improvement in accuracy of the AIF was not the subject in the present study, future work will evaluate the effect of an AIF with a second pass. In addition, we will test AIF’s calculated using the more accurate ‘double reference tissue method’ developed by Yang et al (25), and AIF’s derived from direct arterial measurements.
The widely used two-compartment model approach requires that the impulse response function is simply a mono-exponential decay. The present results demonstrate that this is the case for muscle. However, this simple model is not consistent with the IRF calculated for most tumor ROIs. The tumor IRF usually included an initial uptake phase followed by a multi-exponential decay. These characteristics of the tumor IRF produce systematic errors when the two compartment model is used to fit raw data. The approach presented here provides solutions to these problems. First, the initial IRF and AIF can be adjusted to produce a local AIF and a more realistic, monotonically decaying IRF. Second, the monotonically decaying IRF can be well approximated by a multi-exponential decay. In general, given the signal-to-noise ratio of typical data, a bi-exponential decay is sufficient. Thus, the present approach leads to calculations of IRF’s that accurately model experimental data. Biases due to use of the two compartment model or other physiologic models are avoided. Although the IRF approach is a ‘model free approach’, it can be analyzed to approximate physiological parameters, as illustrated here.
The present work demonstrates that impulse response analysis has potential advantages compared to more widely used methods. We describe a simple, robust deconvolution method that produces artifact-free IRF’s. An extension of the method allows calculation of a local AIF and provides new parameters that describe delay and dispersion of the contrast media bolus in tumor vasculature, and these parameters may have diagnostic utility. The results of this preliminary study suggest that the parameters of the model for the IRF can provide good discrimination between cancer and normal tissue (see Table 2). Future work will further evaluate the use of this approach to distinguish between benign and malignant cancers, and detect tumor response to therapy.
This work was supported by grants from the Lynn S. Florsheim Foundation, the NIBIB (RO1 EB003108-01), the NCI (RO1CA78803), and the Specialized Program of Research Excellence (SPORE) (1 P50 CA125183 02).