Influence of contrast injection - Wiley Online Library

28 downloads 11783 Views 1MB Size Report
As both the TDM and PKM represented low-frequen- ... administered in a single injection, where high rates only mod- .... the frequency domain, respectively.
Magnetic Resonance in Medicine 59:1111–1119 (2008)

System Identification Theory in Pharmacokinetic Modeling of Dynamic Contrast-Enhanced MRI: Influence of Contrast Injection H.J.W.L. Aerts,1,2 N.A.W. van Riel,1,3 and W.H. Backes4* Optimization of experimental settings of dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI), like the contrast administration protocol, is of great importance for reliable quantification of the microcirculatory properties, such as the volume transfer-constant Ktrans. Using system identification theory and computer simulations, the confounding effects of volume, rate and multiplicity of a contrast injection on the reliability of Ktrans estimation was assessed. A new tracer-distribution model (TDM), based on in vivo data from rectal cancer patients, served to describe the relationship between the contrast agent injection and the blood time-course. A pharmacokinetic model (PKM) was used to describe the relation between the blood and tumor tissue time-courses. By means of TDM and PKM in series, the tissue-transfer function of the PKM was analyzed. As both the TDM and PKM represented low-frequency-pass filters, the energy-density at low frequencies of the blood and tissue time-courses was larger than at high frequencies. The simulations, based on measurements in humans, predict that the Ktrans is most reliable with a high injection volume administered in a single injection, where high rates only modestly improve Ktrans. Magn Reson Med 59:1111–1119, 2008. © 2008 Wiley-Liss, Inc. Key words: injection shape; pharmacokinetic model; system identification; DCE-MRI; arterial input function; AIF; tracer distribution model; TDM; tumor perfusion

The development of dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) has led to the noninvasive characterization of tumor tissue (1,2). The volume transfer constant Ktrans, which describes the leakage rate of contrast agents into the tumor tissue, is a pharmacokinetic parameter that provides information on the physiological properties of the tumor, including the microvascular permeability and surface area, the relative blood volume, and the relative volume of the interstitial space (3– 8). DCE-MRI has been recognized as a valuable tool for the evaluation of novel antiangiogenesis therapies (9 –12). Clinically, it is of importance to have a reliable and accurate quantitatively assessment of the progression or regression of tumor angiogenesis, as it may affect continuation of or changes to

1Control Systems Group, Department of Electrical Engineering, Eindhoven University of Technology, Eindhoven, the Netherlands. 2Department of Radiation Oncology (MAASTRO), Research Institute Growth and Development, University of Maastricht, Maastricht, the Netherlands. 3Biomodeling and Bioinformatics Group, Department of Biomedical Engineering, Eindhoven University of Technology, Eindhoven, the Netherlands. 4Department of Radiology, Research Institute Growth and Development, Maastricht University Hospital, Maastricht, the Netherlands. *Correspondence to: Walter H. Backes, Department of Radiology, Maastricht University Hospital, P. Debyelaan 25, 6229HX, Maastricht, the Netherlands. E-mail: [email protected] Received 17 August 2007; revised 7 January 2008; accepted 8 January 2008. DOI 10.1002/mrm.21575 Published online in Wiley InterScience (www.interscience.wiley.com).

© 2008 Wiley-Liss, Inc.

the treatment plan (13–15). The quality of the assessment is reflected in the accuracy of the determined pharmacokinetic parameters quantifying the wash-in and wash-out of the contrast agent. This accuracy is to a large extent determined by the experimental conditions (16). Unfortunately, where efforts have been directed to the standardization of the pharmacokinetic modeling (17), little is known on the most optimal experimental conditions. Previous investigations have focused on influences of sampling time, signal-to-noise ratio (SNR), and experiment duration (18 –22). In the present study it is hypothesized that the accuracy of the estimated pharmacokinetic parameters strongly depends on the dynamic shape of the arterial input function (AIF), which can be controlled by the contrast administration protocol. System identification theory is used to relate the reliability of the estimated Ktrans to the dynamics of the contrast agent injection and to find the optimal protocol for contrast agent administration. This theory states that the number of parameters involved in a model, here the pharmacokinetic model (PKM), is a compromise between complexity of the data and identifiability of the parameters (23). The tradeoff between complexity and identifiability is referred to as “minimal modeling.” The most commonly applied PKM is the two-compartment minimal model (4,5,17), which uses kinetic parameters to describe the exchange of contrast agent between the blood plasma compartment and extravascular-extracellular space (EES). Current approaches for calculating these kinetic parameters from DCE-MRI data show large uncertainties (16,21). Explanations can be found not only in the experimental conditions but also in the intrinsic mathematical structure of the model (22). To find the optimal contrast administration protocol, the relation between the shape of the injection and the blood concentration time-course was determined. Here a new tracer distribution model (TDM) was developed that is based on linear system theory and uses a plug (trapezium) shape of the contrast injection as input and predicts the arterial concentration (i.e., AIF) as output. Using the TDM and the PKM in series, the relation between injection protocol, arterial concentration, and tissue concentration were analyzed in both the time and frequency domain. Parseval’s theorem on the energy density and frequency analysis were used to compare the resulting Ktrans values in response to varying contrast injection protocol. By means of theoretical analysis and computer simulations, the injection protocol was varied in terms of injection volume, injection rate, and multiplicity (i.e., single vs. multiple injections). The results provide information on which injection protocol is optimal for accurate and reliable phar-

