Improved Leakage Correction for Single-Echo Dynamic Susceptibility Contrast Perfusion MRI Estimates of Relative Cerebral Blood Volume in High-Grade Gliomas by Accounting for Bidirectional Contrast Agent Exchange

The authors' hypothesis is that incorporating bidirectional contrast agent transport into the DSC MR imaging signal model will improve rCBV estimates in brain tumors. A unidirectional contrast agent extravasation model (Boxerman-Weisskoff) was compared with a bidirectional contrast agent exchange model. For both models, they compared the goodness of fit with the parent leakage-contaminated relaxation rate curves and the difference between modeled interstitial relaxation rate curves and dynamic contrast-enhanced MR imaging in 21 patients with glioblastoma. The authors conclude that the bidirectional model more accurately corrects for the T1 or T2* enhancement arising from contrast agent extravasation due to blood-brain barrier disruption in high-grade gliomas by incorporating interstitial washout rates into the DSC MR imaging relaxation rate model. BACKGROUND AND PURPOSE: Contrast agent extravasation through a disrupted blood-brain barrier potentiates inaccurate DSC MR imaging estimation of relative CBV. We explored whether incorporation of an interstitial washout rate in a leakage-correction model for single-echo, gradient-echo DSC MR imaging improves relative CBV estimates in high-grade gliomas. MATERIALS AND METHODS: We modified the traditional model-based postprocessing leakage-correction algorithm, assuming unidirectional contrast agent extravasation (Boxerman-Weisskoff model) to account for bidirectional contrast agent exchange between intra- and extravascular spaces (bidirectional model). For both models, we compared the goodness of fit with the parent leakage-contaminated relaxation rate curves by using the Akaike Information Criterion and the difference between modeled interstitial relaxation rate curves and dynamic contrast-enhanced MR imaging by using Euclidean distance in 21 patients with glioblastoma multiforme. RESULTS: The bidirectional model had improved Akaike Information Criterion versus the bidirectional model in >50% of enhancing tumor voxels in all 21 glioblastoma multiformes (77% ± 9%; P < .0001) and had reduced the Euclidean distance in >50% of enhancing tumor voxels for 17/21 glioblastoma multiformes (62% ± 17%; P = .0041). The bidirectional model and dynamic contrast-enhanced-derived kep demonstrated a strong correlation (r = 0.74 ± 0.13). On average, enhancing tumor relative CBV for the Boxerman-Weisskoff model exceeded that for the bidirectional model by 16.6% ± 14.0%. CONCLUSIONS: Inclusion of the bidirectional exchange in leakage-correction models for single-echo DSC MR imaging improves the model fit to leakage-contaminated DSC MR imaging data and significantly improves the estimation of relative CBV in high-grade gliomas.

T he most common DSC MR imaging metric in neuro-oncology is relative CBV (rCBV), 1 which has been used for grading gliomas, 2,3 predicting low-grade to high-grade transformation, 4,5 distinguishing recurrent tumor from pseudoprogression, 6,7 differentiating tumor regression from pseudoresponse, 8 and assessing overall treatment response. 9,10 Relative CBV is typically calculated by integrating the dynamic first-pass change in the transverse relaxation rate (⌬R 2 * ) resulting from bolus injection of a gadolinium-based contrast agent, which transiently causes a dose-dependent change in the magnetic susceptibility of blood. 11 This technique mimics the classic indicator-dilution theory, 12 which assumes intravascular compartmentalization of injected contrast agent "tracer." However, common gadolinium-based contrast agents extravasate in lesions with blood-brain barrier disruption, 13 including malignant gliomas. The exchange of con-trast agent between the intravascular and the extravascular extracellular space, which is the objective measurement in dynamic contrast-enhanced (DCE) MR imaging, [14][15][16] contaminates the desired DSC MR imaging signal, depending on pulse sequence parameters and underlying tumor biology. 17 A popular model-based DSC MR imaging leakage-correction method proposed by Weisskoff and Boxerman 2,18,19 linearly fits measured ⌬R 2 * (t) to 2 constant functions derived from the average relaxation rate in nonenhancing tissue, one of which is permeability-weighted. Deviation from the reference function is used to derive corrected rCBV for each voxel. A limiting assumption of this approach is that contrast agent reflux from the interstitial space back to blood plasma is negligible within the time frame of DSC MR imaging signal acquisition (ϳ2 minutes). However, standard models quantifying contrast agent exchange between blood plasma and the interstitium (ie, DCE MR imaging 14 ) use 2-compartment pharmacokinetics to account for bidirectional transport of contrast agent. We hypothesized that incorporating bidirectional contrast agent transport into the original DSC MR imaging signal model improves rCBV estimates in brain tumors. To test this hypothesis, we compared model-based DSC MR imaging leakage-correction methods with and without consideration of bidirectional transport by using simulations and clinical application to high-grade gliomas.

