Blood–brain barrier permeability of normal‐appearing white matter in patients with vestibular schwannoma: A new hybrid approach for analysis of T 1‐W DCE‐MRI

Purpose To develop and assess a “hybrid” method that combines a first‐pass analytical approach and the Patlak plot (PP) to improve assessment of low blood–brain barrier permeability from dynamic contrast‐enhanced (DCE) magnetic resonance imaging (MRI) data. Materials and Methods Seven patients with vestibular schwannoma were enrolled. T 1‐W DCE imaging was acquired on a 1.5T scanner. Normal‐appearing white matter (NAWM) was divided into four regions of interest (ROIs) based on the magnitude of changes in longitudinal relaxation rate (ΔR1) after gadolinium administration. Kinetic analysis of ROI‐averaged contrast agent concentration curves was performed using both the conventional PP and the hybrid method. Computer simulated uptake curves that resemble those from NAWM were analyzed with both methods. Percent deviations (PD) of the “measured” values from the “true” values were calculated to evaluate accuracy and precision of the two methods. Results The simulation showed that, at a noise level of 4% (a noise level similar to the in vivo data) and using a signal intensity (SI) averaging scheme, the new hybrid method achieved a PD of 0.9 ± 2.7% for v p, and a PD of –5.4 ± 5.9% for K trans. In comparison, the PP method obtained a PD of 3.6 ± 11.3% for v p, and –8.3 ± 12.8% for K trans. One‐way analyses of variance (ANOVAs) showed significant variations from the four WM regions (P < 10−15 for ΔR1; P < 10−6 for K trans; P < 10−4 for v p). Conclusion Both computer simulation and in vivo studies demonstrate improved reliability in v p and K trans estimates with the hybrid method. Level of Evidence: 3 Technical Efficacy: Stage 1 J. MAGN. RESON. IMAGING 2017;46:79–93

contrast concentration time course over the first-pass (FP) transit period of the contrast agent (CA) bolus. 3 This might reflect the weakness of the Patlak method, which depends on linear regression analysis, and is more prone to nonuniform (distorted) noise. In addition, the model does not allow for backflux of CA from the tumor into the plasma, which will be of particular importance if an inappropriate time interval is chosen for analysis.
In this study we propose a new "hybrid" method that combines a first-pass analytical approach 7,8 with the Patlak plot (PP) to improve assessment of BBB permeability. The PP model describes a unidirectional two-compartment system and uses linear regression analysis to estimate K trans and v p . The FP method 7,8 performs an automatic decomposition of the first-pass CA concentration curve into intravascular and interstitial components to allow simultaneous mapping K trans and v p . A leakage-corrected estimate of v p is obtained by integrating the area under the intravascular CA concentration curve over the first-pass of the bolus. When combined with high temporal resolution DCE data, the FP method provides accurate and robust measurements of v p . 2,9 The aim of this study was to develop an easy-to-use and reliable approach for detecting subtle BBB permeability in normal-appearing white matter (NAWM) and to apply the proposed techniques in patients with type II neurofibromatosisassociated vestibular schwannoma (VS). We hypothesize that a hybrid approach combining the benefits of the FP and PP analytical approaches will improve the accuracy of parameter estimates in tissues with low permeability.

Patlak Plot
The Patlak model describes a highly perfused twocompartment tissue assuming unidirectional transport from the plasma into the extravascular extracellular space (EES). The CA concentration in tissue is given by: where C p (t) is the plasma CA concentration time course curve, v p is the fractional plasma volume, and K trans is the volume transfer constant between blood plasma and the leakage space. Dividing both sides of Eq. [1] with C p (t), one obtains: C t ðtÞ C p ðtÞ 5v p 1K trans ð t 0 C p ðsÞd s C p ðtÞ : Equation [2] expresses the PP, where the slope represents K trans and the intercept represents v p . The term on the left side of the equation, C t (t)/C p (t), represents the volume of distribution (v d ) of the CA in brain tissue at the time of sampling, t. 10 If the BBB is intact, then v d should be equivalent to v p . If leakage occurs, then v d becomes larger than the v p , as accessible extravascular compartments are included. 2 The abscissa has the units of time, but this is not laboratory time. It is concentration-stretched time and will be referred to hereafter as t stretch . 3 Figure 1 shows the relationship between the lab time and t stretch calculated using a C p (t) measured from low CA dose, high temporal resolution (LDHT) imaging in a patient with VS. 11 It can be seen that a t stretch interval 80-250 seconds 2 corresponds to a lab time interval 37-157 seconds. The relationship between t stretch and lab time is approximately linear beyond a lab time of $37 seconds, but fluctuates within the first-pass and recirculation phase of the CA bolus. In most studies permeability is evaluated during the steady-state component of the plot, ignoring the initial bolus circulation. An additional problem is that experimental errors are distorted when a nonlinear model is transformed to a linear one. 12 How these distorted errors and the initial timepoints affect the Patlak fitting will be investigated in the current study.
FP Model for Simultaneously Deriving v p and K trans (FP simul ) The FP method 7,8 assumes that backflow during the firstpass transit period of the bolus is negligible, and uses an iterative approach 8 to separate the intravascular and interstitial components of the tissue CA concentration curve, which are then used for v p and K trans estimation, respectively (see details in Appendix A).