1111

1112

Aerts et al.

macokinetic parameter estimation. Hence the aim of the study is to determine the best protocol for contrast agent injection in a clinical setting. THEORY PKM The distribution of a contrast agent between the blood plasma and the tumor tissue was according to Tofts et al (4,5,17) standardized to:



t

C t共t兲 ⫽ v pC p共t兲 ⫹ K trans

C p共u兲e 共⫺K

trans/v 兲共t⫺u兲 e

du,

[1]

0

where Ct(t) and Cp(t) are the concentrations (in mol/liter) in the tissue and blood plasma, respectively. Ktrans is the volume transfer constant (in min–1), ve is the fractional volume of the EES (unitless), and vp is the blood plasma volume per unit volume of tissue (unitless). Ktrans, ve, and vp represent the kinetic parameters that have to be determined from DCE-MRI measurements and the subsequent PKM analysis. The parameter Ktrans depends on the microvascular blood flow through the tissue, the permeability for the contrast agent, and the surface area of the tumor microvasculature. Here the focus will be on Ktrans, as it is generally considered the most important parameter that characterizes the malignancy of the tumor tissue. In the frequency domain the PKM reads: C t共␻兲 ⫽ H PKM共␻兲 䡠 C p共␻兲,

[2]

where ␻ is the angular frequency (in rad/s), and Cp(␻) and Ct(␻) are the input and output concentration functions in the frequency domain, respectively. The tissue-transfer function HPKM(␻) equals: H PKM共␻兲 ⫽

v pi␻ ⫹ K trans共v p/v e ⫹ 1兲 i␻ ⫹ K trans/v e ⫽ vp ⫹

K trans . i␻ ⫹ K trans/v e

[3]

Influences of parameter value changes on HPKM are shown in Fig. 1a. Here the modulus M(f) ⫽ 20log(|HPKM|) (in dB) is displayed as a function of the frequency f (in Hz) on a logarithmic scale to obtain a so-called Bode diagram. Ktrans The frequency fp ⫽ (in Hz) is the pole and fz 2␲ 䡠 ve Ktrans ⫽ 共1/ve ⫹ 1/vp 兲 (in Hz) represents the zero of HPKM. 2␲ The transfer at lower frequencies is completely determined by HPKM 共0兲 ⫽ ve ⫹ vp , which represents the entire tissue space available for contrast material distribution. For higher frequencies, the transfer is determined HPKM 共⬁兲 ⫽ vp , which represents the direct feed-through (i.e., no exchange of contrast material) of the microvascular system. Tumor tissue with “leaky” microvasculature is characterized by relatively high Ktrans values. The corresponding HPKM function will display relatively strong high-frequency components, which represents the rela-

FIG. 1. a: The spectrum of the tissue-transfer function HPKM for Ktrans ⫽ 0.1 and 0.5 min–1. The other pharmacokinetic parameters, ve ⫽ 0.3 and vp ⫽ 0.08, were kept constant. The modulus of the system-transfer function (i.e., 20log(HPKM)) is plotted vs. the frequency (in Hz). The vertical lines represent the frequency of the pole (fp ⫽ 0.88 mHz) and zero (fz ⫽ 4.2 mHz) of HPKM for Ktrans ⫽ 0.1 min–1. The horizontal lines indicate the tissue transfer at f ⫽ 0 and for asymptotically high frequency values. Note that for increasing Ktrans the spectrum contains more high-frequency components. b: The tissue concentration time-course for the two Ktrans values with equal Cp is depicted. Note that when Ktrans increases, the initial tissue wash-in rate and the maximal enhancement also increase.