Patients
We studied 24 sequential patients with histologically proved glioblastoma multiforme (GBM) treated with maximal surgical resection followed by radiation therapy and concurrent temozolomide and both DSC MR imaging and DCE MR imaging performed at initial tumor progression.  20 Tumor sizes ranged from 2.8 to 106.6 mL, with an average enhancing volume of 40.1 Ϯ 28.4 mL. Spheric ROIs of 1.6 mL were also selected in normal-appearing, contralateral white matter for rCBV normalization.

Computation of DSC MR Imaging rCBV
All simulations and calculations were performed in Matlab (MathWorks, Natick, Massachusetts) by using custom scripts. Uncorrected rCBV was calculated from trapezoidal integration of the original DSC MR imaging relaxation rate-time curve, ⌬R * 2 ͑t͒. The whole-brain average relaxation rate for nonenhancing voxels (Equations 3 and 4, all equations are in the Appendix) was used for both the original Boxerman-Weisskoff model 19 (BW model) and the new bidirectional exchange model (bidir model). Linear least-squares optimization was used to determine the free parameters for both the bidir-model (via Equation 7) and the BW model (Equation 5, with k ep ϭ 0) algorithms, and the corrected rCBV was computed from Equation 8. The average run-time per patient in Matlab was 19.5 Ϯ 6.7 seconds for the bidir model and 18.3 Ϯ 6.2 seconds for the BW model (3.2-GHz Intel Core i5, 32 GB RAM). Tumor rCBV for each method was subsequently normalized to median rCBV within the normal-appearing white matter ROI.

Simulation of DSC MR Imaging rCBV
The whole-brain average relaxation rate, ⌬R * 2 ͑t͒, was chosen from a sample patient and corresponds to the curve with K 1 ϭ 1, K 2 ϭ 0, and k ep ϭ 0. K 2 ϭ 0.05 (adding T1-dominant leakage) with k ep ϭ 0 was set to simulate the BW model. A nonzero k ep (0.002 or 0.005) was used to simulate the bidir model of ⌬R* 2 ͑t͒. For k ep ϭ 0.1, the simulation is reflective of the correction of relaxation rate curves at "arterylike" voxels.

Goodness of Fit Analysis
For each enhancing tumor voxel for all patients, we computed the Akaike Information Criterion (AIC) between the leakage-contaminated relaxation rate ⌬R * 2 ͑t͒ (Equation 1) and its model fit (Equation 5) for the BW model and bidir-model: where n is the number of fitted time points (injection to the end of the DSC MR imaging acquisition), RSS is the sum of the squared residuals, and p is the number of free parameters (2 for the BW model, 3 for the bidir-model). 21 Differences in the BW model and bidir model AIC were calculated for all voxels with k ep Ͼ 0. We also computed the Euclidean distance (square root of the sum of the squared differences) between the interstitial leakage relaxation rate curves, ⌬R* 2,E ͑t͒, generated by the BW model and bidir model corrections and the DCE MR imaging signal, in which the DCE MR imaging signal was upsampled from a 5-second resolution to a 1.8-second resolution to match that of the DSC MR imaging data via linear interpolation by using the Matlab function "resample." Because interstitial leakage relaxation rate curves and DCE MR imaging signals have units of 1/s and mM, respectively, both were standardized to an area under the curve equal to unity and were vectorized for computation of the Euclidean distance. Higher AIC and Euclidean distance imply worse fits. Two-sample t tests were used to compare whether the AIC and Euclidean distance measurements were significantly different between the 2 leakage-correction methods.