New Hybrid Method
The new hybrid method contains three parts: 1. a conventional PP for simultaneously deriving K trans and v p ; 2. a modified FP analysis with a known K trans for deriving v p only; 3. a modified PP analysis with a known v p for deriving K trans only. A flow chart of the new hybrid method is presented in Fig. 2. The analysis is performed in three steps: Step 1: CA time course curves are fitted using the conventional PP to derive an initial estimate of K trans and v p values, with a t stretch interval of 85-250 sec used in the fitting; Step 2: The estimate of K trans from Step 1 is used into Eq. [A.4] (Appendix A) as a known value and v p is then calculated using Eqs. [A.5] and [A.6] (Appendix A); Step 3: CA time course curves are then refitted with Eq.
[2] to obtain a refined K trans estimate. During the fitting, the v p in Eq. [2] is fixed using the value obtained in Step 2, leaving K trans as the only free-fitting parameter.
The estimates of v p from Step 2, and K trans from Step 3 are the final output estimates from the analysis.
The conventional PP performed at Step 1 provides an adequate estimate of K trans for use as the initial value 3 in the modified FP method enabling subsequent calculation of v p , corrected for leakage. In low permeability tissues such as NAWM the intravascular CA concentration is very high compared to the amount of transendothelial leakage. This means that errors in K trans estimation from Step 1 will have only very small effects on the calculation of the integral of the intravascular CA time course curve in Step 2. Consequently, the modified FP approach becomes less sensitive to covariance errors that result from simultaneous measurement of K trans and v p . The resulting improved estimate of v p is then used in the modified PP method, leaving K trans as the only free-fitting parameter. The stepwise estimation of K trans and v p may be expected to improve measurement accuracy and enable more reliable estimation in tissue with very low levels of endothelial permeability.

Patients
Seven patients with type 2 neurofibromatosis (NF2), with a total of 20 tumors (14 VS and 6 meningiomas), were recruited into the study. DCE-MRI studies were performed in the absence of any treatment. All subjects gave informed consent and the Local Research Ethics Committee approved the study (reference number O2-051).

MRI
Patients were imaged on a 1.5T whole body scanner (Philips Achieva, Philips Medical Systems, Best, Netherlands) using an 8-channel head coil. LDHT DCE-MRI data were collected as part of the dual temporal resolution technique, ICR-DICE, as described previously. 11 Prior to the LDHT DCE series, four consecutive 3D axial fast gradient recalled echo acquisitions (GRE) with variable flip angle (VFA; a 5 28, 88, 158, and 208) were performed for native longitudinal relaxation rate (R1 N 5 1/T1 N ) mapping. The fourth sequence was then repeated (n 5 300) to produce a high temporal resolution (Dt 5 1.03 sec) T 1 -W dynamic dataset with a low dose of CA (gadoterate meglumine; Dotarem, Geurbet, Roissy, France) given following the 30 th dynamic scan. Contrast agent (a fixed volume of 3 ml, $0.02 mmol/kg depending on body weight) was administered by power injector as an intravenous bolus at a rate of 3 ml/s, followed by a chaser of 20 ml of 0.9% saline administered at the same rate. Two sets of high spatial resolution T 1 -W 3D images with the voxel size of 1 3 1 3 1 mm, ie, isotropic voxels, were acquired before (T 1 -weighted [W] without contrast) and after the DCE MRI (T 1 -W 1 contrast).

