A functional form for injected MRI Gd-chelate contrast agent concentration incorporating recirculation, extravasation and excretion

A functional form for the vascular concentration of MRI contrast agent after intravenous bolus injection was developed that can be used to model the concentration at any vascular site at which contrast concentration can be measured. The form is based on previous models of blood circulation, and is consistent with previously measured data at long post-injection times, when the contrast agent is fully and evenly dispersed in the blood. It allows the first-pass and recirculation peaks of contrast agent to be modelled, and measurement of the absolute concentration of contrast agent at a single time point allows the whole time course to be rescaled to give absolute contrast agent concentration values. This measure of absolute concentration could be performed at a long post-injection time using either MRI or blood-sampling methods. In order to provide a model that is consistent with measured data, it was necessary to include both rapid and slow extravasation, together with excretion via the kidneys. The model was tested on T1-weighted data from the descending aorta and hepatic portal vein, and on T*2-weighted data from the cerebral arteries. Fitting of the model was successful for all datasets, but there was a considerable variation in fit parameters between subjects, which suggests that the formation of a meaningful population-averaged vascular concentration function is precluded.


Introduction
Measurement of the time course of vascular concentration of injected contrast agent is important in a number of applications, including dynamic susceptibility contrast (DSC) cerebral perfusion imaging and dynamic contrast-enhanced (DCE) MRI (Calamante et al 2007, O'Connor et al 2007. In the former, a susceptibility-weighted sequence is used, with the increase in R * 2 being proportional to the contrast agent concentration ([Gd]). In the latter, T 1 -weighted sequences are generally used, and the increase in R 1 is also proportional to [Gd]. Accurate assessment of the vascular concentration allows physiological parameters, such as perfusion or vascular permeability, to be evaluated.
However, accurate measurement of the contrast agent concentration is not straightforward. First, after the initial injection, large changes in concentration occur over a period of a few seconds (Larsson et al 1994), so that if the detail of the time course is to be captured, data must be acquired with good time resolution, of the order of one image per second. Secondly, the blood that carries the contrast agent is flowing rapidly, leading to a non-steady state of magnetization and problems in quantification, particularly in T 1 -weighted imaging. Thirdly, the concentration must often be measured in small vessels, where a single pixel's area may be larger than the vessel, leading to partial-volume dilution of the estimated concentration.
Some methods of estimating perfusion parameters rely on estimating the area under the first-pass concentration-time curve, which can be done by fitting a gamma-variate function to the early part of the curve, avoiding the recirculation phase as much as possible. However, this approach is not without problems, since recirculation occurs after a very short period of the order of 10-20 s. A better approach would be to fit a more complete model for the circulating contrast agent, from which the first-pass parameters could then be extracted.
The motivation for our work was threefold. First, we wished to provide a description of contrast agent bolus recirculation, dispersion and redistribution that is based on physiological factors, rather than heuristics. Secondly, we aimed for a model that is compatible with existing simple models. Thirdly, we wanted a model that could either be fitted to individual subjects' data or could potentially be used to describe a data from a population of subjects, as appropriate. The advantage of fitting a functional form is that the fit then can be more accurately interpolated, integrated or physiological factors extracted.
In this paper, we present such a model, based on physiological factors and previously observed parameters. Along with recirculation of the contrast agent, it includes extravasation and excretion via the kidneys that have been experimentally observed previously for Gd-DTPA (Weinmann et al 1984), and in addition it takes into account a postulated rapid extravasation that occurs after injection, predominantly into the visceral organs. Our model is compatible with a previous bi-exponential model of contrast agent concentration that pertains for times longer than a few minutes, when the contrast agent is evenly distributed in the blood (Tofts and Kermode 1991). The model is formulated in a novel way, that enables the recirculation to be decoupled from extravasation and excretion, so that these two aspects of contrast agent redistribution can be examined separately, and then simply combined at the final stage.
In a similar vein, Parker et al (2006) developed an empirical function and fitted it to a dataset which was from the average of 113 separate measurements of arterial concentration data in 23 subjects. The model was chosen without any theoretical considerations of the expected form, and consisted of a mixture of two Gaussians plus an exponential modulated with a sigmoid function, and captured the main features of the first-pass and recirculation peaks, and a slow decay at long times. Their aim was to create a high-temporal-resolution model AIF that could be used when AIF measurements from individual patients were not possible. The function was not fitted to individual subjects' data, whereas one of the aims of the current work was to develop a functional form that has a smaller number of variable fit parameters so that the fitting procedure would be more reliable when used with individual subjects' data.
Previous work has considered compact functional forms, such as those presented by Orton et al (2008), using exponential, gamma-variate and cosine descriptions of the concentration curve. While these forms provide high computational efficiency, we show in the current work that a higher fidelity to the measured AIF is achieved with our model. We also demonstrate that the model can be applied to both the arterial side of the perfused organ, and to the venous drain.

Recirculation model
We first consider a model of a closed circulatory system where the contrast agent is injected intravenously, as is commonly the case. After the injection is complete, if there were no leakage through the endothelial wall or excretion, the total amount of contrast agent within the vasculature would remain constant. The site for intravenous injection is commonly the antecubital vein, and we assume that the contrast agent is injected as a short bolus. From there, the contrast passes through the cardiopulmonary circulation, and is pumped through the arterial system to a point where it may be measured using MRI methods before the blood enters the organ of interest. The venous system draining the organ then feeds back to the heart as part of the circulatory cycle. From the ascending aorta, there are multiple routes that the blood may take on its way to perfuse the different limbs and organs, and of particular note is the branch that includes the portal circulation to the liver.
Each part of this complex route involves dispersion of the bolus of contrast agent as it is diluted in the blood and as it follows the many possible paths around the circulation. Normally, after a period of time of the order of a couple of minutes, the contrast agent is distributed with a fairly uniform concentration throughout the circulation. The duration of the contrast agent bolus injection is usually between 3 and 10 s, which can be considered quite long in comparison to the time scale of dispersion (see later). Therefore, the injection procedure can be considered as another initial mixing process that occurs before the first pass.

Circulation in a closed system
The first part of the model involves the path from the venous injection site to the arterial measurement point via cardiopulmonary circulation. It was shown by Davenport (1983) that a path that involves multiple serial opportunities for mixing and dispersion may be modelled as a series of 'theoretical mixing chambers', resulting in dispersion that is described by a gamma-variate function: where is the gamma function, the extension of the factorial function to non-integer values.
The area under the above curve is k, with k being dependent on the dose of contrast agent given and the volume of plasma accessible to the contrast agent. In the expression above, the coefficient α reflects the number of 'theoretical mixing chambers' and β is the ratio of the volume of a theoretical mixing chamber to the volume flow rate. In the 'first pass' of contrast agent, before any recirculation is considered, a simple gamma-variate function therefore applies.
From the point of measurement, the blood then passes through the organ and venous system before recirculating through the cardiopulmonary system to the measurement point. The recirculation will generate further dispersion of the bolus which can be described by another gamma-variate function. There are, of course, many possible routes that the circulation can take through different organs and limbs. In general, each of these routes will have a different recirculation time and different coefficients α and β in the gamma-variate function for that path. However, for simplicity, we assume here that this complex superposition of recirculation paths can still be described by a single gamma-variate function. Furthermore, in order to reduce the number of model parameters, and to allow the straightforward inclusion of recirculation, we also assume that α and β have the same value as for the initial first-pass phase. Since β depends on the time required for a theoretical mixing chamber to drain, this may not be an unrealistic assumption, since the 'mixing chamber' is an abstract concept, and β is the ratio of the volume of the chamber to the rate of flow through it. If the bulk of the mixing takes place in a capillary bed, and capillaries are of a fairly uniform length, then β is simply inversely proportional to the flow velocity. For capillary beds, we assume that the flow velocity through an individual capillary is relatively independent of the tissue type, although we recognize that it may be highly variable and dependent on many physiological factors. However, it is not practical to incorporate this level of detail into the model. Also, if the bulk of dispersion occurs in the path from the injection site to the point of measurement, it is reasonable to assume the same values for α and β during the recirculation as for the first pass. Davenport (1983) showed how recirculation in a closed loop can be described by setting the input to the loop in the first recirculation equal to the output of the loop at the end of the first passage around the loop. This can be formulated mathematically by convolving the gamma-variate function by a time-shifted version of the same function, with the time shift being equal to the minimum recirculation time. The result of this operation is added to the original function. Each subsequent recirculation is described by a further convolution with the time-shifted gamma variate, with the time shift being equal to the number of recirculations multiplied by the minimum recirculation time. The sum of all such convolved functions means that a uniform concentration in the blood is eventually reached.
If the minimum time required for blood to recirculate is τ , then we have For the first recirculation, the contrast agent present is that which remains from the first pass, plus that which has recirculated: Here, * is the convolution operator. The result of the convolution of two gamma-variate functions with the same β parameter is another gamma-variate function (Davenport 1983) such that equation (3) can be rewritten as In general, for N recirculations (N 0) After several circulations, the contrast agent is evenly distributed throughout the vasculature. We can use the final value theorem of the Laplace transform to obtain the steady-state concentration: where F(s) is the Laplace transform of C p (t). Taking the Laplace transform of equation (5) gives As N → ∞, Thus, the steady-state concentration is

Extravasation and excretion
Expression (5) does not fit the observed data for MRI contrast agent injected as a bolus in the antecubital vein. In particular, the area under the first-pass gamma-variate function is considerably greater than that under the first and subsequent recirculations. In a closed system, the area under each of the gamma-variate functions in equation (5) is constant, and any variation in the area must be due to loss of contrast from the vascular tree due to extravasation and excretion. For times greater than that for a uniform distribution of contrast agent to pertain, Weinmann et al (1984) used blood sampling and measured a bi-exponential decay of the plasma concentration of Gd-DTPA in normal healthy volunteers, which was attributed to extravasation and excretion. Tofts and Kermode (1991) fitted a bi-exponential curve to Weinmann's data and found that the plasma concentration (in mM) could be expressed as where D is the dose of contrast agent in mmol per kg of body weight and t is the time (in seconds) from contrast injection. If the contrast agent were uniformly distributed straight after contrast injection, this would give an initial plasma concentration of 8.77D mM. With a standard dose of contrast agent (0.1 mmol kg −1 ), we would therefore expect an initial plasma concentration of approximately 0.877 mM. However, taking the example of a healthy adult male with a weight of 75 kg, a whole blood volume of 5 L and a haematocrit value of 0.45, the dose given would be 7.5 mmoles, and if all contrast agents remain in the blood the initial plasma average concentration (at t = 0) is expected to be 2.72 mM. This discrepancy has also been noted in animal studies using the radiopharmaceutical 99m Tc-DTPA, which has a pharmacokinetic profile that is nearly identical to that of Gd-DTPA (Wedeking et al 1990). In a multi-compartment model in the rat, the compartment assigned to plasma was approximately three times the true plasma volume. However, since blood sampling started after 1 min, any very rapid extravasation could not have been captured. Clearly, for these studies, there must have been considerable initial extravasation of contrast agent over the short time between injection and the start of the measurement of arterial concentration. This is seen particularly in abdominal contrast-enhanced MR angiography, where there is loss of image contrast after the first pass, not only due to a reduction in vessel concentration but also because of obscuration of the vessels as the contrast agent rapidly extravasates into the surrounding tissue. We therefore include an extra (fast) compartment into which the contrast agent extravasates, which we expect to result in a tri-exponential decay of the average plasma concentration (Upton 2004). We expect that the amplitude of this extra component will be considerably greater than the amplitudes of the other two seen by Weinmann, and that the rate constant will be much shorter than the other two.
Although the amount of extravasation will vary according to the location in the vasculature, we do not include that level of detail in the model. We simply assume that any loss of contrast agent is uniformly distributed around the system, thereby allowing the rate of loss to be dependent only on the mean concentration of contrast agent in the vasculature. Clearly, contrast agent will extravasate into multiple extravascular compartments of different sizes and with different rates, but we do not have sufficient quality of data to include more than one extra exponential extravasation term, and therefore we model all processes for contrast leakage from the vasculature with two rate constants. Also, different parts of the bolus will pass through different organs at different times, so that extravasation will occur at different rates for different parts of the bolus.
With the assumption of uniformly distributed contrast agent loss from the vasculature, any part of the bolus that passes the measurement point will therefore have been subject to a dilution that depends only on the time since injection which, for a tight bolus injection, is the same for all parts of the contrast agent bolus. The recirculation is thereby decoupled from the extravasation and excretion, and we can simply multiply the solution for the closed system (equation (5)) by a tri-exponential decay to obtain Nτ <t <(N+ 1)τ. (11) In our model, the amplitude and time constant for the first (rapidly decaying) component (A 1 and m 1 ) are determined from the data, whilst the other two are those given in equation (10), obtained from Tofts and Kermode (1991).

Accounting for systematic errors in the relaxation rate
In order to evaluate organ perfusion or vascular permeability, the vascular concentration of contrast agent in the vessel feeding the organ of interest must be evaluated. One of the problems of measuring this arterial input function (AIF) from MRI data is that the feeding vessels are often of a size comparable to the image pixel size. Therefore, it may be necessary to measure the AIF from single pixels, or from small regions that contain significant amounts of signal from neighbouring tissues. This leads to a 'dilution' of the estimated concentration of contrast agent, and it is difficult to obtain quantitative concentration values. Further pulse sequenceand flow-related systematic errors in concentration estimation are also likely. All these types of error can result in an estimated contrast agent concentration that is not proportional to the true concentration. For example, the partial volume dilution would result in a greater percentage underestimate of concentration at high concentrations.
However, if we assume, to a first approximation, that the proportionality between estimated concentration and true concentration is maintained, this issue can be addressed by scaling the AIF by a normalization factor, k N , such that for times longer than that required for the contrast agent to become evenly distributed in the blood, and the rapid extravasation phase is over, the concentration matches that found by Weinmann et al (1984). The uniform plasma concentration without extravasation or excretion is given by equation (9), and so it is necessary that or The final expression for the plasma concentration is then with t 0 corresponding to the time at which contrast agent is first seen in the vessel, and C p (t) = 0 for t < t 0 . In practice, equation (11) is fitted to the concentration-time curve, and then the concentration values are rescaled by multiplying by the normalization factor k N calculated according to equation (13). This rescaling ensures that, after the initial rapid extravasation stage and when the contrast agent is uniformly distributed, the concentration matches that observed by Weinmann. Equation (14) can then be used to describe the concentration-time curve.

Materials and methods
All data were acquired at 1.5 Tesla using two different MRI scanners, both General Electric Signa machines (GE Medical Systems, Waukesha, WI, USA). The model was tested on data acquired in the abdomen using a single-slice T 1 -weighted sequence, and in the head using a multi-slice T * 2 -weighted EPI sequence. The contrast agent was injected at a dose 0.1 mmol per kg of body weight via the antecubital vein in both cases using a power injector at a rate of either 3 mL s −1 (for the abdominal study) or 5 mL s −1 (for the head), with a saline flush, corresponding to bolus injection durations of around 5 or 3 s, respectively. The contrast agent used for the abdominal study was gadolinium diethylenetriamine penta-acetic acid (Magnevist, Bayer Schering Pharma AG, Germany) and for the brain study was gadoterate meglumine (Dotarem, Guerbet, Roissy CDG, France). All subjects gave informed consent for the MRI procedure which was approved by the respective local research ethics committees.

T 1 -weighted imaging of the liver
Our first test of the applicability of the model was in the liver, where a modified T 1 -weighted saturation-recovery snapshot-FLASH sequence was used to quantify the relaxation rate (Black et al 2007). The sequence acquires separate images at two recovery times, suitable for measuring contrast concentration in the artery and in the tissue, so the dynamics of contrast agent uptake can be measured as part of an assessment of organ vascular permeability (Tofts et al 1999). Contrast agent was injected approximately 20 s after the start of the scan. The sequence parameters were as follows: TR = 4.8 ms; TE = 1.3 ms; flip angle = 12 • ; 256 × 128 matrix, slice thickness = 10 mm; parallel imaging factor = 2; centric phase-encode ordering. Two independently oriented images per heartbeat were acquired for 200 heartbeats. After a 90 • non-selective saturation pulse, the first of each pair of images was acquired after a recovery time (T rec ) of 50 ms, and was oriented to cut through the aorta and portal vein. The second was acquired after a recovery time of 150 ms, and was oriented axially to cut through the liver. Data from four healthy volunteers were acquired; see Upton (2004) for further details of the acquisition protocol and table 1 for demographic details.
The contrast agent concentration was measured by placing a region-of-interest (ROI) slightly within the walls of the aorta (area approximately 335 mm 2 ), and another within the . For the period after contrast agent injection, signal intensities were converted to changes in the relaxation rate using the expression where S(t) is the signal intensity measured during the dynamic run, S 0 is the signal intensity measured without the saturation pulse and with full relaxation and R 10 is the relaxation rate measured during the period before contrast agent injection. R 10 was estimated in the same fashion using S 0 and the data acquired as part of the dynamic run, but before contrast agent had been injected.
T * 2 -weighted imaging of the brain The second test of the applicability was to brain perfusion data, where a T * 2 -weighted EPI sequence was used to quantify the relaxation rate. Thirteen patients with low-grade glioma (WHO grade II) underwent perfusion MRI as part of a longitudinal multimodal MRI research study (Danchaivijitr et al 2008). Twelve subjects were selected for analysis based on there being one or more scan planes covering vessels from which the AIF could be measured. A gradient-echo EPI sequence was used (TR = 1200 ms; TE = 40 ms; flip angle = 20 • ; FOV = 26 cm; matrix = 96 × 128; slice thickness = 5 mm; six slices). A low flip angle was used to reduce the degree of signal enhancement due to T 1 -weighting. Contrast agent was injected after approximately two scans had been acquired.
The contrast agent concentration was measured using small ROIs manually placed in the circle of Willis or the middle cerebral arteries, being careful to avoid areas of tumour. Signal intensities were converted to changes in relaxation rate values using the expression where S(t) is the signal intensity measured during the dynamic run and S 0 is the average signal intensity measured for the period of the scan before contrast agent injection.

Contrast agent concentration
For both T 1 -weighted and T * 2 -weighted imaging, the contrast agent concentration was estimated from the expression where R is the change in relaxation rate (either R 1 or R * 2 ) and ρ is the appropriate relaxivity. For the T 1 -weighted data, ρ 1 was assumed to be 4.33 s −1 mM −1 (Fritz-Hansen et al 1996). For the T * 2 -weighted data, the relaxivity ρ 2 * is unknown, and was assigned an arbitrary value of 1.0 s −1 mM −1 such that absolute quantification of plasma concentration was not attempted. Furthermore, the contrast agent used in the brain study was Gd-DOTA, for which the relative amplitudes and time constants in equation (10) are not known; however, we assumed them to be the same as for Gd-DTPA (Weinmann et al 1984), since the two compounds have a similar molecular weight and structure. For the liver data, the blood haematocrit was assumed to be 0.5 for male subjects and 0.45 for females (Purves et al 2004).
The model equation (11) was then fitted to the data using the Levenberg-Marquardt iterative least-squares minimization algorithm (Press et al 1992). This algorithm requires initial guesses to be set for each of the fitted variables, and in order to permit successful fitting, it was necessary to first estimate t 0 by finding the first time point at which 20% of the maximum change in the contrast agent concentration from the initial (pre-injection) value occurred. The initial guess for t 0 was then taken as the previous time point. All other fit variables (α, β, m 1 , A 1 ) were initialized to physically plausible values, with the same initial value for all subjects. Based on the observation of the recirculation and extravasation behaviour, the initial values set were α = 4, β = 10 s, m 1 = 15 s −1 and A 1 = 16 kg L −1 . k N was estimated after fitting using equation (13) and used to rescale the concentration values.
In order to compare the applicability of the new model to current alternatives, we fitted two existing functions to the aorta data. We did not fit these to the portal vein data, since the two models tested were only intended to be applied to arterial data. The first was 'Model 3' used in Orton's paper, which is defined piece-wise and basically models the initial first-pass peak using a raised cosine function, and with an exponential decay at long times (Orton et al 2008). In addition to the time of first contrast appearance, t 0 , the model has four free fit parameters. The second was the one used by Parker to create the high temporal resolution population-averaged AIF (Parker et al 2006). This model has ten free fit parameters. Note, however, that in Parker's paper, the model was fitted to a single dataset created by averaging data from many subjects rather than fitting to individual subjects' data. To test whether the new model provided a better description of the data compared to these existing model, the F-statistic was used, as presented by Boxenbaun et al (1974). F-statistic values were compared to F-table values to determine whether a reduction in the summed-squared error between the fitted model and the data was statistically significant between models.

Results
The quality of data was found to be insufficient to provide a stable estimate of τ , the minimum blood recirculation time, for the data from the portal vein, due to the slower passage and large overlap of the recirculation peaks. τ has been estimated previously to be approximately 5 s in rats and 7.4 s in dogs using radioactive tracer methods (Tothill andMacPherson 1980, Henze et al 1982). For the four datasets from the aorta, the mean fitted value of τ was 10.6 s; τ was therefore fixed at this value and resulted in visually acceptable fits in all cases. With good initial guesses for the remaining fit parameters, the fitting algorithm successfully converged on an optimum fit for all datasets on which it was tested.
We first illustrate the need to include another compartment (a third exponential decay term) as described in equation (11). Figure 1 shows a 'reduced' model fitted to the data from the abdomen of the first subject (aorta and portal vein) study. This model incorporates only the two exponential decay terms of the model by Tofts (equation (10)). In order to better fit to the long-time concentration data, the initial first-pass peak has been considerably underestimated, Figure 1. Illustration of the need to include a third exponential decay term to describe the extravasation of contrast agent. A recirculation model with the bi-exponential extravasation and excretion amplitudes and time constants estimated by Tofts and Kermode (1991) from Weinmann's data (Weinmann et al 1984) has been fitted to the plasma concentration values. First-pass amplitude is underestimated, long-time concentration overestimated and the multiple recirculation peaks have much larger amplitude and persist for a much longer time than is indicated in the data.
suggesting that there is a substantial loss of contrast agent from the plasma over the time-scale of a few recirculations. Imposing this constraint has also forced a shift of the first recirculation peak to the right. Figure 2 shows estimated [Gd] data from the aorta and portal vein of four patients, together with the tri-exponential recirculation model fit to the data (equation (11)), but without rescaling the concentration to match the long-time behaviour as detailed in equation (13). In figure 3, the fitted curves have been rescaled by multiplying by the normalization factor calculated using equation (13) so that their long-time behaviour matches the bi-exponential form of Tofts and Kermode (1991); all aorta measurements are superimposed on one plot, and all portal vein data on another. It would appear by comparing figures 2 and 3 that the plasma concentration is being overestimated by a factor of approximately 2. Also, in figure 3, the time origin has been shifted so that t 0 coincides with t = 0. The fitted model parameters are shown for the four patients in table 1.
Note the more consistent peak amplitude for the aorta during the first-pass of contrast agent when data are rescaled according to the expected concentration at long times. The model can accommodate data with the very different forms that can be seen when comparing aortic data with that from the portal vein. There is also considerable patient-to-patient variation in the original data, particularly noticeable for the portal vein, which is also reflected in the fitted forms and fit parameters shown in table 1.
For the images acquired for the brain perfusion study, because of the short TR of the EPI sequence, we used a low flip angle excitation pulse to reduce T 1 -weighting. Nevertheless, with contrast agent present, the signal intensities were still higher than would have been the case for a pure T * 2 -weighted sequence, and the contrast agent concentration was underestimated, which was particularly apparent during the late phase of contrast passage for some patients (Calamante et al 2007). Thus, quantification of the contrast agent concentration and rescaling to true [Gd] were not possible, and analysis was restricted to fitting the unscaled model, equivalent to equation (11), which is shown for all 12 subjects in figure 4. The means (SD) of the fitted model parameters were α = 3.98 (1.42), β = 2.64 (0.66) s, m 1 = 0.075 (0.030) s −1 , and A 1 = 69.9 (31.2) kg L −1 .
For each of the four datasets from the aorta, the new model better described the simple functional form of Orton's 'Model 3' (p < 0.0001). This was also the case for Parker's model (p < 0.0001). However, in all cases, Parker's model failed to converge to a satisfactory fit, assessed by visual inspection, despite the fit parameters being initialized to the values of the population-averaged model.

Discussion
The proposed model for the concentration of circulating MRI contrast agent is designed to take account of bolus dispersion, recirculation, extravasation and excretion whilst maintaining compatibility with previous work that only modelled the longer time behaviour (greater than 1 min) (Weinmann et al 1984, Tofts andKermode 1991). In cases where an accurate relaxivity of the contrast agent is not known, or where the relaxation rate cannot be accurately measured, the model should allow an accurate determination of absolute concentration during the firstpass and subsequent recirculations, as long as the long-time concentration is known. This assumes that the estimated relaxation rate is proportional to the true relaxation rate in the blood, and in turn that the true relaxation rate is proportional to the blood concentration of contrast agent. The former may not be a realistic assumption when there are inflow effects, partial volume dilution or when there is poor calibration of RF pulse flip angles for MRI data acquisition. There is still potential for error in estimating the high concentration seen during the first pass because of the nonlinear relationship between relaxation rate and signal intensity seen for MRI pulse sequences when the relaxation rate is high. The errors in estimates of R 1 generally increase at high concentrations, although the pulse sequence we used for the abdominal study was particularly robust in this respect (Black et al 2007).
In our data from the aorta and portal vein, we have used the model to establish a connection between arterial and venous estimates, by ensuring consistency with a previous bi-exponential model under the assumption that at long times the contrast agent is evenly distributed throughout the blood. Any differences between the estimates of the contrast agent concentration in the aorta and portal vein at long times must be artefactual and possibly related to the different flow velocities in these vessels or vessel sizes (leading to different partial volume dilution). At long post-injection times, the model is constrained to have known concentration as determined by Weinmann et al (1984), and thus is capable of correcting for these types of systematic error, and providing absolute estimates of concentration. Although not readily apparent from figure 3, the rescaling brings the arterial and venous concentrations in line with each other at long post-injection times; the very slow passage through the gastric circulation means that this is not reached in the portal vein until several minutes have elapsed. Patientto-patient variation from the bi-exponential model could partly be accounted for making an accurate measurement of [Gd] at a single time point. This could correct for differences in the injected dose and overall volume accessible to the contrast agent, but could not take account of any differences in the slow extravasation and excretion time constants and relative amplitudes. This accurate measurement could either be by blood sampling, or by using a lengthier but more accurate relaxation time measurement procedure after the dynamic study had been performed.
Further examination of figure 3 shows corrected peak plasma concentrations of around 5-6 mM in the aorta, which are in line with those measured by Fritz-Hansen using direct arterial sampling (Fritz-Hansen et al 1996). However, they are considerably lower than the population-averaged model constructed by Parker and colleagues (Parker et al 2006) (figure 3) by MRI methods, where peak concentrations are similar to the uncorrected values as shown in figure 2. This difference might be attributable to the type of pulse sequence that was used both in our study and in Fritz-Hansen's, a 2D snapshot FLASH acquisition that might be more susceptible to inflow effects than the 3D acquisition used by Parker. However, there was a considerable variation in Parker's individual subjects' measurements, possibly due to the lower temporal resolution of 4.97 s (compared to ∼1.0 s for our study and ∼3.0 s for Fritz-Hansen's) being unable to capture the peak concentration for all subjects. Note that for the short post-injection time shown in figure 3, Parker's population-average model has concentrations that are higher than those measured by Weinmann, but that the values converge after approximately 6 min.
The fitting procedure was not able to estimate τ , the minimum recirculation time, reliably. This is most likely because the start of the first recirculation occurs right under the peak of the first pass, and the time from the start of the first recirculation to the first recirculation peak is strongly influenced by the fit parameters α and β, making the estimate of τ sensitive to them. Despite careful initialization of the fit parameters, Parker's heuristic model failed to fit satisfactorily to the aorta data, the fitting procedure appearing to converge on a false (local) minimum. However, the model was not intended to be fitted to individual subject's data, but rather to a dataset derived by averaging data from many subjects. Thus, the fact that the fit residuals were smaller for our model than for Parker's should not be taken as an indication of the superiority of our model. However, it does indicate that fitting of a model with fewer fit parameters, which still adequately describes the data, is more reliable. The new model provided a significantly better description of the concentration-time curve than did Orton's four-parameter model, with only one extra fit parameter.
Accurate and realistic modelling the vascular concentration function is advantageous because it allows accurate interpolation of the concentration-time curve, extrapolation beyond the end of the measured data and, if desired, analytical integration of the area under the curve. Furthermore, the first pass of concentration can be extracted from the model parameters. Thus, since some analysis procedures require the first pass to be isolated from recirculated contrast agent, this can be done accurately by fitting to the whole dataset, rather than, for example, truncating the data after a short time, and fitting a simple gamma-variate to this short time segment of data. To illustrate this, we fitted a gamma variate function to the brain AIF data after truncating it at the time of the trough in concentration between the first-pass peak and first-recirculation peak. The mean fitted α and β values were 2.67 and 3.17 respectively for the simple gamma-variate fit, compared to 4.03 and 2.65 respectively for the full model. Since the plasma contains recirculating contrast agent after the minimum recirculation time, which we estimated to be approximately 10 s, there is a 'contamination' of the first-pass peak well before the trough in concentration, which typically occurs 20-25 s after the contrast agent is first seen in the vessel. This accounts for the higher α and lower β values, as the 'tail' of the gamma variate function is extended to accommodate this higher concentration.
In order to provide a model that faithfully accounts for the shape of the vascular concentration, it was necessary to include a third exponential decay term to accommodate the rapid extravasation of contrast agent that occurs in the early phase of contrast passage. Table 1 and m 1 values presented for the cerebral data show that this component has a decay time constant of around 15-20 s, which accounts for it not being observed by Weinmann, whose arterial sampling started only at 1 min post-injection, after this component has largely decayed. In a study in rats using 99m Tc-DTPA, a radiopharmaceutical that has a pharmacokinetic profile virtually identical to that of Gd-DPTA, Wedeking et al (1990) showed that a three-compartment model was necessary to adequately describe the measured plasma concentration using blood sampling that started 1 min post-injection. However, the extra compartment had a time constant of the order of 1 min, and the compartment they assigned to plasma had an estimated distribution volume approximately three times that of the true plasma volume, as was also shown by Weinmann. This large plasma distribution volume remained unexplained by both Weinmann and Wedeking; we postulate here that the contrast agent extravasates into a further compartment with a very short time constant, so that it is unobservable by blood sampling methods. The decay due to extravasation in the Tofts and Kermode model is relatively slow, and was assumed to correspond to extravasation to large compartments of the body such as muscle. Our additional decay is much more rapid, is probably related to extravasation in the abdominal viscera, including the kidneys, and is an average of several decay rates relating to all organs, including the gastro-intestinal tract. The amplitude and rate of enhancement in these organs is considerably higher both observationally and experimentally than in muscle (Medved et al 2004).
Weinmann shows normal variation in the plasma levels of injected gadolinium chelates after 1 min of ±15%, assuming fixed dose per unit of body weight (Weinmann et al 1984). However, we found considerable variability for the first minute of enhancement in our normal volunteers, particularly notable in the portal vein fit parameters. Variation in the cardiac output as a fraction of the total blood volume, can account for some of the variability of the first-pass and recirculation peaks in the arterial data. However, the wide variation seen in the portal vein may be because, although the renal and slow decay rates are relatively constant, the rapid decay rate is related to the splanchnic circulation and abdominal viscera. It is well known that eating or drinking caffeinated liquids produces considerable changes in the splanchnic circulation (Thomsen et al 1993), and this may explain the variability in our normal volunteers compared to clinical scanning protocols, where patients are often restricted to water alone for up to 4 h prior to scanning.
Some of the fitted curves in figure 2 have first recirculation peaks that are less pronounced than those in the data (particularly the aorta data from subjects 2 and 3). The most likely cause is that the fit is being strongly influenced by the first-pass peak, which has perhaps been poorly captured in these subjects. In MRI, the error in the estimate of the contrast agent concentration is considerably greater at high concentration.
At the start of this work, it was envisaged that it would be possible to build a typical functional form for the arterial input function from a population average of the fitted model parameters, in the same way as was done by Parker and colleagues with an empirical model (Parker et al 2006). However, the wide variation in fitted model parameters that we have observed suggests that this may not be feasible, and that direct measurement of the AIF is to be preferred. The importance of assessing the AIF on a subject-by-subject basis was underlined by a recent study, which showed a considerable reduction in the coefficient of variation in derived parameters for DCE-MRI for data-derived AIFs compared to a single common model for the AIF (Ashton et al 2008). Nevertheless, in cases where measurement is not possible, a last resort would be to use such a synthetic population-averaged AIF, although it would require a larger sample size and a broader range of subjects than we have used in order to build the model. It is also likely that the model would vary according to the clinical condition, age and perhaps sex of the subjects.
In summary, we have presented a functional form for the vascular concentration of contrast agent after intravenous bolus injection that can be used to model the concentration at any vascular site at which contrast concentration can be measured. It is based on previous models of blood circulation, but includes an extra, rapid extravasation constant, assumed to be related to the visceral organs. The model is consistent with previously measured data at long post-injection times, when the contrast agent is fully and evenly dispersed in the blood.