Sources of systematic error in DCE‐MRI estimation of low‐level blood‐brain barrier leakage

Dynamic contrast‐enhanced (DCE) ‐MRI with Patlak model analysis is increasingly used to quantify low‐level blood‐brain barrier (BBB) leakage in studies of pathophysiology. We aimed to investigate systematic errors due to physiological, experimental, and modeling factors influencing quantification of the permeability‐surface area product PS and blood plasma volume vp, and to propose modifications to reduce the errors so that subtle differences in BBB permeability can be accurately measured.


| INTRODUCTION
Dynamic contrast-enhanced (DCE-) MRI is a frequently used research technique in assessing breakdown of the blood-brain barrier (BBB) in neurological diseases, including dementias, multiple sclerosis, tumors, and stroke. 1 By measuring signal changes following intravenous injection of a gadoliniumbased contrast agent (GBCA), DCE-MRI provides quantitative measurements of GBCA leakage from the brain's blood vessels into the interstitial space as the permeability-surface area product PS, together with other physiological properties including the blood plasma volume fraction v p . DCE-MRI has long been applied to investigate tissues where vascular leakage rates are high, such as high-grade brain tumors. However, the technique is now increasingly used for research into conditions such as cerebral small vessel disease (SVD), 2 Alzheimer disease (AD), 3,4 multiple sclerosis, 5 and aging, 6 where BBB pathology is considered to be significant but the degree of GBCA leakage is nevertheless low.
Despite increasing application of DCE-MRI in this setting, the subtle nature of the leakage and technical limitations of current MRI technology limit the precision and accuracy of BBB leakage estimation. High variability in PS measurements mean that large sample sizes are required to detect biological effects, while large inter-site variability hinders multi-center studies and meta-analyses. Improvements in accuracy and precision, therefore, would enable clinical associations, disease progression, and treatment effects to be established more reliably and efficiently. The significant impacts of noise, signal drift, and pharmacokinetic model selection have been previously investigated as summarized in recent reviews, together with recommendations for image acquisition and analysis. 7,8 For example, there is consensus that pharmacokinetic analysis using the Patlak approach is optimal in this regime; previous simulation studies show that its two key assumptions (high cerebral blood flow relative to the BBB leakage rate and negligible backflux of GBCA) are typically met and that good model fits are obtained with minimum free parameters. [9][10][11] However, previous reviews have also identified other factors with potential to confound quantification (regardless of the model used to predict GBCA distribution) that have not been quantified in this context. For example, it is standard practice in DCE-MRI studies of neurodegenerative diseases to model longitudinal relaxation in the fast water exchange limit (FXL). This assumption, which requires the differences in longitudinal relaxation rates between the tissue compartments to be much smaller than the inter-compartmental water exchange rates, 12,13 might not be valid for subtle GBCA BBB leakage, particularly during and immediately following injection when the vascular longitudinal relaxation rate is high. Therefore, assuming relaxation rates for vascular and extravascular compartments are distinct may provide more accurate estimates of leakage rates. It is also standard practice to use a bolus injection of GBCA, but although fast injection has benefits over constant infusion, 14 the difficulty of resolving rapid first-pass changes in vascular and tissue GBCA concentrations may introduce additional errors. As a consequence of these and other sources of error, and of divergent approaches to study design, image acquisition, and data analysis, reported values of PS (or, equivalently here, the volume transfer constant K Trans ) and findings in relation to pathology vary widely between studies and research sites, 7 limiting the adoption of these parameters as surrogate biomarkers of BBB integrity.
In this work, we quantify the impact of experimental and physiological effects on measurements of subtle BBB leakage using a Monte-Carlo simulation approach. We hypothesized that BBB water exchange, cerebral blood flow, arterial input function (AIF) delay, flip angle (FA) mis-calibration, and B + 1 inhomogeneity would all impact the accuracy of PS and v p estimation in the slow-leakage regime. We further evaluated modifications to the DCE-MRI acquisition and analysis pipelines hypothesized to attenuate these effects, specifically: (i) bolus injection versus slow GBCA injection, (ii) inclusion versus exclusion of first-pass tissue data from the model fitting, (iii) fitting data under the FXL versus no exchange limit (NXL; equivalent to the slow exchange limit SXL) assumption reduced (≤0.56 × 10 −4 min −1 ) by applying modifications to the acquisition and processing pipeline. Processing modifications also had substantial effects on in-vivo normal-appearing white matter PS estimation (absolute change ≤ 0.45 × 10 −4 min −1 ). Conclusion: Measuring subtle BBB leakage with DCE-MRI presents unique challenges and is affected by several confounds that should be considered when acquiring or interpreting such data. The evaluated modifications should improve accuracy in studies of neurodegenerative diseases involving subtle BBB breakdown.