Image Processing
Prior to kinetic analysis of the CA concentration curves observed in NAWM, we performed pixel-by-pixel mapping of postinjection changes in longitudinal relaxation rate (DR1) to show CA distribution due to BBB leakage in the brain tissue. A four-color scheme was applied to the WM-segmented DR1 maps, which divided the WM into four ROIs based on the magnitude of DR1. Kinetic analysis of ROI-averaged CA concentration curves was then performed using both the conventional PP and the hybrid method. Computer simulations were performed to generate CA concentration curves that resemble the in vivo CA curves to evaluate accuracy (bias) and precision (standard deviation) of the new hybrid method.
MEASURING PLASMA CA CONCENTRATION TIME COURSE CURVE (C P (t)). The plasma CA concentration time course curve (C p (t)) was measured from the superior sagittal sinus (SSS) in patient data as shown in Fig. 1. A semiautomatic extraction technique was used in the C p (t) measurement (see Appendix B for details). To reduce the influence of Rician noise, the average of the signal intensity (SI) in the last 10 dynamic frames (around 4.5 min postinjection) was used to calculate the postinjection tissue R1 (R1 post ). Calculation of R1 post based on R1 N , mean preinjection signal intensity (SI pre ), and postinjection signal intensity (SI post ), with a subtraction method (SI post -SI pre ), is described in Appendix D.
Maps of DR1 ( 5 R1 post -R1 N ) were generated. To view the spatial distribution of the residual CA in brain tissues, WMsegmented DR1 maps were displayed in four colors (blue, green, red, and yellow) based on the magnitude of DR1, and associated v p values. The criteria setting for the color-coding (blue: DR1 < 0 or v p < 0.01; green: 0 < DR1 < 0.012 and v p > 0.01; red: 0.012 DR1 < 0.025 and v p > 0.01; and yellow: DR1 ! 0.025 and v p > 0.01; v p was estimated by integrating the area under the firstpass CA concentration curves without leakage correction) was based on in vivo observation from a patient with NF2, and will be explained in more detail in the Results section below.
SEGMENTATION OF GRAY AND WHITE MATTER. SPM2 13 was used for 1) spatial alignment between R1 N VFA, DCE-MRI, and 3D T 1 -W isotropic images, and 2) segmentation of the MRI data into GM, WM, and CSF. The probability maps of GM, WM, and CSF segmented from the T 1 -W isotropic images were realigned and resliced to the space of the 3D individual frames of the DCE-MRI, as well as the 3D R1 N and 3D pharmacokinetic parametric images, ie, K trans and v p . WM masks were generated from the WM probability maps by including only voxels with a probability greater than 0.95 and were used for the subsequent quantitative analysis.

QUANTITATION OF TISSUE BBB PERMEABILITY WITH THE
NEW HYBRID METHOD. ROI-averaged (SI) time course data were used to test the new hybrid method. Each ROI was made as a collection of all the pixels with same color in a WM-segmented image. Both conventional Patlak and the new hybrid methods were applied and compared. In addition, three t stretch intervals (85-250 sec, 85-300 sec, and 0-250 sec) were used in the fitting to evaluate how much the specific time intervals could affect the results.

