Correcting saturation effects of the arterial input function in dynamic susceptibility contrast-enhanced MRI — a Monte Carlo simulation

Available online at www.sciencedirect.com

Magnetic Resonance Imaging 25 (2007) 1300 – 1311

Correcting saturation effects of the arterial input function in dynamic susceptibility contrast-enhanced MRI — a Monte Carlo simulationB Peter Bruneckera,4, Arno Villringera, Jfrg Schultzea,b, Christian H. Noltea Gerhard Jan Jungehqlsinga, Matthias Endresa,b, Jens Steinbrinka a

Berlin NeuroImaging Center (BNIC) and Department of Neurology, Charite´-University Medicine, D-10117 Berlin, Germany Klinik und Poliklinik fu¨r Neurologie, Charite´-Universita¨tsmedizin Berlin, Campus Mitte, Charite´platz 1, D-10117 Berlin, Germany Received 5 October 2006; revised 26 February 2007; accepted 1 March 2007

b

Abstract To prevent systematic errors in quantitative brain perfusion studies using dynamic susceptibility contrast-enhanced magnetic resonance imaging (DSC-MRI), a reliable determination of the arterial input function (AIF) is essential. We propose a novel algorithm for correcting distortions of the AIF caused by saturation of the peak amplitude and discuss its relevance for longitudinal studies. The algorithm is based on the assumption that the AIF can be separated into a reliable part at low contrast agent concentrations and an unreliable part at high concentrations. This unreliable part is reconstructed, applying a theoretical framework based on a transport-diffusion theory and using the bolus-shape in the tissue. A validation of the correction scheme is tested by a Monte Carlo simulation. The input of the simulation was a wide range of perfusion, and the main aim was to compare this input to the determined perfusion parameters. Another input of the simulation was an AIF template derived from in vivo measurements. The distortions of this template was modeled via a Rician distribution for image intensities. As for a real DSC-MRI experiment, the simulation returned the AIF and the tracer concentration-dependent signal in the tissue. The novel correction scheme was tested by deriving perfusion parameters from the simulated data for the corrected and the uncorrected case. For this analysis, a common truncated singular value decomposition approach was applied. We find that the saturation effect caused by Rician-distributed noise leads to an overestimation of regional cerebral blood flow and regional cerebral blood volume, as compared to the input parameter. The aberration can be amplified by a decreasing signal-to-noise ratio (SNR) or an increasing tracer concentration. We also find that the overestimation can be successfully eliminated by the proposed saturation– correction scheme. In summary, the correction scheme will allow DSC-MRI to be expanded towards higher tracer concentrations and lower SNR and will help to increase the measurement to measurement reproducibility for longitudinal studies. D 2007 Elsevier Inc. All rights reserved. Keywords: Perfusion imaging; MRI; Contrast media; Brain; Arterial input function

1. Introduction Dynamic susceptibility contrast-enhanced magnetic resonance imaging (DSC-MRI) is a well-established tool to determine brain perfusion [1]. Especially for the characterization of the acute stage of an ischemic insult, the superior sensitivity for low perfusion states currently hinders a replacement by arterial spin labeling techniques [2,3]. B

This work was supported by the bBundesministerium fqr Bildung und ForschungQ (Berlin NeuroImaging Center, CompetenceNet Stroke). ME was supported by the bVolkswagen StiftungQ (Lichtenberg program). 4 Corresponding author. Berlin NeuroImaging Center (BNIC), Klinik und Poliklinik fqr Neurologie, Charite´ — Campus Mitte, 10117 Berlin, Germany. Tel.: +49 30 450 560 015; fax: +49 30 450 560 952. E-mail address: [email protected] (P. Brunecker). 0730-725X/$ – see front matter D 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.mri.2007.03.011

For maximal signal-to-noise ratio (SNR) and large scanning volumes, DSC-MRI utilizes T 2*-weighted repetitive imaging with a temporal resolution in order of a second [4]. The paramagnetic tracer leads to a concentrationdependent signal decrease described by   k S ðt Þ cð t Þ ¼  log ; ð1Þ TE S ð 0Þ where S(t) is the MR signal, T E the echo time and k a factor depending on the scanning sequence and the MR scanner [4]. Since k is usually unknown without calibration, it was set to T E, and c(t) can consequently only be given in arbitrary units. From this time course, the regional cerebral blood volume (rCBV), the regional cerebral blood flow (rCBF) and mean transit time (MTT) can be calculated [5].

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311

However, these relevant quantities can only be determined if the dissipation of the bolus during its passage from the venous injection to the brain arteries is taken into account, by determining the arterial input function (AIF). Using the concept of the residue function R(t) to describe tracer washout in tissue in combination with a suitable deconvolution approach rCBF can be calculated by solving the following equation [4]: cðt Þ ¼ rCBFd AIFðt Þ  Rðt Þ:

ð2Þ