Postprocessing of DCE MR Imaging
DCE MR imaging biomarkers, k ep and contrast transfer coefficient (K trans) , were derived via a fit to the model of Tofts and Kermode. 14 As described, the temporal resolution of the DCE MR imaging data was upsampled to match the DSC MR imaging data. For the DCE MR imaging analysis, the "whole-brain average" served as the arterial input function for the DCE model fit. This was done to mirror the DSC bidir model analysis, in which the "whole-brain average" effectively serves as the arterial input function. Voxels with highly fluctuating time courses in either the DSC or DCE images were eliminated from the analysis.

Correlation between DSC-and DCE-Derived Imaging Biomarkers
DSC MR imaging biomarkers, k ep and rCBV, were derived as described in the Appendix. Voxelwise Pearson correlation coefficients between the DSC-and DCEderived parameters were performed in Matlab within contrast-enhancing tumor only, for each patient independently. In this study, we report means and SDs of the correlation coefficients from all 21 patients. Figure 1 compares the simulated total leakage contaminated relaxation rate, ⌬R* 2 ͑t͒, (Fig 1A) and the component from interstitial leakage, ⌬R* 2,E ͑t͒, (Fig 1B) for various conditions according to the Tofts and Kermode model, 14 assuming T1-dominant leakage-associated relaxation enhancement. For the BW model, ⌬R* 2,E ͑t͒ rises with time in the absence of washout. For nonzero k ep , there is less rise in ⌬R* 2,E ͑t͒ and closer approximation of the tail of ⌬R* 2 ͑t͒ to ⌬R * 2 ͑t͒, reflecting tumors with different contrast agent pharmacokinetics. For k ep ϭ 0.1, the tail of ⌬R * 2,E ͑t͒ approaches zero, but because the first-pass of ⌬R * 2 ͑t͒ differs from that of ⌬R * 2 ͑t͒, correction of relaxation rate curves at "arterylike" voxels by using K 1 and K 2 is still required to achieve accurate rCBV estimates. Figure 1C plots sample ⌬R* 2 ͑t͒, with T2*-dominant leakageassociated relaxation enhancement for a representative patient, with superimposed BW model and bidir model fit relaxation rate curves. In this example, the BW model overestimates the first-pass curve, underestimates the second and third passes, and overestimates the tail. The bidir model better approximates ⌬R * 2 ͑t͒ over all time points, visually, and has substantially improved the AIC, quantitating an improved fit to the total leakage-contaminated relaxation rate curve. Figure 1D plots standardized DCE MR imaging signal for the tumor voxel used in Fig 1C, with superimposed standardized interstitial leakage relaxation rate curves, ⌬R* 2,E ͑t͒, from the BW model and bidir model. The standardized interstitial leakage re- laxation rate continually rises with time for the BW model, whereas it better tracks standardized DCE MR imaging for the bidir model, with a substantially improved Euclidean distance. Figure 2 plots the percentage of voxels in which the bidir model outperformed the BW model for AIC and Euclidean distance metrics in whole brain and tumor for the 21 patients with GBM. The bidir model had better AIC performance than the BW model in Ͼ50% of whole-brain (mean, 71% Ϯ 6%, P Ͻ .0001) and tumor (mean, 77% Ϯ 9%, P Ͻ .0001) voxels in all patients, and better Euclidean distance performance in Ͼ50% of whole-brain voxels (mean, 80% Ϯ 9%, P Ͻ .0001) for all patients and in tumor voxels (mean, 62% Ϯ 17%, P ϭ .0041) for 17 of the 21 patients. All were statistically significant for a 1-sample t test with null hypothesis of 50%.