Computer Simulation to Evaluate the New Hybrid Method
Based on the results from the in vivo data analysis, K trans 5 0.0074 min 21 and v p 5 0.024 were chosen as "true" values to synthesize tissue CA concentration curves, which resemble the in vivo CA uptake curve in the yellow-coded region. The C p (t) shown in Fig. 1 was used in the CA concentration curve simulation and fitting.
To investigate the effects on parameter estimation of the assumption that no backflux of CA occurs, zero noise tissue uptake curves were synthesized with the modified Tofts model 14,15 (assuming the fractional volume of the extravascular extracellular space, v e 5 0.20) and the unidirectional two-compartment model, respectively. The modified Tofts model with no backflux corresponds to the unidirectional two-compartment model expressed as Eq. [1]. Fitting errors due to ignoring backflux when using the PP and the new hybrid method were compared.
To investigate the effects of noise on parameter estimation, CA concentration curves were simulated with the modified Tofts model and converted into an SI-time curve based on the in vivo mean baseline SI (470; 33 baseline frames); the precontrast T 1 relaxation time (T 10 ), and a literature value of longitudinal relaxivity (4.39 mM 21 sec 21 ). 16 The generated SI-time courses were sampled with a temporal resolution of 1.03 seconds. Rician white noise with noise level (ie, standard deviation / mean baseline signal) of 1%, 2%, 3%, 4%, and 5%, respectively, was added to the simulated SI-time curves. K trans and v p were calculated using the synthetic datasets to produce the so-called "measured" values. Percentage deviations (PD) of the "measured" values from the "true" values were calculated as: PD 5 (measured -true)/true. Both the conventional PP and the new hybrid method were used for kinetic analysis. A total of 20,000 repetitions were performed for each method to produce mean and standard deviation (SD) of PD for each parameter estimates.
To investigate the effects of averaging SI curves on parameter estimation while using the two kinetic analysis methods, 100 individual SI-time curves, simulated as described above, were averaged to resemble the ROI-averaged SI curves observed from the in vivo WM yellow-coded region. 200 repetitions were performed for each method to produce mean and SD of PD for each parameter estimates.
A t stretch interval 85-250 seconds was used in the above Monte Carlo simulations. To investigate the effects of including the initial timepoints in the Patlak fitting, we repeated the above simulations using a t stretch interval 0-250 seconds instead. We also repeated the simulations setting K trans at 0.004, 0.008, 0.012, 0.016, 0.020, 0.025, 0.030, 0.035 min 21 , respectively (other parameters were fixed at v e 5 0.20, v p 5 0.024, and noise level of 4%), to identify what level of back-diffusion (k ep 5 K trans /v e ) will cause the new hybrid method to fail, as indicated by high absolute value of PD mean or high SD of PD.

Further In Vivo Study
CORRELATION BETWEEN DR1 AND K TRANS IN WM. Mean DR1 of each of the four color-coded WM regions were calculated for each slice in the DR1 image volumes. The unidirectional influx constant K trans and v p were derived from the corresponding regional uptake curves using the new hybrid method. The whole t stretch interval (from 0 to $460 sec) was used in the fitting, which was in agreement with the lab time interval used for calculation of DR1. Using longer t stretch interval also benefits the measurement of very low BBB leakage (the red and green-coded WM regions). Linear regression analysis was performed to evaluate the relationship of DR1 and the unidirectional influx constant K trans .

COMPARISON OF THE FOUR COLOR-CODED REGIONS IN
WM. Mean DR1, K trans and v p measured in each of the four color-coded WM regions were compared. Data are expressed as mean 6 SD for the seven patients. Statistical tests were performed with 95% confidence intervals, using one-way analysis of variance (ANOVA) to assess variations of the parameters in relation to the four color-coded WM regions. Post-hoc Tukey's honestly significant difference (HSD) test was performed to determine statistically significant differences among the mean parameter values for each pair of the WM regions.

ROBUST STABILITY AGAINST VARIATION IN THE T STRETCH
INTERVAL. For assessing the degree of variation in K trans and v p while using various t stretch intervals in the kinetic analysis, the coefficient of variation (CoV), defined as the ratio of the standard deviation to the mean, was computed. Robust stability against variation in the t stretch interval was compared between the two kinetic methods.

Statistical Analysis
The accuracy and precision of the two methods were compared by computer simulation of the mean and SD of percentage deviations for K trans and v p estimates assuming the existence of backflux of CA and at varying levels of Rician noise. Stability of the two methods to variation in the t stretch interval was evaluated by computing the CoV in K trans and v p estimates derived from the in vivo data using a range of t stretch intervals. Linear regression analysis was performed to assess the relationship between CA-induced R1 changes and the unidirectional influx constant K trans in WM. Multiple comparisons of the regional means of DR1, K trans , and v p in the four color-coded WM regions were performed using one-way ANOVA, followed by all pairwise comparisons using Tukey's HSD test (alpha 5 0.05). P < 0.05 was considered statistically significant. Figure 3 shows maps of DR1 in WM measured around 4.5 minutes post-CA injection in a patient with NF2 who had bilateral VSs and a frontal meningioma. The maps in Fig. 3 covers slices 37-54, whereas the right VS occurred in slices 14-24, the left VS in slices 15-17, and the meningioma in slices 57-61. From the DR1 maps it can be seen that the residual CA is not homogeneously distributed in WM.