tively fast enhancement of the tumor tissue. Moreover, for higher Ktrans values the output Ct(t) (with the same Cp(t)) will increase faster and reach higher peak levels (see Fig. 1b). Optimal Excitation For optimal parameter identification, all relevant frequencies of the system should be observable in the output Ct(t) (23). This identifiability depends on the excitation, the distribution of energy in the input signal over the frequencies (i.e., the energy density |Cp(␻)|2), and the frequency characteristics of the system HPKM(␻), which is determined by the kinetic parameters. When the frequency range covered by the input matches the frequency characteristics of

System Identification in DCE-MRI

1113

the system, all relevant frequencies for identification are transferred to the output and the model is considered to be optimally excited. For the pharmacokinetic system optimal excitation is problematic because the input Cp cannot be influenced directly and the various frequencies cannot be controlled independently. The actual independent input is the injection of contrast agent into the vascular system, illustrated in the model setup of Fig. 2. When the bolus of contrast material is injected, it disperses through the vascular system, which decreases the modulus more strongly at high frequencies than at low frequencies. For this reason Cp can only be influenced indirectly by the injection protocol. Because in DCE-MRI only sampled versions of the true contrast agent concentration time-courses can be measured, the signals are frequency-band limited. The Nyquist theorem states that the maximum frequency fmax of the sampled versions equals fmax ⫽ fs/2, where fs is the sample frequency. If the true signals Cp and Ct contain frequencies above fmax, these contribute to, and therefore contaminate (i.e., alias), the sampled signal values of the frequencies below fmax. Energy (Parseval) Parseval’s theorem states that the total energy (Ex) of a time signal x(t) across all times t can be calculated, both in the time domain and in the frequency domain, by integrating the squared modulus spectrum (24):

Ex Ⳏ





兩x共t兲兩 2dt ⫽

⫺⬁





兩X共f兲兩 2df,

[4]

⫺⬁

where the signal x(t) can be any time signal (e.g., Cp(t) or Ct(t)). The theory now states that the more energy an injection protocol contains, the higher the overall moduli of the frequencies, and the better the excitation. For reliable pharmacokinetic parameter determination it is desired that the distribution of the energy, i.e., the energy density |X(f)|2, should be matched to the frequency transfer of the model, such that all relevant frequencies of the input are significantly transferred to the output.

where the TDM parameters Ai and Bi were determined from experimental data. The TDM amounts to a low-pass linear filter describing energy dissipation and can predict the arterial concentration for a specific injection shape: C p共␻兲 ⫽ H TDM共␻兲 䡠 l共␻兲.

[6]

METHODS Parameter Estimation The standard least squares (LSQ) method (23) was used to estimate Ktrans and the TDM parameters. The AIF was determined by DCE-MRI measurements in the femoral arteries in humans, as published by De Lussanet et al. (26), in which the signal intensity curves were converted to concentration curves by using the variable flip angle method for fast gradient echo pulse sequences (27). The injection comprised 22.5 ml of 0.5 mol/liter gadolinium diethylenetriamine pentaacetic acid (Gd-DTPA) (Schering, Berlin, Germany) administrated at a rate of 3 ml/s in the antecubital vein, with the use of an automated injector. Parameter estimations of the PKM and TDM was implemented in MATLAB (The Mathworks Inc., Natick, MA, USA), where the function lsqnonlin was used as leastsquares algorithm with a Levenberg-Marquardt optimization algorithm. The boundaries set to the pharmacokinetic parameters were: 0.001 ⬍ Ktrans ⬍ 1.5 min–1, 0.001 ⬍ ve ⬍ 0.5, and 0.001 ⬍ vp ⬍ 0.5. The initial values of the parameters were randomly selected from these ranges. The relative and absolute tolerance of the algorithm were both set to 10– 6. A trapezium-shaped function I(t) was used to represent the injection rate as a function of time. The rates of I(t) were set to ␮l/s (i.e., a multiplication of 1000). The estimation of the TDM parameters resulted in: B3 ⫽ 1.0 ⫻ 10– 6, B2 ⫽ 2.3 ⫻ 10–5, B1 ⫽ 1.0 ⫻ 10– 4, B0 ⫽ 3.4 ⫻ 10– 6, A2 ⫽ 0.25, A1 ⫽ 4.7 ⫻ 10–2, and A0 ⫽ 6.9 ⫻ 10–5. Figure 3 shows the identified TDM and the measured data points. The identified TDM parameters were based on one combination of injection volume and rate to obtain AIFs for other combinations of injection volume and rate, as used in the curves shown in Fig. 2 and the simulations. Simulations

TDM To obtain a quantitative relation between the dynamic shape of the contrast injection I(t) and the blood plasma concentration Cp(t), a new TDM was developed. Here I(t) is represented by the injection rate (in ␮l/s) as a function of time. Stollberger et al. (25) showed that linear system theory is well applicable for the distribution of contrast material through the vascular system. The circulation is described as a “black box” and the model was only selected on how well it describes the relation between the input and output data. This relation was described using a third-order linear transfer function in the frequency domain: H TDM共␻兲 ⫽

