COMPUTER PROCESSING AND MODELING Full Paper
Magnetic Resonance in Medicine 74:280–290 (2015)
Correcting for Large Vessel Contamination in Dynamic Susceptibility Contrast Perfusion MRI by Extension to a Physiological Model of the Vasculature Michael A. Chappell,1* Amit Mehndiratta,1 and Fernando Calamante2,3 Purpose: Dynamic susceptibility contrast (DSC) perfusion images are contaminated by contributions from macro vascular signal arising from contrast agent within the larger arteries that do not contribute directly to the local tissue perfusion. Methods: A vascular model of the DSC perfusion signal was extended by the inclusion of a macro vascular component based on the arterial input function. This was implemented within a Bayesian nonlinear model-fitting algorithm that included automatic model complexity reduction. Results were compared with existing methods that do not correct for the macro vascular contamination as well as an independent component analysis technique. Results: Macro vascular signal was identified in regions corresponding to larger arteries resulting in reductions by 62% within a region of interest identified with high contamination. Whereas visually similar results could be achieved with independent component analysis, it resulted in reductions in global tissue perfusion and was not robustly applicable to patient data. Conclusion: A model-based strategy for correction of macro vascular contamination in DSC perfusion images is feasible, although the model may currently need extending to more accurately account for nonlinear effects of contrast agent in C 2014 large arteries. Magn Reson Med 74:280–290, 2015. V Wiley Periodicals, Inc. Key words: cerebral blood flow (CBF); dynamic susceptibility contrast MRI; vascular model; Bayesian analysis; macro vascular contamination
INTRODUCTION The quantification of dynamic susceptibility contrast (DSC) perfusion-weighted MRI typically relies upon the identification of an arterial input function (AIF) and subsequent voxelwise numerical deconvolution to calculate 1 Institute of Biomedical Engineering, University of Oxford, ORCRB, Old Road Campus, Headington, Oxford, United Kingdom. 2 Florey Institute of Neuroscience and Mental Health, Heidelberg, Australia. 3 Department of Medicine, Austin Health and Northern Health, University of Melbourne, Melbourne, Australia. Grant sponsor: the Wellcome Trust and EPSRC; Grant number: WT088877/ Z/09/Z. *Correspondence to: Michael A. Chappell, D.Phil., Institute of Biomedical Engineering, University of Oxford, ORCRB, Old Road Campus, Headington, Oxford, OX3 7DQ, UK. E-mail: [email protected] Received 10 June 2014; revised 4 July 2014; accepted 7 July 2014 DOI 10.1002/mrm.25390 Published online 8 August 2014 in Wiley Online Library (wileyonlinelibrary. com). C 2014 Wiley Periodicals, Inc. V
perfusion, mean transit time (MTT) or other quantitative parameters (1) A fundamental assumption in DSC-MRI is that the signal in the voxel arises from the Gadolinium contrast agent within the micro vasculature (primarily capillaries, but also arterioles and venules). Thus the residue function estimated by means of deconvolution represents the distribution of transit times within the voxel and the MTT is the mean of this distribution. However, a proportion of voxels will also include a contribution from contrast agent within larger vessel, for example those proximal to the middle cerebral arteries, either because the vessel is within the voxel or due to susceptibility effects of a vessel nearby. Therefore, this extra contribution results in errors to CBF, MTT, and cerebral blood volume (CBV) quantification, because the measured signal does not reflect solely micro vasculature (2). Despite this, the presence of this macro vascular (MV) contribution is typically ignored. Various methods have been proposed to deal with this “unwanted” MV contribution (NB. in fact, its presence is critical to the ability to extract AIF from the data in the first place). For example, Carroll et al (3) proposed that independent component analysis (ICA) could be used as a preprocessing step to identify and remove components from the signal that could be regarded as macro vascular. Subsequently, Reishofer et al (4) extended this approach, by proposing an automated scheme to identify the relevant MV components. Such an approach is however limited by how accurately appropriate MV components can be identified in practice, particularly in patients with severe vascular abnormalities (5). In principle, the contribution of the MV component to the signal should match that of the AIF, ignoring dispersive effects (6) and complex partial volume effects (7); thus, it should be possible to use the existing information about the measured AIF to assist in MV correction. In this work, we propose a model-based solution to the correction for MV contamination inspired by related recent work in Arterial Spin Labeling (8), building upon work previously presented in abstract form (9,10). To do so, we build upon the vascular model of Mouridsen et al (11) for the DSC-MRI tissue signal, by incorporating a separate MV component, which must be estimated simultaneously from the DSC-MRI data. Like the similar methodology used in Arterial Spin Labeling (8), we subject the extra MV component to an automated model complexity reduction procedure to ensure that the extra
280
Correcting Contamination in DSC Perfusion
281
complexity introduced by adding this component is only used where the data justifies it. Thus, avoiding over fitting of the data in regions where no large vessels are present. The method was tested on in vivo data, and compared with the ICA-based approach. THEORY
mate and apply it to the tissue, because there will be some contribution from the macro vascular blood signal. This can be addressed by a modified version of the vascular model, termed hereafter the mVM. In this model, MTT becomes a further parameter to be estimated from the data within the model fitting, with b ¼ MTT/l. A total of four parameters have be estimated in the mVM: rCBF, k, and d and MTT.
Vascular Model The concentration of contrast agent within a voxel, C(t), can be described according to tracer kinetics as the convolution of an AIF, Ca(t), and a tissue residue function, R(t), scaled by the CBF (1): CðtÞ ¼ CBF ðCa ðtÞ RðtÞÞ¼CBF
Z
t
Ca ðtÞ RðttÞdt
[1]
0
In the vascular model of Mouridsen et al (11), the residue function was parameterized using a model of contrast agent transit through the capillary network that followed a gamma distribution: RðtÞ ¼
Z
1
hðt; a; bÞdt t
hðt; a; bÞ ¼
t 1 t a1 exp b ; a; b > 0 ba GðaÞ
[2]
where a and b are the shape and scale parameters of the distribution respectively, which provides a flexible model to describe a large range of feasible residue function shapes. Because the product ab is the mean of the distribution, this equates to the MTT for contrast agent transit. Mouridsen et al (11) used the central volume principle to reparameterize the distribution in terms of relative CBF (rCBF), relative CBV (rCBV) and a shape parameter: a ¼ l, b ¼ rCBV/(k .CBV). A further parameter d was included to account for differences in arrival between the AIF and the voxel concentration time course, i.e., bolus delay. The estimated concentration time course was converted to a predicted signal assuming that the transverse relaxation rate changes linearly with concentration of the contrast agent SðtÞ ¼ Sð0Þer2 CðtÞTE , where S(0) is the baseline signal before bolus arrival and r2 was incorporated into the CBF parameter so that relative CBF was estimated. Thus there were three parameters in the model: rCBF, k and d. These were estimated within a Bayesian nonlinear model-fitting algorithm, where the relative CBV, to be used in the expression for the shape parameter b, had previously been estimated from the data in the normal way by taking the ratio of the area under the concentration time curve and arterial input function. Priors were applied to the parameters using the estimates from a numerical deconvolution of the data, using the singular value decomposition (SVD) approach of (12), to set the mean for rCBF. A potential disadvantage of the original VM implementation was that the model fitting was constrained by the rCBV value obtained from the data, which due to noise could introduce a bias into the model fitting. A further issue arises in the context of modeling both a tissue and macro vascular component within a single voxel, as is the aim of this work. Under these conditions, it is no longer possible to take the standard rCBV esti-
Macro Vascular Component The presence of a macro vascular contribution within the voxel can be modeled by the addition of an extra component with a concentration time course that matches that of the AIF. This assumes that the flow of the contrast through the voxel is sufficiently rapid that there is no accumulation of the agent within the voxel, unlike the capillary signal. This is expected to be true of larger arteries and arterioles (8). This also neglects any dispersion effects as well as interactions between tissue and MV signals (see discussion). The macro vascular model can parameterized by a magnitude that represents the relative arterial blood volume (raBV) and an arrival time, da, that accounts for delay between the AIF measured elsewhere in the brain and its observation as a macro vascular signal in the voxel. The blood volume parameter relates to the magnitude of the AIF, because this is used as the basis for this model component, given that the scaling of this to absolute blood volume is generally not known only a relative blood volume can be estimated. The arrival time need not match that of the tissue if, for example, the blood in the artery were destined for elsewhere in the brain. This extension to the modified version of the vascular model is termed hereafter the mVMþMV, and can be used to correct for MV contamination. Two variants for the inclusion of the MV component in mVMþMV were considered, these are illustrated in Figure 1: 1. Combining the mVM and MV contributions in the concentration time course: CðtÞ ¼ rCBF ðCa ðtÞ RðtÞÞ þ raBV Ca ðtÞ SðtÞ ¼ Sð0Þer2 CðtÞTE
[3]
2. Combining the corresponding individual signal time courses: CðtÞ ¼ rCBF ðCa ðtÞ RðtÞÞ SðtÞ ¼ Sð0Þðer2 CðtÞTE þ er2 raBV CaðtÞTE Þ
[4]
Note that, whereas the rCBF parameter values might be comparable between the two models, the same may not be true of the raBV parameters. METHODS Model Implementation The mVM and mVMþMV were implemented within a Bayesian nonlinear model fitting algorithm (13), this method differs from that originally used by Mouridsen et al (11), but similarly permits prior information to be
282
Chappell et al.
FIG. 1. Illustration of the mVMþVM model showing both the summation of the concentration time curves (a) and the summation of the predicted signals (b) arising from the separate tissue and MV concentration time curves.
incorporated within the analysis. An implementation of the methods used here can be found in the FMRIB Software Library (www.fmrib.ox.ac.uk/fsl/verbena). The priors used for the three models are given in Table 1. Priors for the mVM were chosen to either be noninformative (for the rCBF parameter) or based on physiologically expected values; unlike the VM, no prior information was derived from the data themselves. For the MV component an Automatic Relevancy Determination (or shrinkage) prior was placed on the raBV parameter (14), the role of this prior was to suppress the MV component in voxels where the data did not support its inclusion. This reflects the expectation that MV signal will not be present in all voxels in the brain and thus this component should be removed where it is not supported to avoid over fitting; for more information see Chappell et al (8,13). For comparison, data were also analyzed using a more commonly used SVD-based approach (using the delay insensitive oSVD method from Wu et al) (12). The results
Table 1 Priors Specified for the Model Fitting Given as the Mean and Standard Deviation of a Normal Distribution Parameter rCBF d (sec) log(k) log(MTT) raBV da (sec) S(0)
mVM
mVMþMV
0, 106 0, 5 2.3,1 1.5,0.3 – – 0,106
0, 106 0, 5 2.3,1 2.3,0.3 0, 106 0,5 0,106
from the SVD analysis were also used to initialize the parameters in the VM and mVM model fitting, namely estimated rCBF and MTT were taken from the SVD analysis and used as the initial values for these parameters. Additionally, the S(0) parameter was initialized using the first volume in the data series. For the VM the SVD results were also used in the specification of the priors as in the original implementation. ICA was also used in an attempt to remove MV contamination. ICA was performed using MELODIC (15) from the FMRIB software library (www.fmrib.ox.ac.ul/fsl), and the MV components were identified manually and removed from the data before oSVD analysis. An important confound for quantification of perfusion using DSC is the dispersion of the arterial bolus. To investigate how this might impact the macro vascular contamination correction proposed here we repeated the model-based analysis above including dispersion of the AIF. The effects of dispersion were modeled by the numerical convolution of the data-derived AIF with a vascular transport function (VTF) (16,17). We used a gamma-variate VTF and chose the parameters to simulate a “low” amount of dispersion following (18) (namely shape parameters: p ¼ 1 and s ¼ 2 in Equation 3 of Mehndiratta et al) (18). The dispersion was taken as fixed and no attempt was made to estimate it directly from the data (see the Discussion section). Vascular Model Comparisons: Numerical Simulations First, the mVM was evaluated in comparison to the original VM and oSVD in the absence of a macro vascular component to identify if any errors might be introduced into the analysis by the use of the modified version.
Correcting Contamination in DSC Perfusion
283
FIG. 2. Estimated CBF values from the simulations for selected cases of residue function shape, noise and delay [c.f. (11) Figure 5], showing mean estimated CBF versus true value with error bars showing standard deviation across the separate realizations of the noise.
Comparisons were performed using the same simulation method as in Mouridsen et al (11). In brief: an AIF was simulated using an analytical gamma-variate function. Simulated data were generated for validation of the method with a CBV of 4 mL/100 g (typical gray matter), CBF in the range 10–70 mL/100 g/min in increment of 10 mL/100 g/min, delay of 0 and 65 s, and two residue
functions: exponential & box (approximated using k ¼ 100). Simulated signals were generated from concentration curves with baseline signal S(0) ¼ 100, echo-time TE ¼ 65 ms, the proportionality constant k was chosen such that 40% peak signal drop was attained at a CBF of 60 mL/100 g/min with CBV ¼ 4%, and 60% peak signal drop for the AIF signal. Zero mean Gaussian noise was
284
Chappell et al.
Table 2 Mean Ratios between Estimated and True CBF Values from the Simulations (Mean 6 Standard Deviation) SNR¼100
Exponential Exponential, delay¼5 Box-car
SNR¼20
mVM
SVD
mVM
SVD
1.0260.15 0.9460.24 1.0860.13
0.8660.11 0.8560.11 1.1560.09
0.9360.43 0.7960.37 1.0960.31
0.7560.22 0.8660.25 1.0560.30
added to the signals and AIF to achieve a given SNR, with both SNR of 100 and 20 being simulated. For each combination of CBF, residue function, and delay a total of 100 noisy signals and AIF were generated. The methods were compared by calculating the regression coefficient between the estimated rCBF and the true value over the range considered. In Vivo Data The MV correction method was evaluated on data from two patients. Pat1: left carotid stenosis (scanned postcarotid endarterectomy); Pat2: stroke patient with right middle cerebral artery (MCA) occlusion (scanned 9 hrs postsymptoms onset). DSC-MRI data were acquired at 3 Tesla (T) (Pat1) or at 1.5T (Pat2), using a gradient-echo EPI, TE/TR ¼ 20/1250 ms (Pat1) or 60/1740 ms (Pat2), and a single (Pat1) or double (Pat2) dose of contrast
agent. The AIF was measured on a branch of the contralateral MCA using the semiautomatic method in Penguin (www.cfin.au.dk/software/penguin). Macro Vascular Correction Comparisons To examine the presence of MV contamination in the estimated rCBF images, the ratio of rCBF (CBFratio) to that in a region of interest (ROI) chosen within the deep white matter (without MV effects) was calculated. To quantify the changes in rCBF and parameters of the VM, two further ROIs were defined that prospectively represented areas of substantial and negligible MV contamination, respectively. These ROI were defined as voxels whose CBFratio from the VM analysis was more than 4.5 (high contamination) and less than 1.5 (low contamination). Within these ROIs, the mean change in CBFratio, l, and MTT were calculated.
FIG. 3. Results for CBF quantification from pat1 showing rCBF and CBFratio from all analysis methods considered, and raBV in the mVMþMV methods.
Correcting Contamination in DSC Perfusion
285
FIG. 4. Results for CBF quantification from pat2 showing rCBF and CBFratio from all analysis methods except ICA correction, and raBV in the mVMþMV methods.
RESULTS Vascular Model Comparisons: Numerical Simulations Figure 2 shows estimated rCBF values plotted against the true value across the simulated data. Table 2 gives the mean ratios between estimated and true CBF along with the standard deviation across the different noise realizations. The mVM was generally found to produce more accurate CBF estimates than SVD with wider standard deviation. Comparing to the VM (under a similar simulation protocol) in Mouridsen et al (11) and their Figure 5 and Table 1, very similar performance was observed for the case without bolus delay, including overestimation of CBF for a box-car residue function shape. The mVM appeared to have higher accuracy for data subject to delay. Macro Vascular Correction Comparisons Figure 3 shows the rCBF and CBFratio images in pat1 for all the methods considered, along with raBV from the mVMþMV approaches. Figure 4 shows the same information for pat2, except for the results from ICA correction because it was not possible to identify ICA components that were suitably representative of MV signal in this case. There was a substantial reduction in rCBF and CBFratio observed in regions in which large
arteries would be expected in all of the MV correction approaches. Some residual regions of high CBFratio persisted in the mVMþMV methods. The raBV images had high magnitude in regions consistent with the presence of MCA branches. Regions of high intensity of raBV were more extensive in the mVMþMV method using the summation of components in the signal domain, than summation of the concentration time courses; this was reflected in smaller areas of residual high CBFratio in the same region. The ICA corrected rCBF and CBFratio images were lower intensity than that from mMVþVM. Figure 5 shows for pat1 and Figure 6 for pat2 the full parameters from the different variants of the modified vascular model, including the case where a (fixed) mild amount of dispersion of the AIF was included. The identification of the macro vascular component was similar with and without the inclusion of dispersion. Visual inspection indicated that a greater area of nonzero raBV was identified when dispersion was included particular for the model in which concentrations of the two components were summed. Table 3 shows changes in the vascular model parameters in regions identified as being subject to substantial and minimal MV contamination from the CBFratio of the VM results. (The vasular models paraemters can be seen in Figures 5 and 6.) Changes in model parameters were
286
Chappell et al.
FIG. 5. Vascular model parameters from pat1 including (as the second row for each model variant) the extra effect of (fixed) mild dispersion of the AIF.
large within the prospectively MV contaminated region, but small within the region in which minimal MV contamination might be expected. In the contaminated region, model-based MV correction led to an increase in both lambda and MTT, consistent with the estimation of signal arising from the capillary bed without a fast flowing MV component (i.e., more truly representative of microvascular perfusion). Table 3 includes the same analysis for the analysis including (fixed) dispersion of the AIF, the results of which were broadly consistent with those in the absence of dispersion: substantial differences in parameter values were seen in areas in which
macro vascular contamination was identified, but minimal elsewhere. DISCUSSION In this work a new method for the correction of MV signal contribution to DSC perfusion quantification has been proposed. This has extended the Vascular Modelbased approach of Mouridsen et al (11). The method has been shown to produce anatomically plausible images of (relative) “arterial” blood volume associated with the MV component, with high values in regions that correspond
Correcting Contamination in DSC Perfusion
287
FIG. 6. Vascular model parameters from pat2 including (as the second row for each model variant) the extra effect of (fixed) mild dispersion of the AIF.
to locations of major vessels in the brain, such in the major branches of the MCA. Accompanied with this are plausible reductions in the relative perfusion estimated within these regions along with changes in the parameters of the transit time distribution consistent with a better characterization of capillary flow. Further validation is still required to confirm whether the MV component is being fully and correctly identified using this approach. This validation may be, however, difficult due to the complexities of using a suitable gold standard. Two approaches to the combination of a tissue and MV contribution were considered, either by summing
the concentrations or MR signals. In theory, summing the concentrations may appear to be more accurate (given that concentrations do add up linearly in the voxel). However, DSC-MRI does not measure contrast concentration, but estimates its value from the signal intensity changes. So, we do not actually measure Ca(t), but we measure Sa(t) instead. Therefore, the approach based on summing the concentrations implicitly assumes that the measured concentration in the large artery has no partial volume effect, dispersion, or other sources or error (i.e., it measures the true AIF). Furthermore, it also assumes that the total concentration (calculated as the
288
Chappell et al.
Table 3 Percent Changes in Vascular Model Parameters after MV Correction in High and Low MV Contamination ROIs from Both Correction Choices of Signal or Concentration (conc) Time Curve Adding within the Model DCBF ratio % Contamination ROI High
Pat1 Pat2 Low Pat1 Pat2 Including mild dispersion High Pat1 Pat2 Low Pat1 Pat2
conc
signal
Dl % conc
DMTT % signal
conc
signal
18.6 60.6 0.55 0.03
65.3 74.9 0.48 0.2
263.0 20.5 1.2 18.3
806.0 305.2 2.9 27.1
33.0 60.8 0.17 1.7
36.2 42.7 0.39 3.1
16.3 65.1 19.5 17.1
61.0 68.6 19.3 15.8
281.5 7.8 9.5 8.2
547.1 212.8 9.0 12.2
31.4 62.2 3.8 5.6
38.8 48.5 3.6 7.0
sum of macro vascular and tissue concentrations) can be converted to a total signal intensity time course by a simple exponential relationship (modeling the R2* effect). This neglects, however, the phase terms from each of these contributions. Similarly, the approach based on summing the MR signals also neglects the phase contributions (i.e., it only sums the magnitude signals) and assumes that the true AIF has been measured. Both approaches have therefore the same underlying assumptions, and it is not clear which one should give the best result. It is for this reason that both approaches were tested to assess empirically their relative performance. The summation of MR signals was found to produce a more extensive raBV image, suggesting it may produce an improved correction the MV effect. The reason behind this empirical finding is not clear at present. In reality, the generation of MR signal in a DSC experiment arises also from R2* changes in the extra vascular space from Gadolinium present within the vasculature. Thus neither assumption of linear combination of concentrations or signal contributions represents a complete model of how the DSC-MRI signal would arise in a voxel containing a mix of tissue (which is a distributed capillary bed) and arterial vessels. In reality, the signals from arteries show some angular dependence (19), indicative of a limitation of the simple model considered here. More extensive modeling along the lines of the work in (7,20,21) might shed more light on this issue, as well as determining the level of accuracy the proposed correction method can be expected to achieve. This is, however, beyond the scope of this study. To produce a combined tissue plus MV model a modified version of the original VM has been proposed. This has added a further parameter to the model to be estimated from the data, MTT, instead of this information being derived by means of the CBV from the data. This permits a contribution to the blood volume in the voxel to be attributable to an MV component. Using a simulation study that attempted to match that used for the original VM in (11) the mVM (without an MV component) was found to be as accurate, if not slightly more so. As expected, however, the variability observed with the mVM was larger, consistent with the increased degrees of freedom. Discrepancies between the mVM and those
observed from the VM might also arise due to the different algorithms used. The correction method proposed here relies upon the accurate identification of the AIF, which is widely recognized challenge (22), particularly due to partial volume effects. Recently methods have been proposed that exploit both magnitude and/or phase images (7,20,21) to provide an improved estimate of the AIF. The proposed MV correction method will directly benefit from any such method, because a more accurate AIF should lead to a more accurate MV correction within the proposed Bayesian framework. Aside from a role in AIF detection, the enhancement in detection of arterial blood signal provided by phase data might also potentially be exploited in a more accurate model-based MV correction method. For example, magnitude data combined with a vessel enhancement filter (23) might be exploited to infer the orientation of the MV component within the voxel (24), which could then be incorporated into a more rigorous modeling of the way the signal contributions are combined, accounting for the fact that that summation occurs in the complex domain (7,20,21,25,26). In this study, we assumed a linear relationship between the change in relaxation rate and the contrast agent concentration for both the tissue and the AIF, despite the known quadratic dependency inside large vessels (15). It should be noted, however, that the AIF is usually measured from a branch of the middle cerebral artery (in our case, using an automatic AIF algorithm) (27). Given the size of these arterial branches (relative to the voxel size), the presence of partial volume is unavoidable, and it has been suggested that the best AIF estimates are actually measured from outside the vessel (20–22,26). Therefore, the relevant contrast agent relationship for the AIF is, in practice, the linear relationship. The situation would be different if the AIF is measured, for example, in the internal carotid artery (15), in which case the quadratic relationship should be used. The approach proposed here relies on the AIF as a model for the MV component. In reality the effects of dispersion will mean that the MV component distant from the location where the AIF was measured may differ in shape (16), leading to some error in MV
Correcting Contamination in DSC Perfusion
component quantification. In this study we investigated how the introduction of dispersion into the AIF used for analysis affected the results. Overall macro vascular contamination was identified in the same locations, with a slightly broader extent in the case of dispersion, particularly for the concentration-summing variant of the model. This might be consistent with the presence of a dispersed AIF within the vasculature and thus better identification of the macro vascular component in the data. However, the application of a dispersion kernel to the AIF may also make it more similar to the tissue signal, particularly in tissue with short MTT. This ambiguity means that it is unlikely to be possible to build a full model that includes dispersion, for example by a parameterized vascular transport function whose parameters are estimated as part of the analysis, a difficulty that has also been observed in combination with nonparametric deconvolution (18). An alternative approach to that taken here would be to attempt to explicitly model the AIF locally using a parameterized model to represent both the macro vascular contribution in any given voxel and the input function for the tissue. This is essentially the approach taken in (28) for a comparable ASL scenario. However, in ASL it is more common to define, as in that case, a dispersed AIF based on the idealized AIF created by the labeling. In DSC perfusion the direct measurement of the AIF is preferred, because there is likely to be greater variability associated with individual physiology and the injection profile for the contrast agent. In Lee et al (29), a parameterized AIF was used for DSC perfusion data to estimate a local AIF and that form could be used within the framework proposed here and the parameters estimated voxelwise. Whereas that work attempted to fit for the AIF in each voxel, unlike this work, no attempt was made to identify a separate macro vascular component simultaneously. A comparison has been made to the alternative ICA based approach for MV correction proposed by (3). An ICA approach, being “model-free,” is appealing because it can potentially remove MV components that deviate from the AIF. However, accurate and consistent identification of MV components within ICA components remains challenging, particularly in patients with severe vascular abnormalities (5), and we found it to not be possible on one of the datasets considered here. We do not, therefore, believe the ICA based approach provides a general solution to the problem, although it will give good results in some cases. The use of gradient-echo (GE) and spin-echo (SE) DSCMRI data (ideally acquired simultaneously (30) could be of great interest for the method proposed in this study. Due to their different sensitivity to vessel size (31) the macro vascular contribution is expected to be much larger for GE data, with SE data only sensitive, in theory, to microvascular effects. These data could, therefore, play a role, for example, to validate the proposed correction method. However, it should be noted that most DSC-MRI studies are based on an EPI-type acquisition, which has a nonnegligible T2* contribution, even for the case of SE-EPI (32). Therefore, the proposed method could in fact be used to assess the residual macro vascu-
289
lar contribution in SE DSC-MRI data. This is somewhat analogous to our previous work on macro vascular correction in arterial spin labeling (8,28), where a similar modeling was applied to flow crushed arterial spin labeling data to assess the residual macro vascular contribution due to incomplete dephasing of macro vascular signal by the crusher gradients. In summary, an extended vascular model was proposed to correct for MV contamination in DSC perfusion imaging within a Bayesian framework. Using a MV component based on the AIF, it appears to be possible to substantially remove MV contamination from DSC CBF estimates, thus more truly representative of microvascular perfusion. An implementation of the methods used here can be downloaded from the FMRIB Software Library (www.fmrib.ox.ac.uk/fsl/verbena). ACKNOWLEDGMENTS We thank Prof. Geoffrey A. Donnan and Dr. Amy Brodtmann (The Florey Institute of Neuroscience and Mental Health, Melbourne, Australia) for providing the MRI data. F.C. is grateful to the National Health and Medical Research Council (NHMRC) of Australia, the Australian Research Council (ARC), and the Victorian Government’s Operational Infrastructure Support Grant for their support. M.C. was employed by The Centre of Excellence in Personalized Healthcare funded by the Wellcome Trust and EPSRC. REFERENCES 1. Calamante F, Thomas DL, Pell GS, Wiersma J, Turner R. Measuring cerebral blood flow using magnetic resonance imaging techniques. J Cereb Blood Flow Metab 1999;19:701–735. 2. Willats L, Calamante F. The 39 steps: evading error and deciphering the secrets for accurate dynamic susceptibility contrast MRI. NMR Biomed 2012;26:913–931. 3. Carroll T, Haughton V, Rowley H, Cordes D. Confounding effect of large vessels on MR perfusion images analyzed with independent component analysis. Am J Neuroradiol 2002;23:1007–1012. 4. Reishofer G, Koschutnig K, Enzinger C, Ischebeck A, Keeling S, Stollberger R, Ebner F. Automated macrovessel artifact correction in dynamic susceptibility contrast magnetic resonance imaging using independent component analysis. Magn Reson Med 2011;65:848–857. 5. Calamante F, Mïrup M, Hansen LK. Defining a local arterial input function for perfusion MRI using independent component analysis. Magn Reson Med 2004;52:789–797. 6. Calamante F, Gadian DG, Connelly A. Delay and dispersion effects in dynamic susceptibility contrast MRI: simulations using singular value decomposition. Magn Reson Med 2000;44:466–473. 7. Bleeker EJW, Webb AG, van Walderveen MAA, van Buchem MA, van Osch MJP. Evaluation of signal formation in local arterial input function measurements of dynamic susceptibility contrast MRI. Magn Reson Med 2012;67:1324–1331. 8. Chappell MA, MacIntosh BJ, Donahue MJ, G€ unther M, Jezzard P, Woolrich MW. Separation of macrovascular signal in multi-inversion time arterial spin labelling MRI. Magn Reson Med 2010;63:1357– 1365. 9. Chappell MA, Mehndiratta A, Payne SJ, Calamante F. Bayesian model-based correction for macro vascular signal in dynamic susceptibility contrast perfusion MRI. In Proceedings of the 20th Annual Meeting of ISMRM, Melbourne, Australia, 2012. Abstract 1950. 10. Mehndiratta A, MacIntosh BJ, Crane DE, Payne SJ, Chappell MA. Improved vascular model-based analysis for DSC-MRI perfusion quantification. In Proceedings of the 20th Annual Meeting of ISMRM, Melbourne, Australia, 2012. Abstract 1948.
290 11. Mouridsen K, Friston K, Hjort N, Gyldensted L, Østergaard L, Kiebel S. Bayesian estimation of cerebral perfusion using a physiological model of microvasculature. Neuroimage 2006;33:570–579. 12. Wu O, Østergaard 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–174. 13. Chappell M, Groves A, Whitcher B, Woolrich M. Variational Bayesian inference for a nonlinear forward model. IEEE Trans Signal Process 2009;57:223–236. 14. Mackay D. Probable Networks and plausible predictions - a review of practical Bayesian methods for supervised neural networks. Netw Comput Neural Syst 1995;6:469–505. 15. Beckmann CF, Smith SM. Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans Med Imaging, 2004;23:137–152. 16. Calamante F. Bolus dispersion issues related to the quantification of perfusion MRI data. J Magn Reson Imaging 2005;22:718–722. 17. Calamante F, Yim PJ, Cebral JR. Estimation of bolus dispersion effects in perfusion MRI using image-based computational fluid dynamics. Neuroimage 2003;19:341–353. 18. Mehndiratta A, Calamante F, MacIntosh BJ, Crane DE, Payne SJ, Chappell MA. Modeling and correction of bolus dispersion effects in dynamic susceptibility contrast MRI. Magn Reson Med 2014;72:1762– 1774. 19. Van Osch M, Vonken E, Viergever M, Van Der Grond J, Bakker C. Measuring the arterial input function with gradient echo sequences. Magn Reson Med 2003;49:1067–1076. 20. Bleeker EJW, van Osch MJP, Connelly A, van Buchem MA, Webb AG, Calamante F. New criterion to aid manual and automatic selection of the arterial input function in dynamic susceptibility contrast MRI. Magn Reson Med 2011;65:448–456. 21. Bleeker E, Van Buchem M, Van Osch M. Optimal location for arterial input function measurements near the middle cerebral artery in firstpass perfusion MRI. J Cereb Blood Flow Metab 2009;29:840–852.
Chappell et al. 22. Calamante F. Arterial input function in perfusion MRI: a comprehensive review. Prog Nucl Magn Reson Spectrosc 2013;74:1–32. 23. Frangi AF, Niessen WJ, Vincken KL, Viergever MA. Multiscale vessel enhancement filtering. Med Image Comput Comput Assist Interv 1998;1496:130–137. 24. Straka M, Newbould RD, Sramek M, Albers GW, Bammer R. Improving DSC-MRI by orientation-corrected phase-based AIF and VOF. In Proceedings of the 18th Annual Meeting of ISMRM, Stockholm, Sweden, 2010. Abstract 622. 25. Kjolby B, Mikkelsen I, Pedersen M, Ostergaard L, Kiselev V. Analysis of partial volume effects on arterial input functions using gradient echo: a simulation study. Magn Reson Med 2009;61:1300–1309. 26. 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:514–523. €stergaard L. Automatic 27. Mouridsen K, Christensen S, Gyldensted L, y selection of arterial input function using cluster analysis. Magn Reson Med 2006;55:524–531. 28. Chappell MA, Woolrich MW, Petersen ET, Golay X, Payne SJ. Comparing model-based and model-free analysis methods for QUASAR arterial spin labeling perfusion quantification. Magn Reson Med 2013;69:1466–1475. 29. Lee JJ, Bretthorst GL, Derdeyn CP, Powers WJ, Videen TO, Snyder AZ, Markham J, Shimony JS. Dynamic susceptibility contrast MRI with localized arterial input functions. Magn Reson Med 2010;63: 1305–1314. 30. Schmiedeskamp H, Straka M, Newbould RD, Zaharchuk G, Andre JB, Olivot J-M, Moseley ME, Albers GW, Bammer R. Combined spin- and gradient-echo perfusion-weighted imaging. Magn Reson Med 2012;68: 30–40. 31. Weisskoff R, Zuo C, Boxerman J, Rosen B. Microscopic susceptibility variation and transverse relaxation: theory and experiment. Magn Reson Med 1994;31:601–610. 32. Birn RM, Bandettini PA. The effect of T2’ changes on spin-echo EPIderived brain activation maps. In Proceedings of the 10th Scientific Meeting of ISMRM, Honolulu, Hawaii, USA, 2002. p. 1324.
Magnetic Resonance in Medicine 74:280–290 (2015)
Correcting for Large Vessel Contamination in Dynamic Susceptibility Contrast Perfusion MRI by Extension to a Physiological Model of the Vasculature Michael A. Chappell,1* Amit Mehndiratta,1 and Fernando Calamante2,3 Purpose: Dynamic susceptibility contrast (DSC) perfusion images are contaminated by contributions from macro vascular signal arising from contrast agent within the larger arteries that do not contribute directly to the local tissue perfusion. Methods: A vascular model of the DSC perfusion signal was extended by the inclusion of a macro vascular component based on the arterial input function. This was implemented within a Bayesian nonlinear model-fitting algorithm that included automatic model complexity reduction. Results were compared with existing methods that do not correct for the macro vascular contamination as well as an independent component analysis technique. Results: Macro vascular signal was identified in regions corresponding to larger arteries resulting in reductions by 62% within a region of interest identified with high contamination. Whereas visually similar results could be achieved with independent component analysis, it resulted in reductions in global tissue perfusion and was not robustly applicable to patient data. Conclusion: A model-based strategy for correction of macro vascular contamination in DSC perfusion images is feasible, although the model may currently need extending to more accurately account for nonlinear effects of contrast agent in C 2014 large arteries. Magn Reson Med 74:280–290, 2015. V Wiley Periodicals, Inc. Key words: cerebral blood flow (CBF); dynamic susceptibility contrast MRI; vascular model; Bayesian analysis; macro vascular contamination
INTRODUCTION The quantification of dynamic susceptibility contrast (DSC) perfusion-weighted MRI typically relies upon the identification of an arterial input function (AIF) and subsequent voxelwise numerical deconvolution to calculate 1 Institute of Biomedical Engineering, University of Oxford, ORCRB, Old Road Campus, Headington, Oxford, United Kingdom. 2 Florey Institute of Neuroscience and Mental Health, Heidelberg, Australia. 3 Department of Medicine, Austin Health and Northern Health, University of Melbourne, Melbourne, Australia. Grant sponsor: the Wellcome Trust and EPSRC; Grant number: WT088877/ Z/09/Z. *Correspondence to: Michael A. Chappell, D.Phil., Institute of Biomedical Engineering, University of Oxford, ORCRB, Old Road Campus, Headington, Oxford, OX3 7DQ, UK. E-mail: [email protected] Received 10 June 2014; revised 4 July 2014; accepted 7 July 2014 DOI 10.1002/mrm.25390 Published online 8 August 2014 in Wiley Online Library (wileyonlinelibrary. com). C 2014 Wiley Periodicals, Inc. V
perfusion, mean transit time (MTT) or other quantitative parameters (1) A fundamental assumption in DSC-MRI is that the signal in the voxel arises from the Gadolinium contrast agent within the micro vasculature (primarily capillaries, but also arterioles and venules). Thus the residue function estimated by means of deconvolution represents the distribution of transit times within the voxel and the MTT is the mean of this distribution. However, a proportion of voxels will also include a contribution from contrast agent within larger vessel, for example those proximal to the middle cerebral arteries, either because the vessel is within the voxel or due to susceptibility effects of a vessel nearby. Therefore, this extra contribution results in errors to CBF, MTT, and cerebral blood volume (CBV) quantification, because the measured signal does not reflect solely micro vasculature (2). Despite this, the presence of this macro vascular (MV) contribution is typically ignored. Various methods have been proposed to deal with this “unwanted” MV contribution (NB. in fact, its presence is critical to the ability to extract AIF from the data in the first place). For example, Carroll et al (3) proposed that independent component analysis (ICA) could be used as a preprocessing step to identify and remove components from the signal that could be regarded as macro vascular. Subsequently, Reishofer et al (4) extended this approach, by proposing an automated scheme to identify the relevant MV components. Such an approach is however limited by how accurately appropriate MV components can be identified in practice, particularly in patients with severe vascular abnormalities (5). In principle, the contribution of the MV component to the signal should match that of the AIF, ignoring dispersive effects (6) and complex partial volume effects (7); thus, it should be possible to use the existing information about the measured AIF to assist in MV correction. In this work, we propose a model-based solution to the correction for MV contamination inspired by related recent work in Arterial Spin Labeling (8), building upon work previously presented in abstract form (9,10). To do so, we build upon the vascular model of Mouridsen et al (11) for the DSC-MRI tissue signal, by incorporating a separate MV component, which must be estimated simultaneously from the DSC-MRI data. Like the similar methodology used in Arterial Spin Labeling (8), we subject the extra MV component to an automated model complexity reduction procedure to ensure that the extra
280
Correcting Contamination in DSC Perfusion
281
complexity introduced by adding this component is only used where the data justifies it. Thus, avoiding over fitting of the data in regions where no large vessels are present. The method was tested on in vivo data, and compared with the ICA-based approach. THEORY
mate and apply it to the tissue, because there will be some contribution from the macro vascular blood signal. This can be addressed by a modified version of the vascular model, termed hereafter the mVM. In this model, MTT becomes a further parameter to be estimated from the data within the model fitting, with b ¼ MTT/l. A total of four parameters have be estimated in the mVM: rCBF, k, and d and MTT.
Vascular Model The concentration of contrast agent within a voxel, C(t), can be described according to tracer kinetics as the convolution of an AIF, Ca(t), and a tissue residue function, R(t), scaled by the CBF (1): CðtÞ ¼ CBF ðCa ðtÞ RðtÞÞ¼CBF
Z
t
Ca ðtÞ RðttÞdt
[1]
0
In the vascular model of Mouridsen et al (11), the residue function was parameterized using a model of contrast agent transit through the capillary network that followed a gamma distribution: RðtÞ ¼
Z
1
hðt; a; bÞdt t
hðt; a; bÞ ¼
t 1 t a1 exp b ; a; b > 0 ba GðaÞ
[2]
where a and b are the shape and scale parameters of the distribution respectively, which provides a flexible model to describe a large range of feasible residue function shapes. Because the product ab is the mean of the distribution, this equates to the MTT for contrast agent transit. Mouridsen et al (11) used the central volume principle to reparameterize the distribution in terms of relative CBF (rCBF), relative CBV (rCBV) and a shape parameter: a ¼ l, b ¼ rCBV/(k .CBV). A further parameter d was included to account for differences in arrival between the AIF and the voxel concentration time course, i.e., bolus delay. The estimated concentration time course was converted to a predicted signal assuming that the transverse relaxation rate changes linearly with concentration of the contrast agent SðtÞ ¼ Sð0Þer2 CðtÞTE , where S(0) is the baseline signal before bolus arrival and r2 was incorporated into the CBF parameter so that relative CBF was estimated. Thus there were three parameters in the model: rCBF, k and d. These were estimated within a Bayesian nonlinear model-fitting algorithm, where the relative CBV, to be used in the expression for the shape parameter b, had previously been estimated from the data in the normal way by taking the ratio of the area under the concentration time curve and arterial input function. Priors were applied to the parameters using the estimates from a numerical deconvolution of the data, using the singular value decomposition (SVD) approach of (12), to set the mean for rCBF. A potential disadvantage of the original VM implementation was that the model fitting was constrained by the rCBV value obtained from the data, which due to noise could introduce a bias into the model fitting. A further issue arises in the context of modeling both a tissue and macro vascular component within a single voxel, as is the aim of this work. Under these conditions, it is no longer possible to take the standard rCBV esti-
Macro Vascular Component The presence of a macro vascular contribution within the voxel can be modeled by the addition of an extra component with a concentration time course that matches that of the AIF. This assumes that the flow of the contrast through the voxel is sufficiently rapid that there is no accumulation of the agent within the voxel, unlike the capillary signal. This is expected to be true of larger arteries and arterioles (8). This also neglects any dispersion effects as well as interactions between tissue and MV signals (see discussion). The macro vascular model can parameterized by a magnitude that represents the relative arterial blood volume (raBV) and an arrival time, da, that accounts for delay between the AIF measured elsewhere in the brain and its observation as a macro vascular signal in the voxel. The blood volume parameter relates to the magnitude of the AIF, because this is used as the basis for this model component, given that the scaling of this to absolute blood volume is generally not known only a relative blood volume can be estimated. The arrival time need not match that of the tissue if, for example, the blood in the artery were destined for elsewhere in the brain. This extension to the modified version of the vascular model is termed hereafter the mVMþMV, and can be used to correct for MV contamination. Two variants for the inclusion of the MV component in mVMþMV were considered, these are illustrated in Figure 1: 1. Combining the mVM and MV contributions in the concentration time course: CðtÞ ¼ rCBF ðCa ðtÞ RðtÞÞ þ raBV Ca ðtÞ SðtÞ ¼ Sð0Þer2 CðtÞTE
[3]
2. Combining the corresponding individual signal time courses: CðtÞ ¼ rCBF ðCa ðtÞ RðtÞÞ SðtÞ ¼ Sð0Þðer2 CðtÞTE þ er2 raBV CaðtÞTE Þ
[4]
Note that, whereas the rCBF parameter values might be comparable between the two models, the same may not be true of the raBV parameters. METHODS Model Implementation The mVM and mVMþMV were implemented within a Bayesian nonlinear model fitting algorithm (13), this method differs from that originally used by Mouridsen et al (11), but similarly permits prior information to be
282
Chappell et al.
FIG. 1. Illustration of the mVMþVM model showing both the summation of the concentration time curves (a) and the summation of the predicted signals (b) arising from the separate tissue and MV concentration time curves.
incorporated within the analysis. An implementation of the methods used here can be found in the FMRIB Software Library (www.fmrib.ox.ac.uk/fsl/verbena). The priors used for the three models are given in Table 1. Priors for the mVM were chosen to either be noninformative (for the rCBF parameter) or based on physiologically expected values; unlike the VM, no prior information was derived from the data themselves. For the MV component an Automatic Relevancy Determination (or shrinkage) prior was placed on the raBV parameter (14), the role of this prior was to suppress the MV component in voxels where the data did not support its inclusion. This reflects the expectation that MV signal will not be present in all voxels in the brain and thus this component should be removed where it is not supported to avoid over fitting; for more information see Chappell et al (8,13). For comparison, data were also analyzed using a more commonly used SVD-based approach (using the delay insensitive oSVD method from Wu et al) (12). The results
Table 1 Priors Specified for the Model Fitting Given as the Mean and Standard Deviation of a Normal Distribution Parameter rCBF d (sec) log(k) log(MTT) raBV da (sec) S(0)
mVM
mVMþMV
0, 106 0, 5 2.3,1 1.5,0.3 – – 0,106
0, 106 0, 5 2.3,1 2.3,0.3 0, 106 0,5 0,106
from the SVD analysis were also used to initialize the parameters in the VM and mVM model fitting, namely estimated rCBF and MTT were taken from the SVD analysis and used as the initial values for these parameters. Additionally, the S(0) parameter was initialized using the first volume in the data series. For the VM the SVD results were also used in the specification of the priors as in the original implementation. ICA was also used in an attempt to remove MV contamination. ICA was performed using MELODIC (15) from the FMRIB software library (www.fmrib.ox.ac.ul/fsl), and the MV components were identified manually and removed from the data before oSVD analysis. An important confound for quantification of perfusion using DSC is the dispersion of the arterial bolus. To investigate how this might impact the macro vascular contamination correction proposed here we repeated the model-based analysis above including dispersion of the AIF. The effects of dispersion were modeled by the numerical convolution of the data-derived AIF with a vascular transport function (VTF) (16,17). We used a gamma-variate VTF and chose the parameters to simulate a “low” amount of dispersion following (18) (namely shape parameters: p ¼ 1 and s ¼ 2 in Equation 3 of Mehndiratta et al) (18). The dispersion was taken as fixed and no attempt was made to estimate it directly from the data (see the Discussion section). Vascular Model Comparisons: Numerical Simulations First, the mVM was evaluated in comparison to the original VM and oSVD in the absence of a macro vascular component to identify if any errors might be introduced into the analysis by the use of the modified version.
Correcting Contamination in DSC Perfusion
283
FIG. 2. Estimated CBF values from the simulations for selected cases of residue function shape, noise and delay [c.f. (11) Figure 5], showing mean estimated CBF versus true value with error bars showing standard deviation across the separate realizations of the noise.
Comparisons were performed using the same simulation method as in Mouridsen et al (11). In brief: an AIF was simulated using an analytical gamma-variate function. Simulated data were generated for validation of the method with a CBV of 4 mL/100 g (typical gray matter), CBF in the range 10–70 mL/100 g/min in increment of 10 mL/100 g/min, delay of 0 and 65 s, and two residue
functions: exponential & box (approximated using k ¼ 100). Simulated signals were generated from concentration curves with baseline signal S(0) ¼ 100, echo-time TE ¼ 65 ms, the proportionality constant k was chosen such that 40% peak signal drop was attained at a CBF of 60 mL/100 g/min with CBV ¼ 4%, and 60% peak signal drop for the AIF signal. Zero mean Gaussian noise was
284
Chappell et al.
Table 2 Mean Ratios between Estimated and True CBF Values from the Simulations (Mean 6 Standard Deviation) SNR¼100
Exponential Exponential, delay¼5 Box-car
SNR¼20
mVM
SVD
mVM
SVD
1.0260.15 0.9460.24 1.0860.13
0.8660.11 0.8560.11 1.1560.09
0.9360.43 0.7960.37 1.0960.31
0.7560.22 0.8660.25 1.0560.30
added to the signals and AIF to achieve a given SNR, with both SNR of 100 and 20 being simulated. For each combination of CBF, residue function, and delay a total of 100 noisy signals and AIF were generated. The methods were compared by calculating the regression coefficient between the estimated rCBF and the true value over the range considered. In Vivo Data The MV correction method was evaluated on data from two patients. Pat1: left carotid stenosis (scanned postcarotid endarterectomy); Pat2: stroke patient with right middle cerebral artery (MCA) occlusion (scanned 9 hrs postsymptoms onset). DSC-MRI data were acquired at 3 Tesla (T) (Pat1) or at 1.5T (Pat2), using a gradient-echo EPI, TE/TR ¼ 20/1250 ms (Pat1) or 60/1740 ms (Pat2), and a single (Pat1) or double (Pat2) dose of contrast
agent. The AIF was measured on a branch of the contralateral MCA using the semiautomatic method in Penguin (www.cfin.au.dk/software/penguin). Macro Vascular Correction Comparisons To examine the presence of MV contamination in the estimated rCBF images, the ratio of rCBF (CBFratio) to that in a region of interest (ROI) chosen within the deep white matter (without MV effects) was calculated. To quantify the changes in rCBF and parameters of the VM, two further ROIs were defined that prospectively represented areas of substantial and negligible MV contamination, respectively. These ROI were defined as voxels whose CBFratio from the VM analysis was more than 4.5 (high contamination) and less than 1.5 (low contamination). Within these ROIs, the mean change in CBFratio, l, and MTT were calculated.
FIG. 3. Results for CBF quantification from pat1 showing rCBF and CBFratio from all analysis methods considered, and raBV in the mVMþMV methods.
Correcting Contamination in DSC Perfusion
285
FIG. 4. Results for CBF quantification from pat2 showing rCBF and CBFratio from all analysis methods except ICA correction, and raBV in the mVMþMV methods.
RESULTS Vascular Model Comparisons: Numerical Simulations Figure 2 shows estimated rCBF values plotted against the true value across the simulated data. Table 2 gives the mean ratios between estimated and true CBF along with the standard deviation across the different noise realizations. The mVM was generally found to produce more accurate CBF estimates than SVD with wider standard deviation. Comparing to the VM (under a similar simulation protocol) in Mouridsen et al (11) and their Figure 5 and Table 1, very similar performance was observed for the case without bolus delay, including overestimation of CBF for a box-car residue function shape. The mVM appeared to have higher accuracy for data subject to delay. Macro Vascular Correction Comparisons Figure 3 shows the rCBF and CBFratio images in pat1 for all the methods considered, along with raBV from the mVMþMV approaches. Figure 4 shows the same information for pat2, except for the results from ICA correction because it was not possible to identify ICA components that were suitably representative of MV signal in this case. There was a substantial reduction in rCBF and CBFratio observed in regions in which large
arteries would be expected in all of the MV correction approaches. Some residual regions of high CBFratio persisted in the mVMþMV methods. The raBV images had high magnitude in regions consistent with the presence of MCA branches. Regions of high intensity of raBV were more extensive in the mVMþMV method using the summation of components in the signal domain, than summation of the concentration time courses; this was reflected in smaller areas of residual high CBFratio in the same region. The ICA corrected rCBF and CBFratio images were lower intensity than that from mMVþVM. Figure 5 shows for pat1 and Figure 6 for pat2 the full parameters from the different variants of the modified vascular model, including the case where a (fixed) mild amount of dispersion of the AIF was included. The identification of the macro vascular component was similar with and without the inclusion of dispersion. Visual inspection indicated that a greater area of nonzero raBV was identified when dispersion was included particular for the model in which concentrations of the two components were summed. Table 3 shows changes in the vascular model parameters in regions identified as being subject to substantial and minimal MV contamination from the CBFratio of the VM results. (The vasular models paraemters can be seen in Figures 5 and 6.) Changes in model parameters were
286
Chappell et al.
FIG. 5. Vascular model parameters from pat1 including (as the second row for each model variant) the extra effect of (fixed) mild dispersion of the AIF.
large within the prospectively MV contaminated region, but small within the region in which minimal MV contamination might be expected. In the contaminated region, model-based MV correction led to an increase in both lambda and MTT, consistent with the estimation of signal arising from the capillary bed without a fast flowing MV component (i.e., more truly representative of microvascular perfusion). Table 3 includes the same analysis for the analysis including (fixed) dispersion of the AIF, the results of which were broadly consistent with those in the absence of dispersion: substantial differences in parameter values were seen in areas in which
macro vascular contamination was identified, but minimal elsewhere. DISCUSSION In this work a new method for the correction of MV signal contribution to DSC perfusion quantification has been proposed. This has extended the Vascular Modelbased approach of Mouridsen et al (11). The method has been shown to produce anatomically plausible images of (relative) “arterial” blood volume associated with the MV component, with high values in regions that correspond
Correcting Contamination in DSC Perfusion
287
FIG. 6. Vascular model parameters from pat2 including (as the second row for each model variant) the extra effect of (fixed) mild dispersion of the AIF.
to locations of major vessels in the brain, such in the major branches of the MCA. Accompanied with this are plausible reductions in the relative perfusion estimated within these regions along with changes in the parameters of the transit time distribution consistent with a better characterization of capillary flow. Further validation is still required to confirm whether the MV component is being fully and correctly identified using this approach. This validation may be, however, difficult due to the complexities of using a suitable gold standard. Two approaches to the combination of a tissue and MV contribution were considered, either by summing
the concentrations or MR signals. In theory, summing the concentrations may appear to be more accurate (given that concentrations do add up linearly in the voxel). However, DSC-MRI does not measure contrast concentration, but estimates its value from the signal intensity changes. So, we do not actually measure Ca(t), but we measure Sa(t) instead. Therefore, the approach based on summing the concentrations implicitly assumes that the measured concentration in the large artery has no partial volume effect, dispersion, or other sources or error (i.e., it measures the true AIF). Furthermore, it also assumes that the total concentration (calculated as the
288
Chappell et al.
Table 3 Percent Changes in Vascular Model Parameters after MV Correction in High and Low MV Contamination ROIs from Both Correction Choices of Signal or Concentration (conc) Time Curve Adding within the Model DCBF ratio % Contamination ROI High
Pat1 Pat2 Low Pat1 Pat2 Including mild dispersion High Pat1 Pat2 Low Pat1 Pat2
conc
signal
Dl % conc
DMTT % signal
conc
signal
18.6 60.6 0.55 0.03
65.3 74.9 0.48 0.2
263.0 20.5 1.2 18.3
806.0 305.2 2.9 27.1
33.0 60.8 0.17 1.7
36.2 42.7 0.39 3.1
16.3 65.1 19.5 17.1
61.0 68.6 19.3 15.8
281.5 7.8 9.5 8.2
547.1 212.8 9.0 12.2
31.4 62.2 3.8 5.6
38.8 48.5 3.6 7.0
sum of macro vascular and tissue concentrations) can be converted to a total signal intensity time course by a simple exponential relationship (modeling the R2* effect). This neglects, however, the phase terms from each of these contributions. Similarly, the approach based on summing the MR signals also neglects the phase contributions (i.e., it only sums the magnitude signals) and assumes that the true AIF has been measured. Both approaches have therefore the same underlying assumptions, and it is not clear which one should give the best result. It is for this reason that both approaches were tested to assess empirically their relative performance. The summation of MR signals was found to produce a more extensive raBV image, suggesting it may produce an improved correction the MV effect. The reason behind this empirical finding is not clear at present. In reality, the generation of MR signal in a DSC experiment arises also from R2* changes in the extra vascular space from Gadolinium present within the vasculature. Thus neither assumption of linear combination of concentrations or signal contributions represents a complete model of how the DSC-MRI signal would arise in a voxel containing a mix of tissue (which is a distributed capillary bed) and arterial vessels. In reality, the signals from arteries show some angular dependence (19), indicative of a limitation of the simple model considered here. More extensive modeling along the lines of the work in (7,20,21) might shed more light on this issue, as well as determining the level of accuracy the proposed correction method can be expected to achieve. This is, however, beyond the scope of this study. To produce a combined tissue plus MV model a modified version of the original VM has been proposed. This has added a further parameter to the model to be estimated from the data, MTT, instead of this information being derived by means of the CBV from the data. This permits a contribution to the blood volume in the voxel to be attributable to an MV component. Using a simulation study that attempted to match that used for the original VM in (11) the mVM (without an MV component) was found to be as accurate, if not slightly more so. As expected, however, the variability observed with the mVM was larger, consistent with the increased degrees of freedom. Discrepancies between the mVM and those
observed from the VM might also arise due to the different algorithms used. The correction method proposed here relies upon the accurate identification of the AIF, which is widely recognized challenge (22), particularly due to partial volume effects. Recently methods have been proposed that exploit both magnitude and/or phase images (7,20,21) to provide an improved estimate of the AIF. The proposed MV correction method will directly benefit from any such method, because a more accurate AIF should lead to a more accurate MV correction within the proposed Bayesian framework. Aside from a role in AIF detection, the enhancement in detection of arterial blood signal provided by phase data might also potentially be exploited in a more accurate model-based MV correction method. For example, magnitude data combined with a vessel enhancement filter (23) might be exploited to infer the orientation of the MV component within the voxel (24), which could then be incorporated into a more rigorous modeling of the way the signal contributions are combined, accounting for the fact that that summation occurs in the complex domain (7,20,21,25,26). In this study, we assumed a linear relationship between the change in relaxation rate and the contrast agent concentration for both the tissue and the AIF, despite the known quadratic dependency inside large vessels (15). It should be noted, however, that the AIF is usually measured from a branch of the middle cerebral artery (in our case, using an automatic AIF algorithm) (27). Given the size of these arterial branches (relative to the voxel size), the presence of partial volume is unavoidable, and it has been suggested that the best AIF estimates are actually measured from outside the vessel (20–22,26). Therefore, the relevant contrast agent relationship for the AIF is, in practice, the linear relationship. The situation would be different if the AIF is measured, for example, in the internal carotid artery (15), in which case the quadratic relationship should be used. The approach proposed here relies on the AIF as a model for the MV component. In reality the effects of dispersion will mean that the MV component distant from the location where the AIF was measured may differ in shape (16), leading to some error in MV
Correcting Contamination in DSC Perfusion
component quantification. In this study we investigated how the introduction of dispersion into the AIF used for analysis affected the results. Overall macro vascular contamination was identified in the same locations, with a slightly broader extent in the case of dispersion, particularly for the concentration-summing variant of the model. This might be consistent with the presence of a dispersed AIF within the vasculature and thus better identification of the macro vascular component in the data. However, the application of a dispersion kernel to the AIF may also make it more similar to the tissue signal, particularly in tissue with short MTT. This ambiguity means that it is unlikely to be possible to build a full model that includes dispersion, for example by a parameterized vascular transport function whose parameters are estimated as part of the analysis, a difficulty that has also been observed in combination with nonparametric deconvolution (18). An alternative approach to that taken here would be to attempt to explicitly model the AIF locally using a parameterized model to represent both the macro vascular contribution in any given voxel and the input function for the tissue. This is essentially the approach taken in (28) for a comparable ASL scenario. However, in ASL it is more common to define, as in that case, a dispersed AIF based on the idealized AIF created by the labeling. In DSC perfusion the direct measurement of the AIF is preferred, because there is likely to be greater variability associated with individual physiology and the injection profile for the contrast agent. In Lee et al (29), a parameterized AIF was used for DSC perfusion data to estimate a local AIF and that form could be used within the framework proposed here and the parameters estimated voxelwise. Whereas that work attempted to fit for the AIF in each voxel, unlike this work, no attempt was made to identify a separate macro vascular component simultaneously. A comparison has been made to the alternative ICA based approach for MV correction proposed by (3). An ICA approach, being “model-free,” is appealing because it can potentially remove MV components that deviate from the AIF. However, accurate and consistent identification of MV components within ICA components remains challenging, particularly in patients with severe vascular abnormalities (5), and we found it to not be possible on one of the datasets considered here. We do not, therefore, believe the ICA based approach provides a general solution to the problem, although it will give good results in some cases. The use of gradient-echo (GE) and spin-echo (SE) DSCMRI data (ideally acquired simultaneously (30) could be of great interest for the method proposed in this study. Due to their different sensitivity to vessel size (31) the macro vascular contribution is expected to be much larger for GE data, with SE data only sensitive, in theory, to microvascular effects. These data could, therefore, play a role, for example, to validate the proposed correction method. However, it should be noted that most DSC-MRI studies are based on an EPI-type acquisition, which has a nonnegligible T2* contribution, even for the case of SE-EPI (32). Therefore, the proposed method could in fact be used to assess the residual macro vascu-
289
lar contribution in SE DSC-MRI data. This is somewhat analogous to our previous work on macro vascular correction in arterial spin labeling (8,28), where a similar modeling was applied to flow crushed arterial spin labeling data to assess the residual macro vascular contribution due to incomplete dephasing of macro vascular signal by the crusher gradients. In summary, an extended vascular model was proposed to correct for MV contamination in DSC perfusion imaging within a Bayesian framework. Using a MV component based on the AIF, it appears to be possible to substantially remove MV contamination from DSC CBF estimates, thus more truly representative of microvascular perfusion. An implementation of the methods used here can be downloaded from the FMRIB Software Library (www.fmrib.ox.ac.uk/fsl/verbena). ACKNOWLEDGMENTS We thank Prof. Geoffrey A. Donnan and Dr. Amy Brodtmann (The Florey Institute of Neuroscience and Mental Health, Melbourne, Australia) for providing the MRI data. F.C. is grateful to the National Health and Medical Research Council (NHMRC) of Australia, the Australian Research Council (ARC), and the Victorian Government’s Operational Infrastructure Support Grant for their support. M.C. was employed by The Centre of Excellence in Personalized Healthcare funded by the Wellcome Trust and EPSRC. REFERENCES 1. Calamante F, Thomas DL, Pell GS, Wiersma J, Turner R. Measuring cerebral blood flow using magnetic resonance imaging techniques. J Cereb Blood Flow Metab 1999;19:701–735. 2. Willats L, Calamante F. The 39 steps: evading error and deciphering the secrets for accurate dynamic susceptibility contrast MRI. NMR Biomed 2012;26:913–931. 3. Carroll T, Haughton V, Rowley H, Cordes D. Confounding effect of large vessels on MR perfusion images analyzed with independent component analysis. Am J Neuroradiol 2002;23:1007–1012. 4. Reishofer G, Koschutnig K, Enzinger C, Ischebeck A, Keeling S, Stollberger R, Ebner F. Automated macrovessel artifact correction in dynamic susceptibility contrast magnetic resonance imaging using independent component analysis. Magn Reson Med 2011;65:848–857. 5. Calamante F, Mïrup M, Hansen LK. Defining a local arterial input function for perfusion MRI using independent component analysis. Magn Reson Med 2004;52:789–797. 6. Calamante F, Gadian DG, Connelly A. Delay and dispersion effects in dynamic susceptibility contrast MRI: simulations using singular value decomposition. Magn Reson Med 2000;44:466–473. 7. Bleeker EJW, Webb AG, van Walderveen MAA, van Buchem MA, van Osch MJP. Evaluation of signal formation in local arterial input function measurements of dynamic susceptibility contrast MRI. Magn Reson Med 2012;67:1324–1331. 8. Chappell MA, MacIntosh BJ, Donahue MJ, G€ unther M, Jezzard P, Woolrich MW. Separation of macrovascular signal in multi-inversion time arterial spin labelling MRI. Magn Reson Med 2010;63:1357– 1365. 9. Chappell MA, Mehndiratta A, Payne SJ, Calamante F. Bayesian model-based correction for macro vascular signal in dynamic susceptibility contrast perfusion MRI. In Proceedings of the 20th Annual Meeting of ISMRM, Melbourne, Australia, 2012. Abstract 1950. 10. Mehndiratta A, MacIntosh BJ, Crane DE, Payne SJ, Chappell MA. Improved vascular model-based analysis for DSC-MRI perfusion quantification. In Proceedings of the 20th Annual Meeting of ISMRM, Melbourne, Australia, 2012. Abstract 1948.
290 11. Mouridsen K, Friston K, Hjort N, Gyldensted L, Østergaard L, Kiebel S. Bayesian estimation of cerebral perfusion using a physiological model of microvasculature. Neuroimage 2006;33:570–579. 12. Wu O, Østergaard 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–174. 13. Chappell M, Groves A, Whitcher B, Woolrich M. Variational Bayesian inference for a nonlinear forward model. IEEE Trans Signal Process 2009;57:223–236. 14. Mackay D. Probable Networks and plausible predictions - a review of practical Bayesian methods for supervised neural networks. Netw Comput Neural Syst 1995;6:469–505. 15. Beckmann CF, Smith SM. Probabilistic independent component analysis for functional magnetic resonance imaging. IEEE Trans Med Imaging, 2004;23:137–152. 16. Calamante F. Bolus dispersion issues related to the quantification of perfusion MRI data. J Magn Reson Imaging 2005;22:718–722. 17. Calamante F, Yim PJ, Cebral JR. Estimation of bolus dispersion effects in perfusion MRI using image-based computational fluid dynamics. Neuroimage 2003;19:341–353. 18. Mehndiratta A, Calamante F, MacIntosh BJ, Crane DE, Payne SJ, Chappell MA. Modeling and correction of bolus dispersion effects in dynamic susceptibility contrast MRI. Magn Reson Med 2014;72:1762– 1774. 19. Van Osch M, Vonken E, Viergever M, Van Der Grond J, Bakker C. Measuring the arterial input function with gradient echo sequences. Magn Reson Med 2003;49:1067–1076. 20. Bleeker EJW, van Osch MJP, Connelly A, van Buchem MA, Webb AG, Calamante F. New criterion to aid manual and automatic selection of the arterial input function in dynamic susceptibility contrast MRI. Magn Reson Med 2011;65:448–456. 21. Bleeker E, Van Buchem M, Van Osch M. Optimal location for arterial input function measurements near the middle cerebral artery in firstpass perfusion MRI. J Cereb Blood Flow Metab 2009;29:840–852.
Chappell et al. 22. Calamante F. Arterial input function in perfusion MRI: a comprehensive review. Prog Nucl Magn Reson Spectrosc 2013;74:1–32. 23. Frangi AF, Niessen WJ, Vincken KL, Viergever MA. Multiscale vessel enhancement filtering. Med Image Comput Comput Assist Interv 1998;1496:130–137. 24. Straka M, Newbould RD, Sramek M, Albers GW, Bammer R. Improving DSC-MRI by orientation-corrected phase-based AIF and VOF. In Proceedings of the 18th Annual Meeting of ISMRM, Stockholm, Sweden, 2010. Abstract 622. 25. Kjolby B, Mikkelsen I, Pedersen M, Ostergaard L, Kiselev V. Analysis of partial volume effects on arterial input functions using gradient echo: a simulation study. Magn Reson Med 2009;61:1300–1309. 26. 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:514–523. €stergaard L. Automatic 27. Mouridsen K, Christensen S, Gyldensted L, y selection of arterial input function using cluster analysis. Magn Reson Med 2006;55:524–531. 28. Chappell MA, Woolrich MW, Petersen ET, Golay X, Payne SJ. Comparing model-based and model-free analysis methods for QUASAR arterial spin labeling perfusion quantification. Magn Reson Med 2013;69:1466–1475. 29. Lee JJ, Bretthorst GL, Derdeyn CP, Powers WJ, Videen TO, Snyder AZ, Markham J, Shimony JS. Dynamic susceptibility contrast MRI with localized arterial input functions. Magn Reson Med 2010;63: 1305–1314. 30. Schmiedeskamp H, Straka M, Newbould RD, Zaharchuk G, Andre JB, Olivot J-M, Moseley ME, Albers GW, Bammer R. Combined spin- and gradient-echo perfusion-weighted imaging. Magn Reson Med 2012;68: 30–40. 31. Weisskoff R, Zuo C, Boxerman J, Rosen B. Microscopic susceptibility variation and transverse relaxation: theory and experiment. Magn Reson Med 1994;31:601–610. 32. Birn RM, Bandettini PA. The effect of T2’ changes on spin-echo EPIderived brain activation maps. In Proceedings of the 10th Scientific Meeting of ISMRM, Honolulu, Hawaii, USA, 2002. p. 1324.