CA-Induced R1 Changes in WM
It was found that about one-fifth of the segmented WM voxels showed a negative DR1 value, and the WM voxels with negative DR1 tended to show low v p estimated with the FP technique. In the above patient, WM voxels with negative DR1 had v p 0.0116 6 0.0107, while those with positive DR1 had v p 0.0216 6 0.0150; WM voxels with v p < 0.01 had DR1 0.0028 6 0.0141, while those with v p > 0.01 had DR1 0.0161 6 0.0163. In Fig. 3, the blue is the area of negative DR1. To reflect the low v p attribute of the blue-coded region, we included all WM voxels with v p < 0.01 into the blue-coded region.
ROI-averaged SI-time curves, which were obtained from the blue-, green-, red-, and yellow-coded regions, respectively, demonstrated variations in baseline level, peak height of first pass of CA bolus, and enhancement extent (Fig. 4).

Quantitation of Tissue BBB Permeability
With PP and the New Hybrid Method Figure 5 shows the fittings of the four regional tissue CA uptake curves using the two kinetic approaches and three t stretch intervals, respectively. For the t stretch interval 85-250 seconds, the results from PP were quite close to those derived from the hybrid method, although the values of v p estimated with PP were generally higher, and the values of K trans were generally lower, compared with their hybrid counterparts. Comparing results with various t stretch intervals, the hybrid method showed better robust stability against variation in the t stretch interval.
Computer Simulation to Evaluate the New Hybrid Method Figure 6 shows that both PP and the hybrid model yielded parameter estimates equal to the true values when there is neither noise nor backflux (Fig. 6a,b). The Patlak analysis overestimates v p (PD 5 4.2%) and underestimates K trans (PD 5 -8.1%) when there is backflux (Fig. 6c). The impact of the backflux was reduced when the new hybrid method was employed (v p : PD 5 0.0%; K trans : PD 5 -5.4%) (Fig. 6d). Figure 7a-d compares PD (mean and SD) in v p and K trans derived from PP and the hybrid methods, at varying noise levels. The conventional Patlak analysis overestimates v p and underestimates K trans . The hybrid method estimates v p accurately, and also reduced the bias in K trans estimation. With the new method, the SD of PD for both v p and K trans estimates is greatly reduced. In addition, averaging the SItime curves simulated under the same conditions produces further improved the precision of both v p and K trans for both methods. The simulation shows that, at a noise level of 4% (a noise level similar to that in the current in vivo study), it is theoretically possible to reach a PD of 0.9 6 2.7% for v p , and a PD of -5.4 6 5.9% for K trans using the new hybrid method combined with the SI averaging scheme. In comparison, the PP method could obtain a PD of only 3.6 6 11.3% for v p , and a PD of -8.3 6 12.8% for K trans . These PD means are similar to the PD values shown in Fig. 6c,d. This means that, under simulation conditions, the bias (means of PD) is caused by the existence of backflux of CA, and the precision (standard deviations of PD) is dominated by the noise. The new hybrid method is shown to be more robust to effects from both the existence of backflux of CA and nonuniform (distorted) noise.
With a t stretch interval 0-250 sec, the simulations showed remarkable improvement in the conventional Patlak analysis, but little effect on the hybrid method. That means that including the initial data points effectively reduced the effects of distorted errors on the Patlak analysis. The accuracy of both K trans and v p estimates and precision of K trans were approaching those of their hybrid counterparts. The precision of v p , from Patlak, however, remained inferior to the estimate produced by the hybrid analysis. On the other hand, in vivo data fitting with a t stretch interval 0-250 seconds (the right columns in Fig. 5a,b) showed good agreement between the two methods for uptake curves from the red-, green-, and blue-coded regions; however, a strong overestimate of v p and underestimate of K trans were seen for the yellow-coded range using the Patlak method. This implies that some other factors, in addition to noise, affect the initial shape of the in vivo uptake curve, especially for the yellow-coded region, and that the hybrid method is more robust in these conditions. Figure 7e shows that, at a noise level of 4% with SI averaging, the conventional PP has started to yield K trans estimates with negative PD beyond -0.20 and v p estimates with positive PD beyond 0.20 when K trans > 0.02 min 21 (k ep > 0.1 min 21 ). In contrast, the new hybrid method started to yield K trans estimates with negative PD beyond -0.20 only when the "true" K trans > 0.030 min 21 (k ep > 0.15 min 21 ), while PD mean for v p estimates remains less than 0.20 even when K trans 5 0.035 min 21 . The SD of percent deviation for K trans estimates reduces when the "true" value of K trans increases, while the SD of PD for v p estimates was not much influenced by the increase in the "true" value of K trans .
Correlation Between DR1 and K trans in WM Figure 8 shows scatterplots of DR1 and K trans measured from the four color-coded WM regions in a representative image slice for each of the seven patients. Linear regression analysis performed on the positive DR1 (yellow, red, and green) regions showed that DR1 and K trans are well correlated, with R 2 5 0.991 6 0.011 and P 5 0.04 6 0.04 (n 5 7). The blue-coded region contained voxels with negative DR1 and voxels with v p less than 0.01 (their DR1 are generally low but may be positive); these were not included in the regression analysis. The correlation between DR1 and K trans indicated that WM-segmented DR1 maps could be used to distinguish between regions in the WM with different K trans levels. However, the variation in the slope and the intercept of the linear regression among the seven patients (slope 5 0.1740 6 0.0548; intercept 5 -0.0012 6 0.0011) implied that DR1 only could not determine the absolute value of K trans . Table 1 illustrates multiple comparisons of mean values of DR1, K trans , and v p measured from four color-coded WM regions in a representative image slice for each of the seven patients as presented in Fig. 8 (right column). One-way ANOVAs showed significant variations from the four WM regions (P < 10 215 for DR1; P < 10 26 for K trans ; P < 10 24 for v p , respectively). Posthoc Tukey's HSD tests indicated statistical differences between the individual regions, except for the blue and green pair on K trans and the green and red pair on v p . Table 1 quantitatively revealed the regional inhomogeneity of BBB permeability in NAWM. It also revealed a trend of positive correlation between the K trans and v p in NAWM. Table 2 lists the K trans and v p derived from fitting the yellow-or the red-coded regional curves as shown in Fig. 5a,b using four t stretch intervals, respectively. Parameter mean from the four t stretch intervals and corresponding CoV are also listed in Table 2 for the Patlak method, and for the hybrid method, respectively. The green-and blue-coded regions were excluded from the CoV analysis due to existence of the negative K trans values in their fitting results. In agreement with the simulation results, in vivo data also showed the hybrid method superior to the Patlak in the robust stability against variation in the t stretch interval.