B 3␻ 3 ⫹ B 2␻ 2 ⫹ B 1␻ ⫹ B 0 , ␻ 3 ⫹ A 2␻ 2 ⫹ A 1␻ ⫹ A 0

[5]

To investigate the influence of different experimental variables, Monte Carlo simulations were performed. Using the PKM and TDM (Eqs. [3] and [5]) the response of the arterial and tissue concentration can be calculated (Eqs. [2] and [6). By changing the injection protocol, the (i.e., true) Ktrans value (range ⫽ 0.05– 0.5 min–1) (e.g., Refs. 22 and 28), signal-to-noise ratio (SNR), and sample time, the simulated arterial and tissue concentrations were calculated. Gaussian noise with zero mean was added to the concentration course Ct(t) to mimic measurement noise. The SNR was defined as the ratio of the power of the signal Ct divided by the power of the noise, both averaged over time for one injection protocol (injection volume ⫽ 45 ml, injection rate ⫽ 5 ml/s) and one set of pharmacokinetic parameter values (Ktrans ⫽ 0.5 min–1, ve ⫽ 0.3, and vp ⫽ 0.08). The absolute noise level was rescaled from this set of

1114

Aerts et al.

FIG. 2. Schematic representation of the model setup. The injection profiles I are transferred to blood concentrations Cp by the TDM. The resulting blood plasma concentration Cp serves as input for the tumor PKM which has the tissue concentration course Ct as output. The TDM, the PKM, and the energy densities of subsequent input and output functions are also displayed in the frequency domain. The low-pass filtering of the vascular system and tumor tissue are shown in the Bode diagram of the TDM and PKM. An injection profile with an injection rate of 2 ml/s and a volume of 20 ml is shown (shape A, in black). To show the influence of the volume, an injection with an increased volume of 40 ml and the same rate is depicted (shape B, in dashed blue). The influence of the rate is shown with (shape C, in dash-dotted red), a rate of 4 ml/s and the same volume as shape A (20 ml). Note that there is a larger influence of the volume than of the rate on the arterial and tissue concentration curves in both the time and frequency domain. Shapes A and C overlap in the energy plots of Cp and Ct.

System Identification in DCE-MRI

1115

15 min for all simulations. The standard settings of the sample frequency and SNR were 0.1 Hz and 10, respectively. To investigate the influence of the contrast injection on the Ktrans estimates, the injection volume and rate were varied between 5 and 45 ml and 0.5 and 5 ml/s, respectively, which are clinically relevant. Here the true Ktrans parameter was varied with constant experimental conditions. To assess the effect of sample frequency, fs was varied from 10–3 to 102 Hz. To investigate the influence of the measurement noise, the SNR was varied between 1 and 100. Finally, to assess the effect of multiple injections, single, double, and triple injections were analyzed and compared as a function of the injection volume. Also the number of injections (up to 20) was further increased, and the timing between the injections was varied (range ⫽ 0.5– 60 s).

FIG. 3. The arterial concentration time-course of a DCE-MRI patient study, injected with a volume of 22.5 ml 0.5 M gadolinium-chelate and rate of 3 ml/s. The circles represent the measurement points and the black curve the fitted result of the TDM.

injection values and pharmacokinetic parameter values, to assure the same absolute noise level for other injection protocols and to comply to the MR SNR characteristics. For varying sampling times, the SNR was rescaled by the square root of the sample time (in seconds), to correspond with MRI measurements. From these data sets, for which the true Ktrans values were known, Ktrans values were predicted for the simulated variations in contrast agent injection. The other kinetic parameters were kept constant at ve ⫽ 0.3 and vp ⫽ 0.08 within their physiologic range (22,28). The duration of the acquisition was kept constant at

RESULTS Model Setup Figure 2 shows for varying injection volumes and rates the time courses of the injection profiles I(t), the arterial concentration time courses Cp(t) resulting from the TDM, and the tissue concentration time courses Ct(t) resulting from the PKM. A trapezoid shape was assumed to mimic the contrast injection of a power injector. The integral of the injection profiles I(t) were equal to the injected volumes. Figure 2 also shows the energy densities in the frequency domain of the time courses. Consider the two injection shapes with a volume of 20 ml (shape A) and 40 ml (shape B), both at a rate of 2 ml/s. The total energy of injection shape A is 39 (units ml2/s) and for shape B it is 79 (see Eq. [4]). This increase in energy due to the volume results in a strong increase of the low-frequency contents

FIG. 4. Multiple injections in the time (a) and frequency domain (b), with a constant injection rate of 2 ml/s and total injection volume of 20 ml. Single (black solid), double (black dashed), and triple (gray dashed) injections are shown. The amplitude of the injection profile in the frequency domain was squared to show the energy density spectrum. Although multiple injections excite higher frequency components (f ⱖ 0.1 Hz), the loss in energy density at low frequencies (f ⬍ 0.1 Hz) is much stronger.

1116

Aerts et al.

FIG. 5. Simulation results showing the influence of the injection volume and rate on the Ktrans estimates. The data points represent the results for a true value of Ktrans of 0.5 (black dots), 0.2 (gray dots), and 0.05 (black crosses) ml/min/cm3. The black lines indicate the ideal case, in which the estimated Ktrans values are equal to the true Ktrans. a: Varying volume with a constant rate of 1 ml/s. b: Varying volume with a constant rate of 5 ml/s. c: Varying rate with a constant volume of 10 ml. d: Varying rate with a constant volume of 45 ml. The other kinetic parameters were kept constant at ve ⫽ 0.3 and vp ⫽ 0.08. The sample frequency and SNR were constant at 0.1 Hz and 10, respectively. A larger volume and rate always resulted in more precise Ktrans values. However, the influence of the volume is larger than the rate. Due to effects of the sample frequency, a bias for the calculated Ktrans is visible. For example, an injection with a rate of 5 ml and volume of 30 ml, with a true Ktrans of 0.5 min–1, will result in an accuracy (bias) of 16% and precision (scatter) of 40%.

of the energy spectra of the injection profile, the arterial concentration, and the tissue concentration (Fig. 2). When the injection volume is kept constant and the injection rate is varied, the total energy of shape C (in Fig. 2) becomes 75 for a volume of 20 ml and a rate of 4 ml/s. Although for faster injection the total energy is higher compared to shape A, and the distribution of this energy in the spectrum of the injection profile contains considerably more high-frequency components, the influence on the arterial and tissue concentration spectra is small. This is due to the low-frequency-pass characteristics of both the TDM and PKM functions. As a result, the optimal energy distribution, i.e., the energy density, of the input strongly increases with the volume and only slightly with the rate. With multiple injections the total injection volume is divided in multiple different fractions injected at separated time points (Fig. 4). Due to the quadratic dependence of the energy on the amplitude, multiple injections are less energetic than a single injection with the same total volume and rate (Fig. 4a). Energy is lost due to the extra rising and falling slopes in the multiple injection scheme. Double and triple injections show modestly higher frequency components for frequencies above 0.1 Hz, but foremost show strongly decreased energy densities at frequencies

below 0.1 Hz (Fig. 4b). The small increase at the higher frequency components for multiple injections has no strong effect on Cp due to the low-pass filtering of the TDM. Ktrans Simulations Figure 5 shows the estimated Ktrans for varying injection volumes and rates. For high values of the true Ktrans (⬎0.2 min–1), the injection protocol becomes increasingly important and higher injection volumes and rates are beneficial. For low values of Ktrans (⬍0.2 min–1), the errors are much smaller and injection volumes and rates are less critical. A larger volume will result in more precise estimates for both low (Fig. 5a) and high rates (Fig. 5b). When the volume is low (⬎15 ml), all estimates show a large variation independent of the rate (Fig. 5c). For a high volume, the influence of the rate is larger (Fig. 5d). A low rate (ⱕ2 ml/s) will result in less precise results. Under all conditions a high true Ktrans value can be estimated less precise, where the highest volume and rate provide the highest precision on the calculated Ktrans. Using two different injection profiles, the effect of the sample frequency has been investigated (Fig. 6a). Due to the rescaling of the SNR by the sample frequency, corre-

System Identification in DCE-MRI

1117

FIG. 6. Simulation results showing the influence of the sample frequency (a) and SNR (b) on the Ktrans estimates using two different injection profiles: 1) a volume of 45 ml and a rate of 5 ml/s (black dots); and 2) a volume of 10 ml and rate of 1 ml/s (gray crosses). The ve and vp values were set to 0.3 and 0.08, respectively, with a constant SNR of 10 (a) and sample frequency of 0.1 Hz (b). The black line indicates the ideal case, in which the estimated Ktrans values are equal to the true Ktrans. Note that the injection profile with the high volume and rate resulted in the best Ktrans values under all sample frequency and noise conditions. No accurate or precise estimates are possible for frequencies below approximately 0.07 Hz, independent of the injection profile. The precision of the estimates improves for a higher SNR, although the inaccuracy (i.e., bias) remains constant.