Various deconvolution algorithms can be routinely used; however, their results — and thus the quality of the calculated parameters — depend on a reliable AIF [6–8]. Other methods used for rCBF quantification that do not need an arterial input function as, for instance, the reference tissue method [9,10], cannot be applied since the prerequisite of a free diffusible tracer is not given for T 2*-weighted DSC-MRI of the brain. To define the AIF, a region of interest (ROI) of arterial voxels has to be selected. For the choice of this arterial ROI, there is a trade-off between using voxels with a medium or a large arterial volume fraction. The former leads to an increased effect by nonarterial vessels surrounding the arteries [11], while the latter can lead to a massive signal drop manifesting by a saturation effect which distorts the AIF [12]. Both mechanisms contribute to the uncertainty of the AIF. However, there are additional effects which can distort the AIF, such as different T 1 and T 2* values in arterial blood as well as in surrounding tissue or different evolution of signal and phase inside and in the vicinity of arteries [13–15]. In summary, all of the effects above leads to a poor accuracy of the AIF and, thus, to an increase of the variance of DSC-MRI. This might be a limitation especially for longitudinal studies. This issue has been recognized, and solutions have been proposed. SNR for arteries and the tissue can be optimized by applying different echo times in distinct slices [16]. In addition, determination of phase images can overcome limitations of magnitude images, especially after applying a correction scheme for partial volume effects [14]. Finally, defining appropriate procedures for automated classification of arterial voxels can reduce variability in absolute quantification [17]. In this work, we propose to treat the damping or saturation of the AIF by pure postprocessing. The main aim of the article is to define and test the reliability of a novel correction algorithm on bolus tracking data. It is based on a parametric description of blood transportation utilizing a transport–diffusion approach. This algorithm separates the AIF into a breliableQ and an bunreliableQ phase, the former at low tracer concentrations, where linearity can be assumed, and the latter at high tracer concentrations at which measurements are potentially distorted due to effects described above. The unreliable phase will be subsequently characterized by a parametric description.

1301

To demonstrate both potential and limitations of the correction procedure, we use a Monte Carlo simulation of the bolus saturation in the tissue. It applies a dedicated template for the AIF function, which allows the second pass and the steady state of the bolus to be considered. A relevant contribution to the simulation is to model the measurement noise via a Rician distribution [18]. We will show that this nonsymmetric distribution is a potential mechanism to describe the observable saturation phenomenon. Further, the tissue response function used in the simulation is an exponential function, spanning a wide range of physiologically meaningful rCBF and rCBV values. By combining the three building blocks (AIF template, noise model, model for tissue response function), a simulation for a DSC-MRI, which closely matches a realistic measurement, is provided. Finally, the simulated data is analyzed using the corrected and the noncorrected AIF. Deriving perfusion parameter and comparing both allows quantification of the improvement. 2. Materials and methods 2.1. The AIF from in vivo MRI experiments For the simulation of a measured tissue response function, an AIF is needed, which closely resembles AIFs determined in vivo. Therefore, six healthy volunteers (mean age, 61.5 years; range, 43–73 years, five females) were examined in an MRI study (1.5T, Vision, Siemens, Germany). Eight transverse slices in anterior–posterior commissure orientation were obtained by a T 2*-weighted echoplanar imaging gradient-echo sequence using the following parameters: field of view =256 mm, T E = 54 ms, T R = 1.3 s, total acquisition time =60 s, acquisition matrix

Fig. 1. The six different ROIs used to define M3 segments of the MCA shown for Subject 2.

1302

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311