Discussion
Low-level BBB permeability has been identified in apparently normal-appearing white and gray matter in a range of diseases including cerebrospinal vessel disease, diabetes, dementia, stroke, multiple sclerosis, and systemic lupus erythematosis, as well as normal aging. 4,[17][18][19][20][21][22][23][24] Recent work has also identified changes in permeability in NAWM following whole brain radiotherapy in patients with cerebral tumors. 25 The increase in BBB permeability induced by whole-brain radiotherapy has been used to improve penetration of therapeutic agents for the treatment of primary central nervous system (CNS) lymphoma and metastatic disease typically depending on simultaneous measurement of cerebrospinal fluid (CSF) and plasma concentrations of the therapeutic agent. 26,27 In patients with glioblastoma and, to a lesser extent, some metastatic cerebral tumors, there is commonly extensive disease spread beyond the enhancing tumor areas identified by conventional imaging. 28 There is therefore a need to effectively deliver drugs across the BBB, which has led to a growing interest in therapeutic manipulation of BBB permeability. 29 The development of techniques for effective measurement of low levels of BBB permeability may provide a potentially valuable imaging biomarker both for prognostic/predictive purposes and for use in clinical trials where BBB permeability manipulation is planned.
We describe a new method, which combines the FP and PP approaches, to measure subtle CA leakage through the BBB using DCE-MRI. Both the FP and PP approaches start from the same equation (Eq. [1]), but use different algorithms for deriving kinetic parameters. The FP approach produces more reliable v p estimates using an integration technique, while the Patlak approach is superior for K trans estimation, collecting more data points beyond the first pass, which is necessary to measure subtle leakage. By combining the advantages of these two techniques, the hybrid method becomes more robust to effects from nonuniform (distorted) noise, variation in t stretch interval choice, and error caused by the existence of backflux of CA.
A wide range of t stretch intervals has been used in PP analysis in previous studies. For example, Larsson et al used intervals of 70-250 seconds and 140-250 seconds in a study with patients with brain tumors and healthy subjects, 2 while Ewing et al performed the linear regression over a t stretch interval of 30 minutes in their study with a rat model. 3 In fact, heterogeneity in brain BBB permeability makes optimization of a t stretch interval choice complicated. Shorter t stretch intervals ensure that the linear regression is performed before major backflux of CA occurs and are appropriate for analysis of higher BBB permeability (yellow-coded WM, eg). On the other hand, the longer t stretch intervals should improve measurements of lower BBB permeability (red-and green-coded WM, eg). A kinetic analysis method with the ability to minimize bias due to backflux and variation in t stretch interval choice, such as we present here, will help in addressing this challenging problem.
This article also proposed the combined use of DR1 mapping and ROI-based pharmacokinetic analysis for BBB permeability analysis in NAWM. The information on spatial distribution of residual CA within NAWM is useful since pixel-by-pixel kinetic analysis of WM tissues remains challenging due to a low signal-to-noise ratio (SNR) in the NAWM uptake curves. The color-coded DR1 maps were used for automatic delineation of ROIs with different DR1 levels. The combined temporal-and spatial-averaging scheme could minimize the effects of low SNR: frameaveraged signal intensities were used to calculate R1 pre and R1 post , producing DR1 maps, and then ROI-averaged tissue uptake curves were used for kinetic analysis. In addition, the   and K trans in NAWM shows that WM-segmented DR1 maps may be used to demonstrate the spatial distribution of relative levels of BBB permeability. The WM-segmented DR1 maps and the pharmacokinetic analysis performed afterwards provided interesting information on the spatially inhomogeneous distribution of BBB permeability in NAWM. We found that a portion of WM voxels showed negative values of DR1 and K trans , generally associated with low v p values (<0.01). The mechanism for the apparent negative DR1 and K trans is not clear. The K trans and v p values we observed in the yellow-or red-coded WM regions were similar to those reported by Larsson et al 2 and Heye et al. 17 In addition, the WM regions with higher BBB permeability tend to have higher v p . These observations have the potential to provide quantitative information about WM microstructure and its changes in pathology.
Low dose (1/5 of the standard Gd-DTPA dose) and high temporal resolution DCE-MRI data were used in this study. Taheri et al used a quarter of the standard Gd-DTPA dose and found that R1 in the vicinity of this reduced dose changed much more rapidly with changes in CA concentration than at higher concentrations. 6 A lower dose of Gd-DTPA entails less possibility of introducing a T Ã 2 effect or truncation of the bolus peak of the arterial input function, 2,30 and also less risk of potential side effects of Gd-DTPA in patients with impairment of renal function. 6,31 Our study supported that a dose as low as only 1/5-1/4 of the standard Gd-DTPA dose was still adequate for performing the measurement of BBB permeability in WM. This would enhance the applicability of the technique.
Our study had some limitations. First, we evaluated the new hybrid method only in a small cohort of patients with NF2. Larger studies will be needed to apply the new method in a larger cohort, including patients and healthy subjects for comparison. Second, the plasma concentration curve we measured from the superior sagittal sinus has the merits of high SNR, distinctive peak of first pass, and is free of artifacts resulting from partial volume error or in-flow signal loss. 32 However, we realize that the ideal place to extract a vessel input function is from the arteries which supply the white matter. Nevertheless, to obtain a reliable measurement of C p (t) from the feeding arteries remains challenging.
In conclusion, this study proposed a new pharmacokinetic analysis method for the measurement of subtle BBB permeability in NAWM. Both computer simulation and in vivo study demonstrated improved reliability in v p and K trans estimates. The combination of the pharmacokinetic analysis with pixel-by-pixel mapping CA-induced T 1 changes provides easyto-apply and reliable imaging methods for evaluation of BBB permeability. The heterogeneous distribution of DR1, v p , and K trans , and the relation between them revealed by MRI, may suggest differences in tissue composition between different locations in WM. How it relates to etiology and tissue destruction will be an interesting field to explore.