Correlation between DSC-and DCE-Derived Imaging Biomarkers
We then performed a voxelwise correlation between the DSCderived imaging biomarkers from the bidirectional leakage-correction algorithm (k ep and rCBV) with the DCE-derived imaging biomarkers (k ep and K trans ). The Pearson correlation coefficient between the 2 k ep measurements was 0.74 Ϯ 0.13 across the 21 patients, with a weak correlation between the Pearson correlation coefficient and tumor size (r ϭ 0.11). Figure 3 demonstrates an example of the correlation between DSC-and DCE-derived k ep . A correlation test was performed between the bidirectional modelderived rCBV and DCE-derived K trans , with a moderate correlation of 0.49 Ϯ 0.22. A moderate correlation was also found between rCBV and plasma volume fraction (vp) at 0.54 Ϯ 0.12. Finally, the correlation between the same rCBV and k ep was r ϭ 0.29 Ϯ 0.26. The average K trans value was 0.0015 Ϯ 0.0018 seconds Ϫ1 (0.09 Ϯ 0.11 minutes Ϫ1 ), DCE K ep was 0.0050 Ϯ 0.0023 seconds Ϫ1 (0.30 Ϯ 0.14 minutes Ϫ1 ), DSC k ep was 0.0057 Ϯ 0.0042 seconds Ϫ1 (0.34 Ϯ 0.25 minutes Ϫ1 ), vp was 0.01 Ϯ 0.01, and rCBV was 1.98 Ϯ 1.24. A closer inspection of the T2*-dominant-versus-T1-dominant voxels (as defined by a negative or positive K 2 , respectively) revealed that the difference between the 2 correction methods in T2*-dominant voxels was 37.7% Ϯ 42.6%, while the same metric for T1-dominant voxels was 5.8% Ϯ 3.4%.

DISCUSSION
By incorporating the Tofts and Kermode model into the singleecho DSC MR imaging relaxation rate equation, we developed an improved postprocessing leakage-correction method accounting for bidirectional contrast agent transport between the intravascular and interstitial spaces that commonly occurs in angiogenic  high-grade gliomas. Our results demonstrate the importance of considering the interstitial washout term, even when modeling the relaxation rate changes during short image acquisitions. For instance, in the simulation, we observed differences between the bidir model and the BW model fits to relaxation rate data in highgrade gliomas in the first-pass curve (as early as 10 -20 seconds after injection). Furthermore, inclusion of a washout term in the bidir model alleviates the error in relaxation rate estimates for arteries and normal brain introduced by conventional models constrained to increasing contrast agent concentration with time in all tissues.
Our results suggest that the conventional BW model undercorrects rCBV, with insufficiently increased and decreased rCBV compared with uncorrected rCBV in T1-dominant and T2*dominant leakage scenarios, respectively. Furthermore, because the low flip angle DSC MR imaging protocol was largely T2*dominant and the largest discrepancies between the bidir model and BW model estimates of rCBV existed for T2* dominant voxels, our results suggest that the bidir model may be particularly advantageous over the BW model for correcting the residual T2* effects frequently encountered in dual-echo gradient-echo acquisitions. This algorithm can be performed without a substantial increase in postprocessing computation time over the unidirectional model; therefore, the bidirectional model can simply replace the previous model in routine clinical work and for evaluating tumor grade, distinguishing pseudoprogression from true progression, and evaluating treatment response.
Several postprocessing leakage-correction techniques have previously been proposed. 22,23 The method by Boxerman-Weisskoff, 2,18,19 which linearly fits measured ⌬R * 2 ͑t͒ to 2 constant functions derived from the average relaxation rate in nonenhancing tissue, can be applied quickly to conventional single-echo (spin-echo or gradient-echo) acquisitions and contrast agent injection schemes. Improved correlation of rCBV with glioma grade compared with uncorrected rCBV 19 provides anecdotal evidence of the benefit of the BW model, which has also been shown to improve correlation of gadolinium-based rCBV measures over those obtained by using the intravascular magnetic iron oxide nanoparticles agent as a criterion standard. 24 Bjornerud et al 25 proposed a method that reduces the sensitivity of rCBV correction to mean transit time that could be combined with the bidir model scheme. Most interesting, Schmiedeskamp et al 23 used a multiecho, gradient-echo, spin-echo acquisition scheme to correct for T1 and T2* leakage by using a backflow term; however, results were highly dependent on literature values for r* 2,E and r* 2,P , the T2* relaxation effects of gadolinium in the extravascular space and plasma, respectively, which can vary quite substantially depending on the literature source. Additionally, Quarles et al 17 suggested that these values could vary from tumor to tumor, depending on physiologic factors such as interstitial, vascular, and cell volume fractions and vessel and cell size. An advantage of the bidir model correction method is the lack of assumptions for r* 2,E and r* 2,P . All of these leakage-correction algorithms aim to isolate the relaxation rate due to the residual intravascular contrast agent by eliminating the T1-and T2*-related contributions to the relaxation rate from the extravasated contrast agent. They do not "add back" T2* relaxation that would have been realized had the extravasated contrast agent not left the plasma space, so "corrected rCBV" may still differ from that computed for a tumor with no vascular permeability, all other parameters (including true blood volume) being equal.
One potential limitation to this study is its retrospective design, which may have yielded a selection bias in the sample. Specifically, all patients were chosen because they failed standard therapy. Another potential limitation is the lack of correlation with a criterion standard, such as histology, or with CBV estimates by using intravascular agents such as iron oxide contrast agents. Moreover, AIC is a unitless quantity, which can compare relative goodness of fit between models but does not have a direct test to determine whether one model is significantly better than the other. Finally, the current study only included patients with glioblastoma; therefore, we were unable to recommend a threshold between low-grade and high-grade gliomas by using the new leakage-correction algorithm.