power injector (Spectris, Medrad, Indianola, PA, USA). Data acquisition was started at the beginning of the contrast agent injection. All examinations were approved by the ethics committee of the Charite´-University Medicine Berlin. In a further step of analysis, several ROIs were marked manually by an anatomically skilled observer. Due to the large number of bifurcations and, hence, branches, the middle cerebral artery (MCA) was chosen for this task, especially its horizontal (M1) and cortical (M3) segment. For each subject, six different ROIs at the M3 level as well as a single ROI at the M1 level of bifurcation were selected (Fig. 1). The number of six M3 segments was chosen as an upper limit for the number of anatomically comparable regions, which can be found in all individuals. Additionally, this number of M3 ROIs was large enough to realize a robust standard deviation. Fig. 2. (A) Signal time course S(t) during bolus passage determined in different branches of the MCA at a comparable bifurcation level (M3) in one subject (#2). (B) Same as (A), but converted into relative concentrations [c(t)] using Eq. (1) and normalized over the second pass and steady state.

128128 and slice thickness 6 mm. To safeguard for Rician-distributed images, the number of averages was set to one. For DSC-MRI, 20 ml Magnevist (0.5 M GdDTPA, Schering, Germany) followed by 20 ml saline were consecutively injected at a rate of 4 ml/s using a

2.2. Identification of the saturated part of AIF The identification of an unreliable (saturated) part of the AIF is the basis of the correction algorithm, which is described in the following section. We propose to separate the AIF in a reliable phase for low concentration and an unreliable phase. The latter will be corrected by a parametric description. The output of this correction algorithm is an improved AIF, which can be used for further analysis toward the determination of perfusion parameters.

Fig. 3. Standard deviations r of the individual AIFs from 6 subjects after normalization over the second pass and the steady state (bold lines). The average equilibrium level r eq is shown in green. The red box marks the unreliable part of the individual AIF defined via the threshold r trust (red line).

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311

1303

When using the DSC-MRI approach, a single arterial region is commonly chosen to determine the AIF, which is subsequently taken as representative for the whole brain. In contrast, we selected six different arterial regions to allow for a determination of the unreliable phase. The regions were manually selected according the following two criteria: 1. 2.

The time course has an early and steep rising edge. Each region should include several (3–8) adjacent voxels.

As shown in Fig. 2A, for voxels including different branches of the MCA, there are remarkable intrasubject variations over the entire bolus passage. After normalization over the second pass and the steady state — meaning, that the average concentration change there is made identical — the concentration time curves become alike, at the expense of higher variations for the first pass (Fig. 2B). These variations were quantified in Fig. 3 for six subjects by showing the standard deviation r of the AIFs after normalization. Clearly, r is much higher during the first pass than the other parts of the AIF. The standard deviation of the 6 AIFs per subject allows to define the bunreliableQ part of the bolus passage. For this purpose, the average standard deviation for the second pass and the steady state is defined as r eq (green line). If the standard deviation exceeds jeq by a factor of 4 or more (r trust, red line) the unreliable phase is identified (red box), and its boundaries in time are defined as t A and t B. Note that the factor of 4 is a first educated guess, which allows including all diverging AIFs. 2.3. Theory of correcting the AIF Once having identified the unreliable phase, a next step is the replacement by two third-order polynomials as depicted in Fig. 4: pA ð t Þ ¼ a3 t 3 þ a2 t 2 þ a1 t þ a0 pB ð t Þ ¼ b3 t 3 þ b2 t 2 þ b1 t þ b0

ð3Þ

As shown in this figure, t (v) stands for the time when the maximum value c (v) is reached. The use of two thirdorder polynomials allows the boundary conditions — (i) equal concentration and (ii) equal first derivative of the concentration at t A, t B and t (v) — to be fulfilled. Using these boundary conditions for the polynomials returns the following eight conditional equations for the coefficients a i and b i : pA jtA ¼ cðtA Þ;

dpA dt

j

¼

tA

pA jtðvÞ ¼ pB jtðvÞ ¼ cðvÞ ;

pB jtB ¼ cðtB Þ;

dpB dt

j

tB

cðtA þ Dt Þ  cðtA Þ Dt

dpA dt

¼

j

t ðvÞ

¼

dpB dt

j

t ð vÞ

¼0

cðtB Þ  cðtB  Dt Þ Dt

ð4Þ

Fig. 4. The third-order polynomials p A (t) (dash-dotted line) and p B (t) (dashed line) replace the AIF for the unreliable part between t A and t B. The maximum c (v) of the resulting curve and its position t (v) can be freely chosen without violating continuity and differentiability at t A and t B. For t b t A or t N t B (solid line), a measured AIF is chosen.

Thus, by using the postulation for continuity and differentiability at t A, t B and t (v), the number of free parameters is reduced to t (v) and c (v). Both need to be derived in the next and final step of the correction scheme. For this step, the correction algorithm uses a parametric convolution approach by defining a model for the tissue transport function, which is characterized by three parameters k 0, k 1 and k 2: sffiffiffiffiffiffiffiffiffiffiffiffiffiffi   k1 3 k1 2 t  k exp  ð Þ hðt Þ ¼ k0 : 1 2k2 pt 3 2k2 t

ð5Þ

The model is described and motivated in the Appendix A. It can be applied for a tracer transportation into a sufficiently large volume of tissue. Given the model for the AIF and the tissue transport function, Eq. (2) can now be specified as:   cðt Þ ¼ AIF t; t ðvÞ ; cðvÞ  h t; k0 ; k1 ; k2 :

ð

Þ

ð6Þ

Thus, the parametrical description of h(t) allows an estimation of the AIF by varying t (v), c (v), k 0, k 1 and k 2. The cost function of the fit is the difference between a reference curve selected in the tissue and its best prediction by Eq. (6). The cost landscape found when solving this problem is relatively flat and noisy. Therefore, a robust optimization procedure has to be chosen, which was realized by selecting the DIRECT algorithm [19]. As a result, all of the parameters needed to describe the unreliable part of the AIF (t (v) and c (v)) are obtained. Fig. 5 depicts the optimal solutions for the six exemplary volunteers. Such a compounded AIF (the parametric description is merged with the original data from the reliable parts) can now serve for any parametric and nonparametric deconvolution techniques to characterize the perfusion state.

1304

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311

3. Simulation The simulation which will now be introduced serves as a platform to test the correction algorithm. To evaluate the stability of the optimization step for a relevant setting, the AIF template and the noise characteristics were chosen close to in vivo measurements.

Fig. 6. (A) AIFs of 6 subjects from ROIs located at M3 segments of the MCA after applying a parametric superposing algorithm. (B) Same as (A) but using M1 segments of the MCA. (C) AIF used for the simulation consisting of two gamma-variate functions describing the first (thin solid line) and second pass (dotted line). The transition to the steady-state level is given by the dashed line. The parameters of the three processes are derived from in vivo measurements from the M1 segment of the MCA shown above (B). The scaling factor r was set to 1.0.

The simulation consists of three building blocks: (i) The AIF was derived from in vivo measurements, and we proposed a parametric description which included the first and second pass and the steady state; (ii) the saturation phenomena was included in the simulation via a Riciandistributed noise model, affecting the raw image signal S(t); and (iii) the residual function of brain tissue was taken into account by applying an exponential tissue model. To decide a systematic interaction between the exponential tissue model and the optimizing kernel (transport–diffusion), all calculations were also repeated with a box-shaped model. 3.1. Parametric template of the AIF

Fig. 5. The corrected AIFs of 6 subjects (bold line) compared to the individual uncorrected AIFs (thin lines).

For the simulation, an artificial AIF was prepared, derived from the data from the six volunteers. Therefore, the measured AIFs were iteratively normalized by adjusting each one to their mean circulation time (by scaling the time axis), mean temporal offset (by introducing a temporal shift) and mean signal intensity (by applying a scaling factor). Fig. 6A and B show the result of this procedure for both M3 and M1 levels of MCA bifurcation. Thereby, an AIF is defined as an average representative of the measured AIFs.

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311

For reasons of practicality and to allow an export to other simulations, a parametrical representation was used. Even though a single gamma-variate function is often applied for such a task [20], it can only account for the first pass. Indeed, the complete behavior can be described by a sum of three components: AIFðt Þ ¼ r½c1 ðt Þ þ c2 ðt Þ þ Hðt Þ:

ð7Þ

Here, c 1(t) and c 2(t) are gamma-variate functions for the first and second pass, and Q(t) is a transformed complementary error function (erfc) simulating the transition to the steady state. The parameter r allows for an additional scaling, which is subsequently used to mimic different partial volumes or tracer concentrations. By fitting Eq. (7) to the average AIF according to the M1 level of the MCA shown in Fig. 6B, the following parametric description was found:

1305

With S 0 being the signal amplitude in absence of noise, the following probability density of S can be derived:     S S 2 þ S02 SS0 pðS Þ ¼ 2 exp  I ; ð10Þ 0 r 2r2 r2 where I 0 is the modified Bessel function of the first kind and zeroth order. While a gaussian distribution can be used as a sufficient model for small r or large S 0, in the opposite case, p(S) is going asymmetric, and its expectation value differs from the position of the maximum. Generally, the former can be analytically derived from the expectation value E(S): rffiffiffiffi     2 p S0 2 S0 2 S E ðS Þ ¼ r exp  2 1 þ 2 I0 0 2 2 4r 2r 4r  2  S 2 S þ 0 2 I1 0 2 ; ð11Þ 2r 4r where I 1 is the modified Bessel function of the first kind and first order. Depending on SNR, this equation leads to a positive bias in observed signal magnitudes [22].

c1 ðt Þ ¼ 7:72d 104 d ðt  10:2Þ6:683   t  10:2  exp  ; if tN10:2; else ¼ 0 1:394 c2 ðt Þ ¼ 7:41d 1093 d ðt  1:6Þ79:36   t  1:6  exp  ; if tN1:6; else ¼ 0 0:4887   t  39:4 Hðt Þ ¼ 0:171d erf c  9:374

ð8Þ

Fig. 6C shows the parametric AIF by a thick line and the three components c 1(t), c 2(t) and Q(t) by thin lines. 3.2. Saturation of the AIF due to Rician-distributed magnitudes As indicated in Fig. 3, AIF measurements can be perturbed by saturation issues. This phenomenon can be explained by a nongaussian noise influencing the magnitude of reconstructed T 2*-weighted images. A saturation-like behavior at the peak concentration is often observed due to a low magnitude, which may result in low SNR [12]. The basis for the saturation is that an additive gaussian noise is found in the k space matrix, which originated from the receiver unit [21]. Due to the invariance of this noise for the Fourier transformation, both the real (S real) and the imaginary part (S imag) of the raw image matrix are polluted by two statistically independent sources of gaussian noise, n real and n imag, each with a variance of r 2. In this case, the image intensity S is mathematically equal to the signal magnitude qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

2ffi S ¼ ðSreal þ nreal Þ2 þ Simag þ nimag ; ð9Þ where a Rician distribution is expected [22].

Fig. 7. The noise model (Rician-distributed) changes the expectation value of the concentration time course, depending on the SNR and the relative tracer concentration r. (A) For a fixed SNR = 10, r was varied, and the resulting curves were scaled to the noiseless case over the steady state. (B) Same as (A) but now with fixed r = 1 and varying SNR. In both cases, a saturation-like behavior can be observed for the first pass. (C) Statistical fluctuations are large around the first pass, and they generate overshoots especially in the saturated part.

1306

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311

Table 1 The parameters used in the Monte Carlo simulation Parameter

Symbol

Unit

Value

Cerebral blood volume Cerebral blood flow Mean transit time Standard deviation of noise Signal-to-noise ratio Prebolus signal Relative bolus dosage Number of simulation runs Tissue-amplifying factor

rCBV rCBF MTT r SNR S(0) r N f tissue

ml/100 ml ml/100 ml/min s 1 1 1 1 1 1

2/3/5 20/45/100 6/4/3 0, 15, 30,. . ., 120 l, 40, 20,. . ., 5 600 0.5, 0.75, 1.0,. . ., 2.0 120 2.5

The flow parameters (line 1–3) were altered to characterize the different types of brain tissue (white matter/mixture/gray matter).

However, overestimation of magnitudes is inverted into an underestimation of tracer concentration by Eq. (1) due to the characteristic of the negative logarithm, which is monotonically decreasing. When the AIF is taken from large vessels, signal drops may lead to a low SNR, and the resulting erroneous values dominate T 2*-weighted DSCMRI. Fig. 7A and B demonstrate the expectation value of the artificial AIF for different noise levels and a variation of the tracer concentration. For an SNR greater then 5 and a relative tracer concentration r between 0.5 and 2, a saturation phenomenon is observed but restricted to the first pass. This is in line with findings in volunteers that saturation leading to deviations from linear model [Eq. (1)] is restricted to first pass only [12]. Fig. 7C shows the 95% confidential interval for a distinct parameter combination (SNR =10, r = 1.0). It is demonstrated that above the expectation values, fluctuations occur in a broader range then below it and can easily exceed it. Therefore, in a single trial, this asymmetrically distributed probability might cover the saturation phenomenon to the observer.

To validate possible interactions between the tissue response function for large ROIs [Eq. (5)] and the choice of the tissue model used here, a box-shaped model was also included in the analysis: Rbox ðt Þ ¼ rCBFð H ðt Þ  H ðMTT  t ÞÞ;

ð13Þ

where H(t) is the Heaviside step function. 3.4. Parameters used in the simulation For the evaluation of the algorithm, the simulation was performed to return the measurement curves of a DSC-MRI

3.3. Modeling the tissue function To simulate a DSC-MRI examination, a model for the residue function R(t) is needed. We used a commonly applied model, which is different from the transport– diffusion model described above: Rexp ðt Þ ¼ rCBFd expð  t=MTTÞ

ð12Þ

This response function was stepwise linearly convolved with the artificial AIF by an algebraic convolution scheme [5]. In addition, an amplifying factor f tissue was set to 2.5 to correct for the small partial volume of arteries. Since the correction algorithm needs AIFs from different arterial branches and therefore different partial volumes, a Gaussdistributed factor K (mean =1.0, variance =0.04) was introduced post scaling the partial volume of every individual arterial branch. In a final step, all generated data (AIFs and simulated tissue) are calculated from concentrations into image magnitudes, consecutively polluted with Riciandistributed noise of different strength and, finally, will be back-computed into relative concentrations [c(t)].

Fig. 8. rCBV, rCBF and MTT determined using a truncated SVD for two different arterial signal strengths r. An SNR of 10 was assumed. The thin line indicates the value obtained in absence of noise. A single AIF (S) and a corrected AIF (C) is compared. Significant deviations from the expected results are depicted with asterisks. *P = .05; **P = .005.

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311

experiment. First, noise-free AIFs where generated using Eq. (7). These AIFs where convoluted with the tissue response function to allow for a calculation of the signal time curves in the tissue [Eqs. (12), (13)]. Furthermore, a measured AIF was simulated by disturbing the noise-free AIF with the noise model [Eq. (10)]. The main input of the simulation were the perfusion parameters representing blood flow in gray matter, a mixed compartment and white matter (Table 1). Furthermore, parameters such as the number of runs and the prebolus signal had to be fixed. For each parameter combination of SNR and r signal-time curves were generated. To be in line with common practice in DSC-MRI, a 33 uniform filter was simulated by calculating and consecutively averaging nine voxels for each run and parameter combination. The six AIFs were generated by simulation of individual regions, which differed in their noise patterns and in their partial volumes. In addition, as for a realistic measurement, each AIF was derived from five simulated voxels modeling an arterial input ROI. Finally, to simulate the reference region consisting of a mixture of different tissues types, 100 voxels were generated with rCBV varying from 2 to 5 ml/ 100 ml and MTT varying from 3 to 6 s. For the final evaluation described in the next section, the results of the simulations were fed to the correction algorithm, which returns a corrected AIF. To be able to quantify the effect on the determined perfusion parameter, the concentration time course in the tissue and the AIF with

1307

and without corrections were analyzed. This was done by performing a deconvolution analysis based on a truncated singular value decomposition (SVD, threshold = 0.20). 4. Results The simulation data were analyzed exactly as it would have been done for a real experiment. The focus of the evaluation was the comparison between the perturbed and the corrected AIF. Fig. 8 shows rCBV (upper plots), rCBF (middle plots) and MTT (lower plots) determined from simulated data. Two different relative tracer dosages (r = 1.0, 2.0) were applied. The main comparison was performed by using the AIF from a single arterial branch (bSQ) vs. applying the proposed correction algorithm (bCQ). Independent of the exact parameterization, some trends were observable: for single AIFs, rCBV and rCBF are always overestimated in comparison with the expectation value (thin horizontal line: MTT = 4 s, rCBV= 3 ml/100 ml, rCBF = 45 ml/100 ml/min), whereas MTT was underestimated in all cases. The overestimation became more relevant for rCBF and rCBV if the relative concentration is r = 2.0. This is more likely to happen in branches larger than the MCA at M1 level or by using a higher tracer dosage. Especially, rCBF and rCBV were significantly overestimated with P b.005 and .05, respectively. MTT reductions were in the same order of magnitude in both cases. For all cases, the results generated by using a corrected version

Fig. 9. Relative alteration in rCBV (A and B) and rCBF (C and D) caused by the saturation effect for different SNR ( y-axis) and r (x-axis). An MTT of 4 s, rCBV of 3 ml/100 ml and rCBF of 45 ml/100 ml/min was assumed. (A and C) In unmodified AIFs, strong deviations of 50% and more can be observed for large values of r or small values of SNR. (B and D) Using corrected AIFs, only small differences to the expected values were found.

1308

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311

Table 2 Original and corrected values of rCBF, rCBV and MTT in cortical gray matter for the six volunteers Subject 1 2 3 4 5 6 MeanFS.E.M.

rCBF (a.u.)

rCBV (a.u.)

MTT (s)

Original

Corrected

Change

Original

Corrected

Change

Original

Corrected

Change

0.073 0.108 0.145 0.188 0.164 0.079 0.126F0.019

0.056 0.062 0.118 0.129 0.142 0.064 0.095F0.016

23% 43% 19% 32% 13% 19% 25F4%*

0.294 0.454 0.509 0.591 0.657 0.329 0.472F0.058

0.254 0.296 0.413 0.480 0.570 0.299 0.385F0.051

14% 35% 19% 19% 13% 9% 18F4%*

4.06 4.32 3.57 3.15 3.87 3.98 3.82F0.17 s

4.53 4.96 3.60 3.75 3.99 4.58 4.24F0.22 s

+12% +15% +1% +19% +3% +15% +11F3%*

Percentage change caused by the correction is also shown. A significant effect by the correction scheme ( P b.05) is marked with an asterisk. a.u., arbitrary unit.

of the AIF were very close to the expected values, and the deviations never reached statistical significance. Finally, the standard deviation of the corrected perfusion parameters were almost smaller than for single AIFs. Fig. 9 shows the deviations to rCBV (upper plots, A and B) and rCBF (lower plots, C and D) caused by the saturation phenomena in the uncorrected measurements (left plots, A and C) and after using the correction scheme (right plots, B and D) over a wide range of r and r values. The tissue parameters used here are MTT= 4 s, rCBV=3 ml/100 ml, rCBF =45 ml/100 ml/min. In the native, uncorrected case large deviations can be found for certain parameter combinations. Especially, large values of r or small values of SNR induce a severe misjudgment of perfusion parameters. Using corrected AIFs, deviations are less marked, and even for adverse combinations of r and r (both r and r are high), smaller differences (rCBV: b 10%; rCBF: b20%) are found with respect to the observed deviations without AIF correction. Note that for SNR values greater than 10 in the corrected case deviations were smaller then 5% and were nearly independent from the relative concentration r. The most general observation is, however, that the correction algorithm strongly increases the span of SNR and concentration values, which can be applied without a negative influence on the perfusion parameters. This observation was valid for the parameters summarized in Table 1. The analysis shown in Figs. 8 and 9 was also performed with the box-shaped tissue function [Eq. (13)], and very similar results were obtained. Table 2 shows the results after applying the correction scheme for the in vivo data for the six volunteers. rCBF, rCBV and MTT are presented for an analysis with (bcorrectedQ) and without applying (boriginalQ) the correction scheme. Note that the correction leads to a significant alteration in percentage change of rCBF, rCBV and MTT. 5. Discussion We were able to show that (a) an arterial input function perturbed by Rician-distributed measurement noise leads to an overestimation of rCBV and rCBF and that (b) a novel correction algorithm can compensate for these artifacts. In the following, the finding on the perturbation of the AIF and the major mechanism of the correction algorithm, the

derivation of the peak of the AIF from the low concentration phase, is discussed. 5.1. Noise leads to saturation In the present work, the saturation of the AIF is modeled by applying Rician-distributed noise, which is a more realistic statistical description of magnitudes in MR images, as compared to the gaussian one [22]. Fig. 7A and B show that this saturation depends on the SNR as well as on the relative tracer concentration r for an artificial arterial input function. A saturation-like behavior is observed for large values of r and small values of SNR. Moreover, Fig. 7C showed that the variance rises dramatically. The SNR evaluated here is lower than the global SNR normally determined for the T 2*-weighted images. However, in the arteries, physiological noise (e.g., pulsation) in combination with noise from the MR receiver unit are usually much larger then in other regions of the brain, and therefore, low SNR values are expected. For this reason, the SNR in the arterial input voxel is a major issue and is much lower than the global SNR values, which are commonly determined using noise from an empty space outside of the skull. An ideal arterial input voxel tries to optimize two contradicting principles: on the one hand, the voxel should include as much arterial volume as possible to decrease the effect of a capillary contamination and partial volume effects [11,23]; on the other hand, the high concentration of contrast agent in the arteries induces the saturation issues discussed here [12,24]. So far, three strategies have been proposed to reduce the later effect. The first strategy was shown by van Osch et al. [23]: in the selected case of arteries orientated perpendicular to the imaging plane, the distortions in magnitude and phase can be calculated, consecutively corrected and used for a correction scheme. In this model, low concentrations showed only weak dependencies on the resulting partial volume effects. A second strategy proposes to use optimal dosages of the contrast agent [24]. Finally, further strategies for selecting appropriate voxels for AIF selection have been suggested [25]. Different from these approaches, our correction scheme can be applied for any choice of the arterial voxel and without a reduction of the dosage of the contrast agent, potentially resulting in SNR issues in the parenchyma.

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311

5.2. Deducing the peak from lower concentrations We postulated that the peak of the AIF can be deduced from the reliable part of the AIF by the proposed correction algorithm. The fact that the saturation phenomena affects only higher concentrations around the first pass is in line with the findings of other groups [13,14,25]. However, in the case of an immense tracer concentration, the second pass and steady state could also show saturation. At the same time, the first pass is most crucial for the success of deconvolution strategies like SVD [26]. The presented approach tries to bypass these problems by reconstructing the first pass of the AIF based on information from nonaffected, reliable parts of the AIF (second pass and steady state). To estimate the unperturbed AIF, a model was needed that could describe the transition of the paramagnetic tracer from the arterial part of the vascular bed to the capillaries. A one-dimensional form of the transport–diffusion equation was chosen for this task. Due to the statistical nature of this equation, it is well suited for an averaged concentration-time course in a large reference region. The unreliable part of the AIF was replaced by a parametric description. There, two polynomials of third degree were utilized. An alternative would be to use other classes of functions. However, this would involve more control parameters potentially overstressing the overall fitting procedure. By varying the polynomials, the relationship between the real concentration time course in the tissue and the one predicted by the transport–diffusion model was compared. In the case of an optimal prediction, the parametrically described first pass was concluded as the most probable. This correction approach should be able to handle saturation issues observable in arteries. We do not solve the issue of an absolute quantification. By itself, DSC-MRI is not quantitative, since the individual arterial tree has a large influence on the resulting flow parameters. This is demonstrated in Table 2, where the perfusion parameters of the six volunteers are shown. Despite a thorough selection of the arterial voxels, deviations in rCBV of more then 100% between individuals can be observed. Nevertheless, our approach has the potential to improve longitudinal studies, since it reduces the variance caused by the saturation of the AIF. In combination with a careful selection of arteries, this could be a successful strategy to improve reliability. To apply this correction scheme in cases of pathological perfusion (tumor, stroke), the user only has to ensure that the arterial signal and the reference tissue show a nonpathological pattern. In the case of ischemia, for instance, the nonaffected hemisphere could be used to obtain these inputs for the algorithm. The approach suggested in this work was not restricted to saturation issues. We propose that it is suitable for all cases where the bhigh concentration sectionQ of the arterial input function has to be discarded, whereas the blow concentration sectionQ of the curve is btrustable.Q This may become

1309

relevant, e.g., when introducing modern contrast agents with a high magnetization and, thus, a strong signal drop or using an improved spatial resolution leading to lower SNR. 6. Conclusions We demonstrate that the first pass of the arterial input function in DSC-MRI can be completely reconstructed from its later parts, with lower tracer concentration and a reference region in the tissue. This observation is relevant for all cases where the definition of the first pass is uncertain, e.g., when SNR is low or due to saturation effects. The reconstruction is achieved by applying a transport–diffusion equation in combination with an optimization procedure. We found good agreement with the expected perfusion parameters in Monte Carlo simulations. Appendix A A.1. Parametric model of blood transportation In this section, the parametric model describing the distribution of an intravascular tracer along a vascular network will be presented, which can be used to correct the peak saturation of the AIF as proposed in the last section. The transport approximation was used since the average tracer curve in a relatively large reference region had to modeled. Regarding perfusion as a process connected with blood transport and dispersion along the vascular system, a Navier–Stokes equation can be applied, which describes the three-dimensional fluid dynamics of incompressible media [27]. With some assumptions, the complexity of the three-dimensional equation can be reduced to the onedimensional random walk model with drift [28] Bc B2 c Bc ¼D 2 u Bt Bx Bx

ð14Þ

The problem is thus reduced to a second-order, onedimensional partial differential equation, thus not relying on

Fig. 10. The parametric response function h(t) for four different values of k 2 (s2), with k 0 = 1 and k 1 = 1 s. The function models a washout process for k 2 N k 1 (black line), a transport-dominated process for k 2 b b k 1 (gray line) or an intermediate behavior (dashed and dotted lines).

1310

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311

any knowledge about vascular architecture. Here, x maps the complex path of the tracer towards a sample volume onto a one-dimensional real axis. D and u are parameters characterizing blood transportation, where u is related to a convective drift dominantly appearing in the larger vessels, and D is linked to the dispersion process in the capillary system as well as from different transport paths towards the same observation volume. For an initial condition expressed by an arriving arterial bolus ( = AIF) associated with the boundary condition c(x = 0,t) = AIF(t) and c(xYl,t)u0, Eq. (14) can be solved to yield   k0 x ð x  ut Þ2 cð x; t Þ ¼ AIFðt Þ  pffiffiffiffiffiffiffiffiffiffiffiffi exp  4Dt 4pDt 3

[5]

[6]

[7]

[8]

In difference to a previously presented solution [29], k 0 is an additional free parameter, which is identified with the partial blood volume in the ROI and is equal to the zeroth moment of the response function h(t). Thus, k 0 is identical to the cerebral blood volume rCBV [30]. With the first moment k 1 and second centralized moment k 2, Eq. (15) and, hence, its higher-order moments, can be completely characterized: k0 ¼

[3] [4]

ð15Þ

¼ AIFðt Þ  hðt; k0 ; x; D; uÞ:

Z

[2]

[9]

[10]

[11]

l

hðt Þdt 0

k1 ¼

Z

[12]

l

thðt Þdt 0

k2 ¼

Z

[13]

l

ðt  k1 Þ2 hðt Þdt:

ð16Þ

0

[14]

Expressing h(t) in terms k 0, k 1 and k 2 results in [15]

sffiffiffiffiffiffiffiffiffiffiffiffiffiffi   k13 k1 2 exp  ð t  k1 Þ hð t Þ ¼ k0 2k2 pt 3 2k2 t

ð17Þ

Fig. 10 shows h(t) for four different parameter combinations. In cases of k 2 b b k 1, h(t) is dominated by a pure transport phenomenon, resulting in a narrow curve with a temporal offset. For large k 2 (with regard to k 1), the perfusion is mainly modeled by a dispersion process (washout curve).

[16]

[17]

[18]

[19]

References [1] Villringer A, Rosen BR, Belliveau JW, Ackerman JL, Lauffer RB, Buxton RB, et al. Dynamic imaging with lanthanide in normal brain:

[20]

contrast due to magnetic susceptibility effects. Magn Reson Med 1988;6:164 – 74. Li T, Chen ZG, astergaard L, Hindrnarsh T, Moseley ME. Quantification of cerebral blood flow by bolus tracking and artery spin tagging methods. Magn Reson Imaging 2000;18:503 – 12. Warach S. Stroke neuroimaging. Stroke 2003;34(2):345 – 7. Barbier EL, Lamalle L, De´corps M. Methodology of brain perfusion imaging. J Magn Reson Imaging 2001;13:496 – 520. astergaard L, Weisskoff RM, Chesler DA, Gyldensted C, Rosen BR. High resolution measurement of cerebral blood flow using intravascular tracer bolus passages. Part I: Mathematical approach and statistical analysis. Magn Reson Med 1996;36:715 – 25. Andersen IK, Szymkowiak A, Rasmussen CE, Hanson LG, Marstrand JR, Larsson HBW, et al. Perfusion quantification using gaussian process deconvolution. Magn Reson Med 2002;48:351 – 61. Wu O, astergaard L, Weisskoff RM, Benner T, Rosen BR, Sorensen AG. Tracer arrival timing-insensitive technique for estimating flow in MR perfusion-weighted imaging using singular value decomposition with a block-circulant deconvolution matrix. Magn Reson Med 2003;50:164 – 74. Calamante F, Gadian DG, Connelly A. Quantification of bolustracking MRI: improved characterization of the tissue residue function using Tikhonov regularization. Magn Reson Med 2003; 50:1237 – 47. Yang C, Karczmar GS, Medved M, Stadler WM. Estimating the arterial input function using two reference tissues in dynamic contrastenhanced MRI studies: fundamental concepts and simulations. Magn Reson Med 2004;52(5):1110 – 7. Yankeelov TE, Cron GO, Addison CL, Wallace JC, Wilkins RC, Pappas BA, et al. Comparison of a reference region model with direct measurement of an AIF in the analysis of DCE-MRI data. Magn Reson Med 2007;57(2):353 – 61. Thornton RJ, Jones JY, Wang ZJ. Correcting the effects of background microcirculation in the measurement of arterial input functions using dynamic susceptibility contrast MRI of the brain. Magn Reson Imaging 2006;24:619 – 23. Ellinger R, Kremser C, Schocke MFH, Kolbitsch C, Griebel J, Felber SR, et al. The impact of peak saturation of the arterial input function on quantitative evaluation of dynamic susceptibility contrast-enhanced MR studies. J Comput Assist Tomogr 2000; 24(6):942 – 8. Chen JJ, Smith MR, Frayne R. The impact of partial-volume effects in dynamic susceptibility contrast magnetic resonance perfusion imaging. J Magn Reson Imaging 2005;22:390 – 9. van Osch MJ, van der GJ, Bakker CJ. Partial volume effects on arterial input functions: shape and amplitude distortions and their correction. J Magn Reson Imaging 2005;22(6):704 – 9. Kiselev VG. On the theoretical basis of perfusion measurements by dynamic susceptibility contrast MRI. Magn Reson Med 2001;46: 1113 – 22. Vonken E, van Ossch MJP, Bakker CJG, Viergever MA. Measurement of cerebral perfusion with dual-echo multi-slice quantitative dynamic susceptibility contrast MRI. J Magn Reson Imaging 1999; 10:109 – 17. Rausch M, Scheffler K, Rudin M, Radq EW. Analysis of input functions from different arterial branches with gamma variate functions and cluster analysis for quantitative blood volume measurements. Magn Reson Imaging 2000;18:1235 – 43. Sijbers J, den Dekker AJ, van Audekerke J, Verhoye M, van Dyck D. Estimation of the noise in magnitude MR images. Magn Reson Imaging 1998;16(1):87 – 90. Gablonsky JM, Kelley CT. A locally-biased form of the DIRECT algorithm. J Global Optim 2001;21(1):27 – 37. Thompson HK, Starmer CF, Whalen RE, McIntosh HD. Indicator transit time considered as a gamma variate. Circ Res 1964;14: 502 – 15.

P. Brunecker et al. / Magnetic Resonance Imaging 25 (2007) 1300 – 1311 [21] Rice SO. Mathematical analysis of random noise. In: Wax N, editor. Selected papers on noise and stochastic processes. New York7 Dover Publications; 1958. p. 133 – 294. [22] Sijbers J, den Dekker AJ, Raman E, van Dyck D. Parameter estimation from magnitude MR images. Int J Imaging Syst Technol 1999;10:109 – 14. [23] van Osch MJP, Vonken E, Viergever MA, van der Grond J, Bakker CJG. Measuring the arterial input function with gradient echo sequences. Magn Reson Med 2003;49:1067 – 76. [24] Smith MR, Lu H, Frayne R. Signal-to-noise ratio effects in quantitative cerebral perfusion using dynamic susceptibility contrast agents. Magn Reson Med 2003;49:122 – 8. [25] Duhamel G, Schlaug G, Alsop DC. Measurement of arterial input functions for dynamic susceptibility contrast magnetic resonance imaging using echoplanar images: comparison of physical simulations with in vivo results. Magn Reson Med 2006;55(3):514 – 23.

1311

[26] Murase K, Shinohara M, Yamazaki Y. Accuracy of deconvolution analysis based on singular value decomposition for quantification of cerebral blood flow using dynamic susceptibility contrastenhanced magnetic resonance imaging. Phys Med Biol 2001; 46:3147 – 59. [27] Calamante F, Yim PJ, Cebral JR. Estimation of bolus dispersion effects in perfusion MRI using image-based computational fluid dynamics. Neuroimage 2003;19(2 Pt 1):341 – 53. [28] Sheppard CW, Uffer MB. Stochastic models for tracer experiments in the circulation, II. Serial random walks. J Theor Biol 1969;22: 188 – 207. [29] Sheppard CW. Stochastic models for tracer experiments in the circulation: parallel random walks. J Theor Biol 1962;2: 33 – 47. [30] Zierler KL. Theoretical basis of indicator-dilution methods for measuring flow and volume. Circ Res 1962;10:393 – 407.