Appendix A
By assuming that backflow during the first-pass of the CA bolus was negligible, Eq. 1 can be used for kinetic analysis of FP MRI data. 7,8 Note that the two terms on the right side of Eq. [1] represent the intravascular (C vas (t)) and interstitial components (C ees (t)), respectively, ie:  first pass of the CA bolus was based on the integral of the input curve, which was termed the "leakage profile" (LP) in previous studies 7, 8 : where T rR is the time of the beginning of the recirculation phase identified from the input function. 7,8 An iterative analysis was proposed for separation of the intra-and extra-components of the enhancement curve 8 :

Appendix B
The 4D LDHT DCE-MRI data were spatially aligned with and resliced to the high spatial resolution series of the ICR-DICE, 11,30 using SPM2. 13 The voxel size of the spatially aligned and resliced HT images is much smaller, ie, 1 3 1 3 2 mm 3 (after reslicing) vs. 2.5 3 2.5 3 6.35 mm 3 (before reslicing). With the high temporal resolution (Dt 5 1s), for catching up the peak of Gd-contrast enhancement, the SSS could be easily outlined in several consecutive slices of the axial orientation image slab, in spite of the low contrast dose. A semiautomatic procedure for C p (t) extraction was used.
Operator interaction is limited to the identification of: 1) The axial slab for searching SSS ROIs. The operator should first identify the time frame that caught the peak of the first pass of CA bolus in the SSS, and then find the axial slab (in the 3D volume identified) over which the SSS spanned across. The two end slices of the predefined axial slab were denoted as SL i-end and SL i1end .
2) The slice where the automatic searching started (SL i ). The starting slice was chosen, whereas the SSS was approximately perpendicular to the image plane, preferred to be near the center slice of the predefined axial slab. A small oval or rectangle ROI of a 3 b was manually drawn on the selected slice (SL i ), just covering the SSS. The ROI associated with the selected slice (SL i ) was denoted as ROI i . The automatic searching for SSS ROIs was performed in the order of SL i-1 , SL i-2 , . . ., SL i-end , and then SL i11 , SL i12 , . . ., SL i1end . These ROIs were predefined as of the same shape and size as ROI i , but might need to be slightly repositioned along the SSS runs. The searching on slice SL i-1 started from the X-Y-coordinates coincided with the center of ROI i , but within a region with extended size (2a 3 2b). The search ended when the ROI i-1 achieved its maximum mean SI values. The search of ROI i-2 on SL i-2 repeated as for ROI i-1 , but in a 2a 3 2b region whose center coinciding with the X-Y-coordinates of the center of ROI i-1 , and so on. As such, candidate voxels were collected along the SSS over the whole imaging volume. The SI time course for each of the collected voxels was extracted from the 4D LDHT.
The area under the SI enhancing curve within 30 seconds of the arrival time (AUC 30 ) was calculated. A mean SI(t) curve was calculated from N (~50) voxels, which had the highest AUC 30 . The averaged mean SI(t) was then converted to the plasma CA concentration time course curve (C p (t)). The number of voxels, N, for the calculation of the final C p (t), can be altered as different data acquisition protocols are used. In principle, the finer of the spatial resolution, the more voxels could be used for the calculation. Bolus arrival time (BAT) was estimated for each tissue uptake curve. The C p (t) measured from the SSS has to be time-shifted to match the BAT for the kinetic analysis.