CONCLUSIONS
The bidir model more accurately corrects for the T1 or T2* enhancement arising from contrast agent extravasation due to blood-brain barrier disruption in high-grade gliomas by incorporating interstitial washout rates into the DSC MR imaging relaxation rate model. To this end, the bidir model may potentially improve patient diagnosis and evaluation of treatment response by more accurately estimating rCBV in DSC MR imaging.

APPENDIX
Following Equation A6 of Boxerman et al, 19 the leakage-contaminated DSC MR imaging relaxation rate-time curve, ⌬R* 2 ͑t͒, equals the intravascular contrast-driven transverse relaxation rate change, ⌬R* 2 ͑t͒ plus ⌬R* 2,E ͑t͒, a tissue-leakage term describing the simultaneous T1 and T2* relaxation effects resulting from gadolinium extravasation: where E 1 ϭ e ϪTR/T 1o , T 1o is the precontrast tissue T 1 , r 1 is the T 1 relaxivity of gadolinium, C E (t) is the concentration of gadolinium in the extravascular extracellular space, and r* 2,E represents the T2* relaxation effects of gadolinium extravasation, as described by Quarles et al 17 and Schmiedeskamp et al. 23 From the original Tofts and Kermode model describing bidirectional contrast agent flux between the intravascular and extravascular compartments, 14 we can estimate the concentration in the extravascular space as: 2) C E ͑t͒ ϭ k trans ͓C p ͑t͒ ϫ e Ϫkept ͔, where k trans and k ep are the transfer coefficients for intra-to extravascular and extra-to intravascular contrast flux, respectively, and C p (t) is the plasma contrast concentration. C p (t) and ⌬R 2 *(t) can be defined as scaled versions of the whole-brain average relaxation rate in nonenhancing voxels, ⌬R * 2 ͑t͒ 19 : 3) C p ͑t͒ ϭ k ϫ ⌬R * 2 ͑t͒ 4) ⌬R* 2 ͑t͒ ϭ K 1 ϫ ⌬R * 2 ͑t͒.
Combining Equations 1-4 yields the following: K 1 , K 2 , and k ep (units of second Ϫ1 ) are the free parameters of Equation 5. In general, K 1 depends on CBV, vessel size, and other physiologic factors, while K 2 is related to vascular permeability. Substituting k ep ϭ 0, which occurs with no backflow of extravasated contrast agent, yields the original Boxerman-Weisskoff leakage-correction algorithm, where K 1 and K 2 are solved by linear least-squares fit to ⌬R* 2 ͑t͒. 19 For the bidir model correction method, a linear least-squares fit to K 1 , K 2 , and k ep can be used with the methodology of Murase, 26 as described by the following equation: Integrating the corrected relaxation rate-time curve yields the following expression for leakage-corrected rCBV: 8) rCBV corr ϭ rCBV ϩ K 2͵ 0 T ͵ 0 t ⌬R * 2 ͑͒ ϫ e Ϫkep͑t Ϫ ͒ ddt.