K E Y W O R D S
blood-brain barrier, DCE-MRI, dementia, gadolinium, small vessel disease for exchange across the BBB, and (iv) B + 1 correction. The impacts of processing modifications were also evaluated in-vivo in a cohort of mild stroke patients with varying degrees of SVD severity, with PS and v p estimated in normal-appearing white matter (NAWM), subcortical gray matter (scGM), and white matter hyperintensities (WMH). We provide an opensource MATLAB application with graphical user interface so that measurement accuracy and precision can be estimated for arbitrary DCE-MRI protocols and modeled physiology. The implications of our findings for planning future studies and for interpreting results are discussed.

| Simulations
DCE-MRI time series were simulated using publically accessible scripts (https://doi.org/10.7488/ds/2997) written inhouse using Matlab (Mathworks, Natick MA, USA). Except where specified, the white matter tissue physiological and MR parameters used to simulate data are provided in Table  1. Spoiled gradient echo (SPGR) time series were simulated using the in-vivo protocol described in Section 2.2.1 (repetition time/echo time [TR/TE] = 3.4 ms/1.7 ms, FA = 15°, 3 pre-and 29 post-injection acquisitions with total duration 1268 s). The simulated protocol follows recent consensus recommendations for measurement of subtle BBB leakage. 8 2.1.1 | Pre-contrast T 1 measurement T 1 measurement using the variable FA (VFA) method (TR = 5.4 ms, FA = 2°, 5°, 12°) was simulated by generating synthetic signals including additive Gaussian noise, with variance adjusted to achieve a test-retest T 1 error comparable to that reported in the literature for 3T MRI (0.7%). 15 The simulated signal was fitted using a non-linear least squares minimization approach to obtain the measured T 1 .

| Simulated DCE-MRI signaltime courses
DCE-MRI simulations were performed using two synthetic AIFs corresponding to bolus injection and slow injection over 2 min (Figure 1). A bolus AIF was generated using a population average function following a similar approach to previous studies, 9,10 while a population average slowinjection AIF was constructed using patient measurements, as described in Section 2.2. Full details are provided in the Supporting Information. A venous vascular input function (VIF), which is frequently used for modeling in-vivo data, was obtained by time-shifting the AIF by 6 s. 16,17 First, the simulated high-temporal-resolution AIF was convolved with the impulse response function for the 2-compartment exchange Model (2CXM) to generate tissue, capillary plasma, and EES concentration curves as described previously. 5,9,10 This model is defined by four parameters: PS, EES volume fraction v e , blood plasma volume fraction v p , and blood plasma flow F p . MRI signal was calculated for steady-state SPGR acquisition; BBB water exchange effects were modeled using the 2-site-1-exchange model (2S1X). 18,19 Second, signal-time courses were sampled at the required temporal resolution (Δt = 39.6 s) and random Gaussian noise added to achieve a temporal signal-to-fluctuation-noiseratio (SFNR) of 230 to match the average measured in-vivo value. SFNR was calculated as the ratio of the baseline precontrast signal to the root mean square of the residuals for the VIF and tissue concentrations were estimated from the simulated MRI signal and T 10 , assuming the FXL as per standard practice in this application. 10 GBCA concentration was estimated analytically using the equation for the SPGR signal 21 and assuming linear dependence of 1/T 1 on the concentration (r 1 = 5.0 mM −1 s −1 ). 22 Concentration in blood plasma c p was determined from the blood concentration c b as c p = c b ∕ (1 − Hct), where Hct is the hematocrit. Finally, the Patlak model 23 was fit to tissue concentration-time curves using an unconstrained multiple linear regression approach 8 to yield estimated v p and PS. PS and v p were also estimated assuming the NXL for water exchange across the BBB. A forward model was defined to predict signal enhancement for a given PS and v p . Pre-contrast compartmental T 10 was calculated assuming the FXL as described previously, 8 since the rate of water exchange across the BBB is likely to be much greater than the corresponding difference in compartmental R 1 . 12,13 GBCA concentrations in the capillary and EES compartments were predicted using the Patlak model, and the corresponding MRI signals were calculated assuming two distinct well-mixed compartments, that is, vascular and extravascular (the combined EES and intracellular compartments) spaces. The two compartmental signals were combined, weighted by their volume fractions, to predict an overall tissue signal. PS and v p were adjusted to minimize the sum-of-square differences between predicted and "measured" enhancements using the Matlab lsqcurvefit function.
All simulations were repeated 1000 times to assess the influence of noise.

| In-silico experiments
Using the simulation framework described in the preceding sections, we assessed the accuracy of PS and v p measurements in white matter and their sensitivity to variations in physiological and experimental parameters. First, we explored the impact of the BBB water exchange rate, using k be = 2.75 s −1 as a typical value for healthy NAWM, as well as half of and double this value, to cover potential variation due to age and subtle BBB pathology. 18 Second, we determined the impact of cerebral blood plasma flow rate F p , using the value 11 mL 100 g −1 min −1 (equivalent to CBF = 19 mL 100 g −1 min −1 ) representing NAWM in subjects with a similar age range to that of our clinical cohort, 24 and lower values (8.25 and 5.5 mL 100 g −1 min −1 ) to simulate conditions of chronic ischaemia. Third, the effect of AIF delay due to unknown contrast arrival time in the brain circulation (for example due to variations in injection timing, path length and cardiac output) was simulated by repeating simulations with the AIF time-shifted by 0-12 s relative to the reference AIF. For each of these three effects, simulations were performed assuming bolus-and slowinjection protocols, and with and without exclusion of the first 3 data points from the start of the injection (corresponding to a period of 2 min) from the model fitting (specifically, data points were excluded for the purpose of calculating the sum of squares difference between the data and Patlak model prediction, but were always included when calculating the model tissue concentration, which requires integration of the VIF).
Finally, we explored the impact of B + 1 variation and uncertainty by introducing a proportional error K FA in the transmitted flip angle FA true for both the VFA T 1 and DCE-MRI acquisitions, such that FA true = FA nom × K FA , where FA nom is the nominal FA. Simulations were performed assuming (i) zero FA error, (ii) equal FA error for tissue and the VIF and (iii) unequal FA error for tissue and the VIF. K FA values were obtained from measurements in our clinical cohort using the DESPOT1-HIFI technique, as outlined in section 2.2; values were first averaged over voxels within the VIF and NAWM masks and second over all patients. Simulated data were analyzed assuming the nominal and actual FAs to determine the effect of B + 1 correction.

| In-vivo experiments
We used data from the first 50 patients with recent nondisabling ischaemic stroke recruited into an ongoing prospective longitudinal study of cerebral SVDs, as per published protocol. 25 MRI was performed ≥1 month post-stroke to minimize acute effects of stroke on regional BBB integrity. The study was approved by the South-East Scotland Research Ethics Committee (REC 18/SS/0044) and informed written consent was obtained for each participant.

| MRI
Structural imaging, T 1 measurement and DCE-MRI were acquired using a MAGNETOM Prisma 3T clinical MRI scanner (Siemens Healthcare GmbH, Erlangen, Germany) with a 32-channel receive head coil (acquisition parameters are provided as supplementary material of Ref. 25).
DCE-MRI was acquired using a 3D sagittal T 1 w SPGR sequence with non-selective RF excitation and whole-brain coverage (TR/TE = 3.4/1.7 ms, FA = 15°, acquisition matrix size 120 × 96 × 96, 2-mm isotropic resolution). A total of 32 volumes were acquired with a temporal resolution of 39.6 s over 21 min. A dose of 0.1 mmol/kg body weight gadobutrol (1 M Gadovist, Bayer AG, Leverkusen, Germany) was injected intravenously using a power injector following acquisition of three pre-contrast volumes, followed by a 20 mL saline flush; the flow rate was adjusted to deliver the required dose volume over a period of 110-130 s.
Tissue masks representing NAWM, scGM, cerebrospinal fluid, WMHs (a main indicator of SVD), and stroke lesions were derived from structural images, as described in the Supporting Information.

| DCE-MRI processing
DCE-MRI images were spatially realigned (SPM 12 https:// www.fil.ion.ucl.ac.uk/spm/) using the brain-extracted (FSL-BET 28 ) volumes. A mean pre-contrast image was obtained by averaging the first three volumes and was used as the target image for transforming and resampling the T 10 maps and binary tissue masks into the DCE-MRI space using FSL FLIRT. 29 As in previous studies of subtle BBB leakage and consistent with recommendations, 8,30 a venous VIF was used, without scaling or transformation, in order to minimize partial volume, head motion and inflow effects that can affect AIF measurement. The VIF was obtained in each patient by manually selecting five voxels from the superior sagittal sinus, preferentially selecting voxels with a high post-injection signal peak and smooth post-contrast signal decay, and obtaining the mean signal-time course of these voxels. This approach was previously found to have good interobserver reproducibility. 31 The median signal at each time point was obtained from voxels within each tissue mask, 10 in order to suppress the potential impact of skewed intensity distributions and outlier voxels. The signal enhancement was then calculated at each time point relative to the average pre-contrast signal, and concentration was estimated analytically based on the SPGR signal equation and assuming linear dependence of 1/T 1 on the concentration. PS and v p were estimated as described in Section 2.1.3 using patient Hct measurements; where no Hct was available from the time of the scan, previous values were obtained from patient records.

| Statistics
The error in simulated PS and v p measurements was defined as the difference between the fitted and ground-truth values. The mean and SD of the errors over all runs of the simulation were recorded to indicate the predicted systematic bias and random error, respectively. We defined the sensitivity s of estimated PS or v p to variation in physiological parameters (ie, k be , F p and AIF delay) for a given acquisition and processing approach as the range of PS or v p estimates across all tested values of the confounding parameter, averaged over all simulation runs and simulated values of PS or v p . s was calculated excluding the notional water exchange cases k be = 0 and k be = 1000 s −1 . In-vivo measurements were summarized using descriptive statistics and differences were assessed using the paired samples t-test.

| Water exchange
To determine the impact of BBB water exchange on BBB leakage measurement, DCE-MRI data were simulated for white matter using the 2CXM and 2S1X models to simulate | 1893 MANNING et Al.

| Cerebral blood flow
To determine the additional impact of cerebral blood flow on BBB leakage measurement, DCE-MRI data were simulated F I G U R E 2 Simulation results showing the effect of variable BBB water exchange rate k be on PS estimated using the standard Patlak analysis assuming fast water exchange (two left-most columns) and by fitting the same model under the assumption of no BBB water exchange (two rightmost columns). Simulations are shown for both bolus and slow injection acquisitions, and with or without exclusion of early data points. Note that a wider y-axis range is used in A. Results are shown for typical (2.75 s −1 ), low (/2), and high (×2) values of k be in NAWM, as well as notional values representing the no-and fast-exchange limits (k be = 0, 1000 s −1 ). Error bars show the mean ± SD errors for 1000 simulations; thus, the lines indicate systematic error (bias), while the error bars indicate random error due to noise. Source code to generate this figure is available to download at https://doi.org/10.7488/ds/2997 as described above, assuming different blood plasma flow rates F p but constant water exchange effects (k be = 2.75 s −1 ) (Figure 3). PS measurements are sensitive to F p but F p sensitivity is greater for bolus versus slow injection (0.39 vs. 0.19 × 10 −4 min −1 for NXL estimation). For both injection protocols, the sensitivity of the fitted parameters to F p is virtually eliminated by excluding the first three post-injection data points from the cost function (F p sensitivity ≤0.03 × 10 −4 min −1 for NXL estimation). Findings were qualitatively similar using the standard FXL fitting approach and for v p estimation ( Figure S2).

| AIF delay
To determine the additional impact of AIF delay due to variable contrast arrival time in the cerebral arteries, further simulations were performed in the same manner, with a variable time delay (0, 4, 8, and 12 s) applied to the AIF (Figure 4). PS measurements are highly sensitive to AIF delay for a bolus injection (delay sensitivity 1.00 × 10 −4 min −1 for NXL estimation), but this is reduced when early data points are excluded (0.24 × 10 −4 min −1 ). For slow injection PS sensitivity is negligible, regardless of whether F I G U R E 3 Simulation results showing the effects of variable blood plasma flow rate F p on estimated PS incorporating water exchange effects (k be = 2.75 s −1 ). PS was estimated using the standard Patlak approach assuming fast water exchange (two left-most columns) and by fitting the same model under the assumption of no BBB water exchange (two right-most columns). Simulations are shown for both bolus and slow injection acquisitions, and with and without exclusion of early data points. Note that a wider y-axis range is used in A. Results are shown for typical (11 mL 100 g −1 min −1 ), low (−25%), and very low (−50%) values of F p in NAWM. Error bars show the mean ± SD estimates for 1000 simulations; thus, the lines indicate systematic error (bias) while the error bars indicate random error due to noise. Source code to generate this figure is available to download at https://doi.org/10.7488/ds/2997 early data points are excluded. Findings were qualitatively similar using the standard FXL fitting approach and for v p estimation ( Figure S3).

inhomogeneity
Further simulations were performed to determine the impact of FA error K FA , assuming typical water exchange effects (k be = 2.75 s −1 ) and with exclusion of early data points from the model fitting ( Figure 5). The impact of spatially uniform FA error is very slight: PS estimates are very similar to those obtained with accurate FAs, with low systematic error in PS for bolus (−0.47 to −0.11 × 10 −4 min −1 ) and slow (−0.30 to −0.17 × 10 −4 min −1 ) injection acquisitions with NXL fitting. For the more realistic case where B + 1 inhomogeneity results in different K FA in blood and tissue, the systematic error is substantially greater for both bolus (−1.87 to −0.19 × 10 −4 min −1 ) and slow (−1.90 to −0.38 × 10 −4 min −1 ) injections with NXL fitting. B + 1 correction results in PS estimates close, but not identical, to those obtained in the absence of B + 1 error. Findings were qualitatively similar using the standard FXL fitting approach and for v p estimation ( Figure S4).

F I G U R E 4
Simulation results showing the effects of AIF delay on estimated PS incorporating water exchange effects (k be = 2.75 s −1 ). PS was estimated using the standard Patlak approach assuming fast water exchange (two left-most columns), and by fitting the same model under the assumption of no BBB water exchange (two right-most columns). Simulations are shown for both bolus and slow injection acquisitions, and with and without exclusion of early data points. Note that a wider y-axis range is used in A. Results are shown for AIF delays of 0, 4, 8, and 12 s. Error bars show the mean ± SD estimates for 1000 simulations; thus, the lines indicate systematic error (bias), while the error bars indicate random error due to noise. Source code to generate this figure is available to download at https://doi.org/10.7488/ds/2997

| In-vivo measurements
After excluding one dataset due to severe motion artefact (N = 1), 49 patients (18 female, age 66.4 ± 9.6 y) had useable DCE-MRI, acquired using the slow-injection protocol described in Section 2.2.1.
We quantified the impact of processing modifications on our in-vivo results (Table 2, Figure 6). There is a substantial negative change in mean PS (−0.447 × 10 −4 min −1 in NAWM) and increase in v p (+0.17%) estimates when a NXL fitting approach is used in place of the standard FXL method. Excluding early data points from the fitting had less impact, resulting in slightly lower PS values (−0.04 × 10 −4 min −1 ) and unchanged v p . B + 1 correction had a substantial impact, increasing both PS (+0.201 × 10 −4 min −1 ) and v p (+0.12%) estimates.
Using the processing approach predicted to yield greatest accuracy (NXL fitting, excluding early data points and B + 1 correction), PS values were greater than zero in all tissues (P < .0077). PS was higher in WMH (P < .0001) and scGM (P < .0001) vs. NAWM, and v p was also higher in WMH (P < .0001) and scGM (P < 10 −5 ) vs. NAWM. There was no significant difference in PS between WMH and scGM (P = .8), although v p was higher in scGM vs. WMH (P < 10 −5 ).

F I G U R E 5 Simulation results
showing the effects of FA error on estimates of PS, incorporating water exchange effects (k be = 2.75 s −1 ). PS was estimated using the standard Patlak approach assuming fast water exchange (two left-most columns), and by fitting the same model under the assumption of no BBB water exchange (two right-most columns); for clarity, results are only shown for analysis with exclusion of early data points. Simulations are shown for both bolus and slow injection acquisitions, and with and without B + 1 correction. Note each row uses a different y-axis range. Error bars show the mean ± SD estimates for 1000 simulations; thus, the lines indicate systematic error (bias), while the error bars indicate random error due to noise. Source code to generate this figure is available to download at https://doi.org/10.7488/ds/2997 | 1897 MANNING et Al.

| DISCUSSION AND CONCLUSIONS
In this work, we quantified the impact of biological and instrumental factors on the accuracy of low-level BBB leakage measurements. Simulations predicted clinically significant systematic errors due to injection protocol, restricted water exchange, cerebral blood flow, AIF delay, and B + 1 inhomogeneity. In some cases, the errors have comparable magnitude to the parameters being measured. Although the recommended T A B L E 2 Mean (SD) of regional PS and v p measured in our clinical cohort (N = 49) and estimated using (i) the standard Patlak approach under the assumption of fast water exchange, (ii) with the assumption of no BBB water exchange (NXL fit), (iii) NXL fit with exclusion of early data points from the cost function, and (iv) NXL fit with exclusion of early data points and B + 1 correction of T 1 measurement and DCE-MRI Patlak analysis method was considered in this work, these confounds are likely to affect subtle BBB leakage measurement using any standard model of GBCA distribution with a comparable MRI protocol. Modifications to the standard acquisition and processing methods were shown to reduce the errors. In-vivo results in minor stroke patients with varying degree of SVD severity reflected previously reported tissue differences; however, parameter estimates changed significantly when the evaluated processing modifications were applied.

| Causes and reduction of systematic errors
First, we showed that the assumption of fast BBB water exchange, ubiquitous in DCE-MRI studies of neurodegenerative diseases, 8 results in substantial bias in PS (up to 4.48 × 10 −4 min −1 ). This is caused by violation of the FXL assumption, with consequent underestimation of GBCA concentration ( Figure S5). Since this effect is strongest during the early part of the time course, GBCA concentration appears to increase more rapidly (or decrease more slowly), leading to overestimation of PS. The importance of water exchange in contrast-enhanced MRI has long been known, 32 but its impact has not been widely recognised in the BBB literature. Cao et al. reported the effect of restricted water exchange on the GBCA transfer constant, although the impact on accuracy was not presented. 33 Our findings are qualitatively similar to those of Paudyal et al, 12 who simulated measurements of faster BBB leakage with the extended Tofts model. Larsson et al. reported a significant impact on cerebral perfusion values measured using the inversion-recovery SPGR sequence. 34 We showed that predicted systematic errors are reduced for a bolus injection by excluding early post-injection data points from the model fit, where violation of the FXL is most pronounced. Bias is also substantially reduced by modeling the MRI signal in the NXL. However, PS estimates remain somewhat sensitive to the value of k be , which is likely to be altered in neurodegenerative pathology. 18 Variation in individual AIF shape is also likely to affect the bias and requires further investigation. Other approaches could include: acquisition protocols that facilitate simultaneous estimation of k be and PS; independent measurement of k be for use in the DCE-MRI signal model; reducing the GBCA dose (at the likely cost of reduced sensitivity to BBB leakage) 35 ; and optimizing the acquisition protocol to reduce k be sensitivity. 32,36 Second, we simulated the impact of cerebral blood flow, which is commonly reduced in aging and neurovascular disease. 37 Previous work has shown that excluding early postinjection time points improves the Patlak model fit, reducing blood flow effects. 10,11 Our simulations results support this approach and, in the case of bolus injection, predict that large systematic PS errors and F p -dependence would otherwise result. Slow-injection DCE-MRI is less sensitive to F p , since the slower change in arterial GBCA concentration leads to similar arterial, capillary, and venous concentrations ( Figure  S6), but exclusion of early data points nevertheless reduces F p sensitivity.
Third, we probed the effect of AIF delay, which can vary in-vivo due to multiple factors including injection timing, cardiac output and path length from the injection site. For a bolus injection, a small delay significantly affects the estimated parameters and a significant error remains even when early data points are excluded, while slow injection of GBCA virtually eliminates the impact. The reason, illustrated in Figure S7, is that the area under the VIF, which approximately determines GBCA extravasation, is more accurately quantified for a slowly changing VIF. For the same reason, additional errors due to variation in bolus shape will likely also be reduced for a slow injection. Furthermore, the high peak concentrations induced by a bolus injection, which may cause T 1 saturation and transverse dephasing, are avoided. This problem of VIF sampling can alternatively be addressed by acquiring data at higher temporal resolution-either during the first-pass 38 or throughout the acquisition-in order to resolve the first-pass profile. However, a compromise must be struck between the temporal resolution, spatial resolution, anatomical coverage and CNR requirements of the application. In this work, we employed a protocol with whole-brain coverage and moderate spatial and temporal resolutions in order to capture pathological changes across multiple brain regions; whole-brain non-selective excitation has the additional benefit of reducing inflow artefact in the large vessels, which can affect VIF quantification. Spatial resolution should nevertheless be adequate to resolve the VIF vessel and the anatomy of interest. A third solution is to substitute a population-averaged VIF during the first-pass period; however, this could lead to errors due to individual variation.
Fourth, we investigated the effect of B + 1 error. Our results confirm that B + 1 effects largely self-cancel for the SPGR sequence, as previously reported, 39 provided that the FA error is the same for tissue and VIF. However, our in-vivo measurements using a 3T MRI scanner revealed a 17% difference in average FA for NAWM and VIF. In this case, substantial additional errors in PS (around −40%) and v p were predicted. The reason is illustrated in Figure S8. Fortunately, B + 1 correction was shown to virtually eliminate this measurement error.

| In-vivo findings
In our patient cohort, NXL fitting yielded PS estimates that were substantially closer to zero in all tissues compared with standard FXL fitting (0.60 vs. 0.16 × 10 −4 min −1 in NAWM), while v p estimates increased by approximately 50%. The incremental effect of excluding early data points was small, as predicted by simulations for the case of a slow injection. B + 1 correction resulted in higher PS and v p , as predicted; however, final estimates remained lower for PS (0.32 vs. 0.60 × 10 −4 min −1 in NAWM) and higher for v p (0.66 vs. 0.38%) compared with standard FXL analysis without exclusion of early data points or B + 1 correction. The wider PS distribution following B + 1 correction may represent biological variation, since predicted B + 1 errors ( Figure 5) become more negative as PS increases and, therefore, would compress the apparent distribution without B + 1 correction. In-vivo data processed with the NXL assumption, exclusion of early data points and B + 1 correction suggested detectable BBB leakage in all three tissues. Leakage was greater in WMH and scGM compared with NAWM, but similar in WMH and scGM. This pattern of tissue differences confirms the findings of our earlier 1.5T DCE-MRI study in a similar patient cohort. 10 However, PS estimates are lower in the present 3T study (0.3 vs. 3.0 × 10 −4 min −1 in NAWM). There are several possible reasons for this. Sham 1.5T DCE-MRI scans without contrast revealed a positive signal drift of approximately 0.08 % min −1 , resulting in predicted overestimation of PS by around 3.0 × 10 −4 min −1 , while 3T signal drift was close to zero. In addition, 1.5T data were processed using the standard FXL approach, which would likely have resulted in further overestimation. The impact of B + 1 inhomogeneity, which was not measured at 1.5T, use of different GBCAs and patient differences could also have contributed.
The greater BBB leakage rate measured in WMH versus NAWM in both studies is consistent with BBB dysfunction as a component of WMH pathophysiology, as seen histologically. 40 However, findings from other SVD studies are varied. 8 PS is determined by capillary surface area as well as permeability 41 ; thus, increased blood volume and potential differences in capillary width 42 could contribute. v p was also higher in WMH compared with NAWM, and higher still in scGM. This also confirms our 1.5T findings, however absolute values are somewhat higher in the present study (0.7 vs. 0.6% in NAWM, 1.5 vs. 1.2% in scGM), presumably due to the factors discussed above. Both sets of values are lower than reference values obtained in individuals of similar age using 15 O positron emission tomography. 24 This may reflect differences in the cohorts (patients with cerebrovascular disease versus healthy volunteers), resolution and tissue sampling (eroded NAWM masks vs. regions of interest).

| Implications
Our findings have significant implications for studies of subtle BBB permeability. Inter-site differences in reported PS measurements have comparable magnitude 7 to the predicted systematic errors and might be partly caused by the experimental and biological aspects considered in this work, in addition to the previously described impact of scanner drift. Reported differences between and within patients also have a similar magnitude to the predicted systematic errors. For example, Zhang et al. compared SVD patients with controls, reporting a PS difference of approximately 0.5 × 10 −4 min −1 43 ; Montagne et al. reported a difference of around 5 × 10 −4 min −1 according to apolipoprotein E status 3 in cognitively normal or early AD participants; finally, Heye et al. reported a difference of 10 −4 min −1 between NAWM and WMH in patients with recent mild stroke. 10 Since variation in BBB water exchange rate, cerebral blood flow, and other factors could potentially mimic or obscure differences in the actual PS, their impact should be estimated for specific MRI protocols and, where possible, reduced when designing studies, for example, using the source code provided by the authors.
Some general inferences can also be made. First, it is essential to exclude first-pass data from the model cost function in the case of a bolus injection, in order to avoid large errors due to BBB water exchange, cerebral blood flow rate, and AIF delay. This is not a new recommendation 11 but implementation or otherwise is often not reported in the literature 8 and our simulations reveal the high impact of omitting this measure. Exclusion is less essential in the case of slow injection but nevertheless reduces sensitivity to blood flow.
Second, the substantial systematic errors resulting from the standard FXL assumption can be reduced by modeling the MRI signal in the NXL, with the caveat that PS estimates remain somewhat sensitive to the AIF shape and to the water exchange rate.
Third, parameters measured using a bolus injection at low temporal resolution remain moderately sensitive to patient-and protocol-dependent variation in AIF timing, for example due to cardiac output or injection start time. These errors can be eliminated via a slow injection approach. Similar errors are likely to result from variation in bolus shape and arteriovenous delay, which should also be reduced with a slow injection. A disadvantage of this approach is that the full GBCA dose is delivered later; thus, a faster injection rate could be more sensitive to GBCA leakage for shorter acquisitions than simulated here. Further research is needed before recommending an approach, in order to determine the inter-patient variation in VIF profiles, which also modulate the impact of water exchange.
Fourth, FA mapping substantially improves the accuracy of subtle BBB leakage. Failure to correct for variation in the transmitted FA results in substantial errors that vary spatially, could confound tissue comparisons, and that may be larger than those predicted here for the relatively homogeneous B + 1 field of our dual-transmit 3T MRI system. B + 1 shimming, where available, should also improve accuracy.

| Strengths and limitations
Strengths of our work include the use of a computational approach to determine the impact of factors that cannot be easily measured and controlled in-vivo. We propose remedies that were evaluated theoretically and demonstrated in-vivo. Our source code is made publically available for researchers to determine the impact of these effects in future studies. A limitation of our work is that simulation results were obtained for specific sets of parameters, and would be quantitatively different for other protocols and tissues; however, the protocol evaluated is representative of those reported in the literature and the supplied software permits alternative settings to be explored. All simulations are also limited by the forward model used to generate the signal, since knowledge of all underlying biophysical properties is limited. For clarity, we did not include signal drift in the simulations, since this effect has been addressed previously 9,10 and is scanner-dependent. AIF-VIF delay, but not dispersion, was simulated, since there is limited literature regarding the latter and previous work indicates similar time courses 16,17,44 ; dispersion will primarily affect the VIF during the period immediately after injection and is likely to affect quantification in a similar manner to variation in AIF timing (Figure 4). The impact should, therefore, be minimized by excluding early data points or using a slow contrast injection.

| CONCLUSIONS
In conclusion, we investigated significant confounds of DCE-MRI subtle BBB leakage measurement and approaches to reduce their impact. Properties such as water exchange and blood flow rates, if correlated with pathology, could mimic or disguise differences in BBB leakage if not addressed. Invivo results in patients with imaging markers of SVD were consistent with simulations and previous findings. Our work should permit improved measurement of BBB leakage in future studies.