sponding with MRI signal characteristics, the noise varied by the square root of the sample frequency. The injection profile with a high injection volume (45 ml) and high rate (5 ml/s) resulted in more precise Ktrans estimates for all sample frequencies compared to the profile with a low volume (10 ml) and rate (1 ml/s). For a true Ktrans of 0.5 min–1, a sample frequency of 0.1 Hz, used for the other simulations, will result in a lower accuracy due to a bias of the Ktrans estimates. Also the SNR appeared important for the precision, as high SNR values (thus with low noise) resulted in more precise estimates (see Fig. 6b). Under all conditions the injection profile with the highest volume and rate, yielded the best results and was less sensitive to noise. Also for this profile the precision of the estimates improves for a higher SNR. However, the accuracy remained the same for all noise conditions due to the sample time.

The effects of a single, double, and triple injection on the Ktrans estimates are shown in Fig. 7. For the whole range of volumes, the single injection had the most precise estimates compared to double and triple injections. The triple injection even resulted in poorer estimates than the double injection. Further increasing the number of injections (up to 20) and varying the timing between the injections did not improve the Ktrans estimates (results not shown). DISCUSSION The aim of this study was to determine the best protocol for contrast agent injection for accurate and precise pharmacokinetic parameter estimation. DCE-MRI combined with pharmacokinetic modeling is increasingly applied to determine the volume transfer constant Ktrans, representing the microvascular leakage rate in tumor tissue. To opti-