Appendix C
The theoretic prediction of the steady-state signal, S, from a transverse-spoiled gradient-echo acquisition, assuming short echo time (TE ( T Ã 2 , T Ã 2 is the effective transverse relaxation time), is given by: where a is the flip angle, TR is the repetition time, R1 is the longitudinal relaxation rate, and M 0 is the equilibrium longitudinal magnetization, which depends on the proton density, the profile of the receiver coil sensitivity, the receiver gain setting, and the scaling parameters used for image reconstruction. 36 Maps of M 0 and R1 N were calculated by fitting the VFA signals with Eq. [C.1]. A nonlinear least squares method was used in the calculation of 3D R1 N and M 0 maps. 7,12 Appendix D There are two ways for calculation of R1 post based on R1 N , mean preinjection signal intensity (SI pre ), and postinjection signal intensity (SI post ): By division: By subtraction 16 : (D. 2) The subtraction method needs both R1 N and M 0 for calculation of R1 post . The M 0 value obtained from fitting of the VFA signal can be used for this purpose only when all the factors affecting the value of M 0 are identical in the VFA and DCE sequences. That can be achieved when the DCE acquisition repeats one of the sequences in the VFA acquisition. In this study, the subtraction method was used for R1 post calculation to reduce computation instability that might occur in case of fluctuation in small values of SI pre . 16,37