FIG. 7. Simulation results showing the influence of the number of injections on the Ktrans estimations. a: Single injection (black dots) compared to double injection (gray crosses). b: Single injection (black dots) compared to triple injection (gray crosses). The injection rate is set to 3 ml/s, and the kinetic parameters Ktrans, ve, and vp are set to 0.5 min–1, 0.3, and 0.08, respectively. The sample frequency and SNR were constant at 0.1 Hz and 10, respectively. The black line indicates the ideal case, in which the calculated Ktrans equals the true Ktrans. A single injection appears favorable above double or triple injections. Note that a bias is visible due to the sample frequency. However, for all injections the precision of the estimates improves for increasing volume.

1118

mize the reliability of the determined Ktrans, the experimental settings of the contrast material injection could be optimized. In this study we showed by theoretical and simulation analysis that the contrast injection volume, rate, and multiplicity may largely influence the predictability and reliability of Ktrans. System identification theory was used to investigate the reliability of the determined Ktrans parameter. To this end, the subsequent input and output functions of the macrocirculation and tumor tissue microcirculation were described in both the time and frequency domain. The distribution of the contrast agent in the blood plasma space as well as the transfer of the contrast agent from the blood plasma space to the tissue appeared to be low-frequencypass systems. This implies that high frequency inputs are strongly attenuated to the output, i.e., the tissue concentration time-course. Increasing the contrast injection rate, which excites higher frequencies, only modestly improved the reliability of the Ktrans determination. Distributing the injection volume over multiple injections instead of one injection increased the high frequency components at the cost of decreasing the lower frequency components and resulted in poorer parameter estimates. Low frequency input, on the contrary, was strongly transferred by both the blood distribution and the extravasation effect. Since the low frequency energy of the input can be strongly controlled by the injection volume of the contrast agent, increasing the volume significantly improved the Ktrans determination. The relation between the injection profile and the accuracy is complex, strongly depending on a number of factors including pharmacokinetic parameter values, SNR, and sample frequencies. For the interpretation of the results, trends and relations were emphasized rather than exact values. Injection volume and rate, had a strong influence on the estimation of the Ktrans parameter. Low injection rates (below 2 ml/s) and low injection volumes (below 15 ml of Gd-DTPA) lead to unreliable Ktrans values and should therefore be avoided (Fig. 5). Henderson et al. (18) also reported that the most accurate Ktrans was found for the highest injection rate. In addition to their study, we evaluated the effect of the actual shape and, perhaps even more importantly, the volume of the injection. Our results show a low accuracy for the Ktrans estimates in the simulation study (Fig. 5–7). This reduction in accuracy was due to the combination of the low sample frequency (0.1 Hz) and high true Ktrans value (0.5 min–1). However, the accuracy improves for higher sample frequencies (Fig. 5), or lower true Ktrans values (Fig. 6a) (22). Injecting the volume in multiple fractions separated over time had a negative effect on the energy content, in agreement with theory, because of loss of energy due to the extra rising and falling slopes. This was supported by both the frequency and simulation analysis. A tradeoff was observed between a relatively strong influence decrease of in the lower frequency and a small increase of the higher frequency components. As the lower frequency range is most essential for accurate identification of the parameters, due to the low-pass filtering of both the TDM and PKM, multiple injections are not recommended. The selection of the TDM was based on AIF data from a DCE-MRI patient study, which resembled the AIF concen-

Aerts et al.

tration curves reported in literature (26). This patient study was performed with the same injection volume and rate as used in the present simulation study. A limitation of the present study is that for other injection protocols, the arterial curve is a prediction of the TDM. Further studies may validate the TDM predictions. However, such an experimental validation is complicated because various different injection protocols, including multiple injections, have to be conducted for every participating patient under identical conditions. Such a validation study would require a difficult clinical setup. Numerous confounding factors can disturb the DCE-MRI measurements, such as the measurement location of both the AIF and the tumor signal and variations in movements and cardiac output of the patient (29). Using a simulation study only the factors under study can be analyzed, and the influence of the confounding factors is absent. Another limitation of this study is the fact that no actual MR-derived signal time-courses were used, but that concentration time-courses were computationally generated without actual measurements. This approach is correct at low concentrations as MRI signal changes relative to the nonenhanced tissue are proportional to the concentration of the contrast medium. However, for high injection rates and/or volumes the signal increase becomes less than would be expected from extrapolations of the low concentration regime (e.g., due to T*2 effects), which effects the Ktrans determination. This study demonstrated the advantage of using high doses (up to 45 ml of 0.5 mol/liter Gd-DTPA) with respect to the accuracy and precision of pharmacokinetic parameter determination. However, it should be mentioned that currently the administration of gadolinium-based contrast media is discouraged in patients with severe renal impairment (i.e., estimated glomerular filtration rate ⬍30 ml/ min) and particularly the use of high doses (30).

CONCLUSIONS The injection volume and rate influence the accuracy and precision of the pharmacokinetic leakage parameter Ktrans. Low injection volumes (below 15 ml of 0.5 M gadoliniumbased contrast material) and small injection rates (below 2 ml/s) cause less precise Ktrans values. Multiple injections are also suboptimal. For practical application, single, fast, and, foremost, high-volume contrast injection is advised to obtain the most reliable Ktrans estimation. For optimal experimental conditions, measurements of the Ktrans parameter should reckon with the low-frequency-pass characteristics of both the contrast agent circulation and tissue leakage characteristics. This implies that the sample frequency should be sufficiently high to reveal the frequency characteristics of the tumor tissue and that the frequency spectrum of the injection profile should match the frequency spectrum of the tissue-transfer function.

ACKNOWLEDGMENT We thank Prof. Dr. Ir. P.P.J. van den Bosch, who provided many useful insights on this topic.

System Identification in DCE-MRI

REFERENCES 1. Jackson A, Buckley DL, Parker GJM. Dynamic contrast-enhanced magnetic resonance imaging in oncology. Berlin: Springer-Verlag; 2005. 311 p. 2. McDonald DM, Choyke PL. Imaging of angiogenesis: from microscope to clinic. Nat Med 2003;9:713–725. 3. Kety SS. The theory and applications of the exchange of inert gas at the lungs and tissues. Pharmacol Rev 1951;3:1– 41. 4. Tofts PS, Kermode AG. Measurement of the blood-brain barrier permeability and leakage space using dynamic MR imaging. 1. Fundamental concepts. Magn Reson Med 1991;17:357–367. 5. Tofts PS. Modeling tracer kinetics in dynamic Gd-DTPA MR imaging. J Magn Reson Imaging 1997;7:91–101. 6. Ludemann L, Wurm R, Zimmer C. Pharmacokinetic modeling of GdDTPA extravasation in brain tumors. Invest Radiol 2002;37:562–570. 7. St Lawrence KS, Lee TY. An adiabatic approximation to the tissue homogeneity model for water exchange in the brain: I. Theoretical derivation. J Cereb Blood Flow Metab 1998;18:1365–1377. 8. St Lawrence KS, Lee TY. An adiabatic approximation to the tissue homogeneity model for water exchange in the brain: II. Experimental validation. J Cereb Blood Flow Metab 1998;18:1378 –1385. 9. Leach MO, Brindle KM, Evelhoch JL, Griffiths JR, Horsman MR, Jackson A, Jayson G, Judson IR, Knopp MV, Maxwell RJ, McIntyre D, Padhani AR, Price P, Rathbone R, Rustin G, Tofts PS, Tozer GM, Vennart W, Waterton JC, Williams SR, Workman P. Assessment of antiangiogenic and antivascular therapeutics using MRI: recommendations for appropriate methodology for clinical trials. Br J Radiol 2003;76(Spec No 1):S87–S91. 10. Evelhoch JL, Gillies RJ, Karczmar GS, Koutcher JA, Maxwell RJ, Nalcioglu O, Raghunand N, Ronen SM, Ross BD, Swartz HM. Applications of magnetic resonance in model systems: cancer therapeutics. Neoplasia 2000;2:152–165. 11. Knopp MV, Giesel FL, Marcos H, von Tengg-Kobligk H, Choyke P. Dynamic contrast-enhanced magnetic resonance imaging in oncology. Top Magn Reson Imaging 2001;12:301–308. 12. O’Connor JP, Jackson A, Parker GJ, Jayson GC. DCE-MRI biomarkers in the clinical evaluation of antiangiogenic and vascular disrupting agents. Br J Cancer 2007;96:189 –195. 13. Padhani AR. Functional MRI for anticancer therapy assessment. Eur J Cancer 2002;38:2116 –2127. 14. Padhani AR. Dynamic contrast-enhanced MRI in clinical oncology: current status and future directions. J Magn Reson Imaging 2002;16: 407– 422. 15. Armitage P, Behrenbruch C, Brady M, Moore N. Extracting and visualizing physiological parameters using dynamic contrast-enhanced magnetic resonance imaging of the breast. Med Image Anal 2005;9:315–329.

1119 16. Buckley DL. Uncertainty in the analysis of tracer kinetics using dynamic contrast-enhanced T1-weighted MRI. Magn Reson Med 2002;47: 601– 606. 17. Tofts PS, Brix G, Buckley DL, Evelhoch JL, Henderson E, Knopp MV, Larsson HB, Lee TY, Mayr NA, Parker GJ, Port RE, Taylor J, Weisskoff RM. Estimating kinetic parameters from dynamic contrast-enhanced T(1)-weighted MRI of a diffusable tracer: standardized quantities and symbols. J Magn Reson Imaging 1999;10:223–232. 18. Henderson E, Rutt BK, Lee TY. Temporal sampling requirements for the tracer kinetics modeling of breast disease. Magn Reson Imaging 1998;16:1057–1073. 19. Li X, Feng D, Wong K. A general algorithm for optimal sampling schedule design in nuclear medicine imaging. Comput Methods Programs Biomed 2001;65:45–59. 20. Horsfield MA, Morgan B. Algorithms for calculation of kinetic parameters from T1-weighted dynamic contrast-enhanced magnetic resonance imaging. J Magn Reson Imaging 2004;20:723–729. 21. Evelhoch JL. Key factors in the acquisition of contrast kinetic data for oncology. J Magn Reson Imaging 1999;10:254 –259. 22. Lopata RG, Backes WH, van den Bosch PP, van Riel NA. On the identifiability of pharmacokinetic parameters in dynamic contrast-enhanced imaging. Magn Reson Med 2007;58:425– 429. 23. Ljung L. System identification: theory for the user. New Jersey: Prentice Hall; 1999. 24. Girod B, Rabenstein R, Stenger A. Signals and systems. New York: John Wiley & Sons; 2001. 592 p. 25. Stollberger R, Aschaur M, Ebner F. Prediction of contrast bolus shape using test-bolus measurements. In: Proceedings of the 13th Annual Meeting of ISMRM, Miami Beach, FL, USA, 2005 (Abstract 377). 26. de Lussanet QG, Backes WH, Griffioen AW, Padhani AR, Baeten CI, van Baardwijk A, Lambin P, Beets GL, van Engelshoven JM, Beets-Tan RG. Dynamic contrast-enhanced magnetic resonance imaging of radiation therapy-induced microcirculation changes in rectal cancer. Int J Radiat Oncol Biol Phys 2005;63:1309 –1315. 27. Haacke EM, Brown RW, Thompson MR, Venkatesan R. Spun density, T1 and T2 quantification methods in MR imaging. Editors: Haacke EM, Brown RW, Thompson MR, VenKatesan R. New York: Wiley; 1999. pp. 637– 667. 28. Roberts C, Buckley DL, Parker GJ. Comparison of errors associated with single- and multi-bolus injection protocols in low-temporal-resolution dynamic contrast-enhanced tracer kinetic analysis. Magn Reson Med 2006;56:611– 619. 29. Port RE, Knopp MV, Brix G. Dynamic contrast-enhanced MRI using Gd-DTPA: interindividual variability of the arterial input function and consequences for the assessment of kinetics in tumors. Magn Reson Med 2001;45:1030 –1038. 30. Kuo PH, Kanal E, Abu-Alfa AK, Cowper SE. Gadolinium-based MR contrast agents and nephrogenic systemic fibrosis. Radiology 2007;242: 647– 649.