Reliability of dynamic contrast-enhanced magnetic resonance imaging data in primary brain tumours: a comparison of Tofts and shutter speed models

Purpose The purpose of this study is to investigate the robustness of pharmacokinetic modelling of DCE-MRI brain tumour data and to ascertain reliable perfusion parameters through a model selection process and a stability test. Methods DCE-MRI data of 14 patients with primary brain tumours were analysed using the Tofts model (TM), the extended Tofts model (ETM), the shutter speed model (SSM) and the extended shutter speed model (ESSM). A no-effect model (NEM) was implemented to assess overfitting of data by the other models. For each lesion, the Akaike Information Criteria (AIC) was used to build a 3D model selection map. The variability of each pharmacokinetic parameter extracted from this map was assessed with a noise propagation procedure, resulting in voxel-wise distributions of the coefficient of variation (CV). Results The model selection map over all patients showed NEM had the best fit in 35.5% of voxels, followed by ETM (32%), TM (28.2%), SSM (4.3%) and ESSM (< 0.1%). In analysing the reliability of Ktrans, when considering regions with a CV < 20%, ≈ 25% of voxels were found to be stable across all patients. The remaining 75% of voxels were considered unreliable. Conclusions The majority of studies quantifying DCE-MRI data in brain tumours only consider a single model and whole tumour statistics for the output parameters. Appropriate model selection, considering tissue biology and its effects on blood brain barrier permeability and exchange conditions, together with an analysis on the reliability and stability of the calculated parameters, is critical in processing robust brain tumour DCE-MRI data.


Introduction
Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is a non-invasive methodology that allows tissue perfusion and permeability to be quantified through analysis of T 1 -weighted MR images acquired before, during and after an intravenous (IV) injection of a gadolinium-based contrast agent (CA).
Nowadays, DCE-MRI finds application in a wide range of oncological studies. In brain tumours, the Tofts model (TM), also known as the standard model, together with its extended version (EMT), are regularly applied [1][2][3][4][5][6]. There are a few examples of the application of shutter speed model (SSM) in brain tumours in the literature [7,8] while it has been mainly implemented in the study of breast cancer [9], prostate cancer [10] and hepatocellular carcinoma [11]. The extended shutter speed model (ESSM), also called second generation shutter speed model, or BALDERO (blood agent level dependent and extravasation relaxation overview), has previously found applications only in hepatocellular carcinoma [12] and simulated data [13]. These four models assume that the CA passes readily between the intravascular compartment and the tissue interstitium. However, this assumption is not valid in the presence of the intact blood-brain barrier (BBB), where there is negligible leakage. In this case the CA will only affect the intravascular T 1 value on first pass of the bolus. In addition, the vascular component of a brain tumour is heterogeneous and certain regions may not be perfused (e.g. necrotic regions). We therefore also implemented a no-effect model (NEM) which assumes negligible effect of the CA on the voxel T 1 properties.
In the last decade, the number of pharmacokinetic studies of DCE-MRI brain data has increased as this quantification technique has been found to be indicative of the malignant grade of brain tumours [14]. However, two main issues need to be taken into account when analysing DCE-MRI data. Firstly, pharmacokinetic studies often show different and discordant results, thus bringing the reliability of this quantification technique into question [15,16]. In fact, DCE-MRI data analysis is affected by (a) the acquisition protocol (trade-off between spatial and temporal resolution) [17] and (b) the quantification procedure. Furthermore, most studies present only results from application of a single pharmacokinetic model and the consequent statistical analysis of averaged values evaluated over the whole tumour volume [18][19][20]. This approach completely ignores the particularly heterogeneous nature of brain tumour vasculature and vascular permeability.
In this prospective study, we investigate the robustness of pharmacokinetic modelling of DCE-MRI brain data. We propose a method to identify reliable DCE-MRI data based on a model selection procedure, building on previous work by Bagher-Ebadian et al. [3], and a stability test: five different pharmacokinetic models (TM, ETM, SSM, ESSM and NEM) are assessed and the Akaike information criteria index is used for model selection; we then evaluate the stability of each parameter extracted from the model of choice, in terms of coefficients of variations (CVs), through a noise propagation procedure.

Patient population
Fourteen patients (7 male, 7 female; aged 23-73 years, mean 40 years) with primary brain tumours were recruited to this study. Ethical approval was given by the local ethics committee and informed consent was obtained from all patients. Patients had an MRI at diagnosis, prior to receiving any treatment. Following surgery, histopathological data showed three patients with WHO grade IV glioblastoma, two patients with WHO grade III astrocytoma, two patients with WHO grade III oligodendroglioma, three patients with WHO grade II astrocytoma, three patients with WHO grade II oligodendroglioma and one patient with a WHO grade I dysembryoplastic neuroepithelial tumour.

DCE-MRI data acquisition
MR images were acquired on a 3-T Siemens Verio MRI system using a 32 channel head coil; including pre-and post-contrast T 1 -weighted images, T 2 -FLAIR images and a DCE sequence with a variable flip angle pre-contrast T 1 map acquisition. The DCE-MRI protocol included five pre-contrast spoiled gradient recalled echo (SPGRE) 3D vibe sequences at five different flip angles (2,8,12,15,20,26), and a dynamic 3D vibe sequence (TR = 3.34 ms, TE = 0.99 ms, flip angle = 26, FOV 240 × 240 mm, acquisition matrix 128 × 128, slice thickness 5 mm, slice gap 1 mm, 80 volumes). To obtain acquisitions before, during and after the injection of the CA, 0.1 mmol/kg body-weight gadolinium-diethylene triaminepentacetate (Gd-DTPA, Gadovist) was injected using a power injector on the fifth acquisition using a flow rate of 3 ml/s immediately followed by 20 ml saline solution. The 3D acquisition allowed us to cover the entire brain; 80 time points were acquired with an average temporal resolution of 2.89 s and a total acquisition time of ≈ 4 min.

DCE-MRI data analysis
Volumes of interest (VOIs) were drawn by a Radiologist and confirmed by a Consultant Neuroradiologist for each patient around T 2 -FLAIR hyperintense regions and on a 2-cmdiameter circular region in normal-appearing contra-lateral white matter. DCE-MRI data were analysed using a semiautomated in-house software written in MATLAB (Mathworks, R2017a). Before the application of any of the aforementioned pharmacokinetic models, we calculated the relaxation rate at baseline (R 10 ) and relaxed signal (M 0 ) as 3D maps, with the Ernst formula (assuming TE < <T 2 *) using the set of SPGRE pre-contrast images acquired at different flip angles [21]: where α is the flip angle having values [ (2,8,12,15,20,26]) and TR is 3.34 ms. Reformulating Eq. 1 as a linear regression system (y = cx + d) following the method described by Liberman et al. in [22], gives: where E ¼ e −R 10 TR .
The slope c = E and intercept d = M 0 (1 − E) can thus be estimated and continuing from [22], R 10 and M 0 can then be obtained through: Then, 4D (x, y, z, t) post-injection longitudinal relaxation rate R 1 (t) maps for each dynamic phase are calculated using signal intensity data from the post-contrast dynamic series: where α = 26, A = S(t) − S(0)/M 0 sin(α), B = (1 − E)/1 − E · cos(α). S(0) and S(t) are the pre-contrast injection signal intensity and the signal at the dynamic phase t, respectively [23]. The longitudinal relaxation rate is determined in order to calculate the concentration of the CA. This is done through a calibration between the concentration of CA [CA] and the measured H 2 O MR signal. This can be modelled by either a linear or nonlinear relationship as described below.

Tofts model (TM)
In applying the Tofts model to brain tumours, we expressed the flux of the tracer across two well-mixed compartments (blood and the extravascular, extracellular space) through the volume transfer constant K trans [24]. The TM, also called the standard model, assumes negligible plasma compartment and can describe weakly vascularized brain tissue. More importantly, it assumes a linear dependence of R 1 on [CA] (that is the equivalent of assuming the equilibrium transcytolemmal water exchange kinetics in the fast exchange limit(FXL)): where r 1 is the CA relaxivity. The extravasation of the contrast from the plasma to the extravascular extracellular space (EES) was accounted by the Kety-Schmidt rate law [25]: where K trans is the first-order rate constant for plasma to interstitium CA transport (min −1 ) and v e is a measure of the EES volume fraction. The ratio between K trans and v e results in the third pharmacokinetic parameter k ep , which is the back flux rate constant (min −1 ). [CA 0 ] and [CA p ] are the concentration of CA in the 'outside' space (the extravascular extracellular space) and in the plasma, respectively. [CA p ] is also called the arterial input function (AIF). The fitting of Eq. 6 was performed using the inbuilt MATLAB fminsearch function, which uses the Nelder-Mead simplex algorithm as described in Lagarias et al. [26]. The minimization procedure is done voxel-wise in order to obtain a 3D map for each pharmacokinetic parameter. We set input values of 0.1 and 0.01 for the K trans and v e , respectively, and run the algorithm with 10,000 iterations and a tolerance of 10 −8 . The fitting procedure was also carried out with the user developed MATLAB function fminsearchbnd [27]. This function takes into consideration boundaries in the output values settings, which were set so as to consider only positive values. A comparison between the two functions' results was done in terms of goodness of fit by estimating the Akaike Information Criteria (AIC) for each method using Eq. 7: where n is the number of data points, k the number of fitted parameters and RSS is the residual sum of squares [28].
In general, when performing model selection using the AIC, the model resulting with the lowest AIC value is the model that represents the best balance between complexity (i.e. the number of parameters) and goodness of fit (i.e. lower RSS). In this case, as the number of parameters is the same, only the goodness of fit is being tested. A comparison between the AIC maps relative to the bounded and unbounded procedure allowed the estimation of final K trans , k ep and v e maps where the value of each voxel was obtained from the fitting procedure with the best fit (lowest AIC) (Fig. 1).

Shutter speed model (SSM)
Both the TM and its extended version (2.3.3) embed the implicit assumption that equilibrium transcytolemmal water exchange (between the intracellular space and extracellular extravascular space) is infinitely fast, or that the system is in, what is called, the fast exchange limit [29]. Water exchange between the intracellular space and the extracellular extravascular space effects the degree of T 1 shortening caused by CA. To account for this effect on the brain MRI signal amplitude, we applied the shutter speed model (SSM), which introduces a new pharmacokinetic parameter, the mean intracellular water molecule lifetime, τ i [29].
In the SSM, Eq. 5 is applied to the distribution of the CA in the blood, without assuming that the equilibrium transcytolemmal water exchange kinetics are in the FXL. The longitudinal relaxation rate is measured as: where b stands for the whole blood, p for the plasma and h the blood haematocrit. However, about half of the water in the blood is intracellular and cannot be accessed directly by the CA molecules [30]. The transport outside the erythrocytes therefore needs to be considered, as described by the two equilibria: The mean water molecule lifetime in commonly used CA is normally < 10 −7 s, and the linear Eq. 5 is suitable for homogeneous solutions. In the case of erythrocytes, Eq. 10 is also considered fast for some commonly measured [CA p ] values [30,31]. After extravasation, the CA is commonly distributed into the interstitial extracellular space, at a rate defined by: where R 1 *(t) is the rate constant of the extravascular water signal, r 10 is the interstitial CA relaxivity and p 0 is the fraction of the extracellular tissue water. The application of Eq. 11 to biological tissues assumes that the interstitium is a homogeneous solution and that the system remains in the fast exchange limit. However, many studies have shown that, even though the equilibrium in Eq. 10 is fast, the FXL assumption is not true for all [CA 0 ] values following a bolus injection [30]. [CA 0 ] depends on the dimensions of the parenchymal cells that are generally much larger than erythrocytes and have a less water-permeable cytolemmae. Furthermore, tissue parenchyma cannot be considered as a single homogeneous solution and a single MRI voxel will constitute a number of compartments. The main result of this compartmentalization is given by: where R 1L (t) is the long relaxation rate constant of the shutter speed model. R 1i is the H 2 O rate constant in the absence of exchange of CA and τ i is the average intracellular lifetime of a water molecule. The SSM was fitted by substituting Eq. 12 in Eq. 6 using the MATLAB functions fminsearch and fminsearchbnd, similarly to the TM. The initial estimates for the SSM K trans and v e were taken as the outputs of the TM, while the initial estimate for τ i was set at 0.1 [13]. The final K trans , k ep , v e and τ i maps were obtained from the fitting procedure with the best fit (lowest AIC value) as described in Fig. 1.  Fig. 1 The fitting procedure for K trans . A bounded and unbounded fitting were calculated together with the Akaike Information Criteria (AIC) map (AIC b and AIC u for the bounded and unbounded procedure, respectively). The final value of K trans , for each voxel of the map, was the one obtained from the function with the lowest AIC (k b when AIC b < AIC u and k u vice versa). The same procedure was carried out for each parameter

Extended Tofts model (ETM)
For highly perfused brain tissue, we applied the extended version of TM, which was introduced by Tofts in 1999 [32]. The ETMfits data with an additional parameter: the fractional plasma volume, v p [32]. This model is able to distinguish the effects due to contrast leakage from those due to intravascular contrast. Eq. 6 becomes: The ETM was fitted by substituting Eq. 13 in Eq. 6 using the MATLAB functions fminsearch and fminsearchbnd, similarly to the TM and SSM. The initial estimates for ETM K trans and v e were again taken from the output of the TM and the initial estimate for v p was set at 0.01 [13]. The final K trans , k ep , v e and v p maps were obtained from the fitting procedure with the best fit (lowest AIC) as described in Fig. 1.

Extended shutter speed model (ESSM)
The ESSM accounts for the contribution of the CA from the brain plasma compartment. This includes both v b , the fractional blood volume and τ b , the intravascular water molecule lifetime [13]. The contribution of the water signal comes from the three compartments (whole blood, EES and intracellular space) and is described by the matrix format of the Bloch equation [13]: where the column vectors are M = (M b , M o , M i ) and C = (M b0 R 1b , M o0 R 1o , M i0 R 1i ) with the 1 H 2 O magnetization vector M ≈ to the signal S. The exchange matrix X is given by: The subscripts b, o and i stand for blood, outside space and intracellular space, respectively. k bo (= 1/τ b ) represents the blood to interstitium transfer of water; k io (= 1/τ i ) the transfer of water from the intracellular space to the interstitium; k ob (proportional to 1/τ o ) the EES to blood transfer and k oi the EES to intracellular transfer [13]. For a SPGR sequence, the solution to Eq. 13 is the matrix form of the Ernst-Anderson relationship [33] which assumes that if the change in [CA] is relatively small during the acquisition, at every discrete data acquisition time point, the relaxation time can be estimated using: I is the identity matrix and S 0 (= (S b0 , S o0 , S i0 )) is the signal at baseline.
The ESSM was fitted by considering the measured signal E(t) as a combination of the signals in the three compartments (blood, outside space and intracellular space) [12] using: Furthermore, Eq. 16 was simplified to: where the column vectors are S = (S b , S o , S i ) and Þ . The model was fitted by substituting Eq. 18 into Eq. 17 using the MATLAB function fminsearch function. The outputs of the SSM model were used as the initial estimates for K trans and v o and one third of the measured signal was used as the initial estimate for S b0 , S o0 and S i0 . Furthermore the initial estimates for k bo , k ob , k oi and k io were taken from literature defined values as 1.2, 1.5, 1.1 and 1.2, respectively [13]. The following parameters were also derived: where f w is tissue volume fraction accessible to mobile aqueous solutes (assumed to be a constant and set to 0.8) and v i = 1 − (v b + v o ) [13].

No-effect model (NEM)
The no-effect model describes the case where [CA] within a brain voxel is so low that there is no permeability or vascular filling. In this situation, the MR signal is assumed to be unperturbed by the injection of the gadolinium-based CA and, as a consequence, the longitudinal relaxation rate does not change from its baseline value (R 10 ). The system is therefore described by this value at all times such that the data is fitted by the constant R 10 .

Arterial input function (AIF)
An additional ROI was drawn around the external carotid artery for the calculation of the image-derived AIF [34] (Fig. 2). Signal intensity curves were converted to R 1 -time curves by using the baseline signal intensity before the first pass of the CA as a reference, setting the haematocrit in the blood to 0.45, and getting the baseline blood T 1 from the T 10 map (Eq. 8).

Model comparison
A model comparison was carried out using the AIC to test for the best model in a given voxel. In particular, in the presence of exchange (where the NEM fails in the description of data), a voxel-wise comparison between models was carried out (with ETM and SSM being an extension of the TM, and the ESSM an extension of the SSM) to indicate which model provided the best fit using the AIC in each voxel. The selection method is shown in the flowchart in Fig. 3.

Stability of pharmacokinetic parameters
Once the model selection had been selected, the stability of each parameter within the selected models was evaluated in a simulation environment. Tissue curves were generated back from the extracted pharmacokinetic parameters and signal intensity curves were calculated with the inverse formula of Eq. 4. White Gaussian noise was added to the signal intensity curves using a signal to noise ratio (SNR) of 20. The SNR value for the simulated data was set by evaluating the SNR of the acquired data from the second and third phase of the dynamic acquisition sequence using the subtraction method [35]. The simulated noisy signal intensity curves were reconverted to noisy tissue concentration curves, and fitted to the selected pharmacokinetic model. This procedure was repeated 500 times for every kinetic parameter and the variability of each parameter was expressed in terms of coefficient of variation (CV): the percentage ratio between the standard deviation and the mean.

Model selection and parameter variability
The behaviour of each model was assessed by studying the quality of fit for each of the models. The input data, together with the fitted curves were normalized by the maximum value of the input data in order to compare results from the different fits. An example of a comparison of fit is shown in Fig. 4. For each tumour, a map with the result of the statistical comparison among models was built (Fig. 5). In this map, each colour represents the model for which the voxel-wise AIC value was lowest: The selected model of choice representing majority of low AIC was NEM (35.5% of voxels), followed by the ETM (32%), TM (28.2%), SSM (4.3%) and ESSM (< 0.1%). Fig. 5 shows the model selection map evaluated for two different lesions. Furthermore, final pharmacokinetic maps, for the K trans , k ep and v e , were built considering the model selection procedure. Within the final pharmacokinetic maps, each voxel was represented by the model with the best fit (lowest AIC) within that voxel. For each of these maps, the stability of each parameter and for each lesion was presented in terms of CV maps. Table 1 shows the results of the stability test on the final K trans map. The total number of voxel of the lesion, for each patient, together with the average K trans value evaluated over the whole tumour, is shown. Furthermore, the map was thresholded with a CV lower than 10, 20 and 50%, and the resultant percentage of preserved voxels (and their average value) is shown. Fig. 6 shows two final K trans maps with their CV map overlaid on them (A and D).

Discussion
In this study, we showed that in over one third of the brain tumours voxels (35.5%), standard model fitting of DCE-MRI data was inconclusive and therefore fitting these models to the data would lead to incorrect perfusion parameters. Considering a CVof 20%, only ≈ 25% of remaining voxels were found to be reliable. The reproducibility of this technique and, as a consequence, its reliability can be improved by improving the main sources of variability in quantitative DCE-MRI (the acquisition method and the quantification process) [16]. Nonetheless, in most studies, the intrinsic heterogeneity of the lesion is ignored by quantifying the perfusion with one single pharmacokinetic model and, more importantly, carrying out statistical analyses on one single whole tumour statistic (usually the average). In this study, we investigated the reliability of DCE-MRI focusing on the quantification analysis. We took into consideration the particularly heterogeneous nature of brain tumour vascular permeability due to the presence of the BBB, as well as existence of necrotic regions and provided a method to identify robust DCE-MRI data based on a model selection procedure and a stability test.

DCE-MRI models
MR scanners usually employ post processing perfusion tools which fit DCE data with the TM. This model (together with its extended version [32]) considers the system in a fast exchange limit [36,37], assuming an infinitely fast transcytolemmal water exchange between the EES and the intracellular space, The stability test consisted in the evaluation of maps of the coefficient of variation (CV). The table shows the stability of the final K trans maps. For each patient, the first two columns show the total number of voxels of the lesion and the average value of K trans evaluated across the whole lesion. Following, the percentage of voxels with a CV lower than 10, 20 and 50% and the relative mean K trans values are shown. The WHO grade and diagnosis of each lesion are also reported Each colour is representative of the model which best fitted the input data. An example of one slice of an enhancing (a, WHO grade IV) and non-enhancing (b, WHO grade II) lesion is shown which does not affect the overall signal decrease [36,37]. Therefore, many studies on the cell membrane water permeability coefficient have shown FXL to be physiologically unreasonable and inconsistent [38]. The shutter speed model was introduced to reflect a more realistic tissue environment. The model accounts for the intercompartmental water exchange effect, modelling this non-infinitely-fast exchange with the mean intracellular water molecule lifetime τ i . In 2005, Li et al. introduced a second generation of the shutter speed model which considers also a non-infinitely fast equilibrium transendothelial water exchange.

Model comparison and stability
The heterogeneity that exists in brain tumours means that one model is insufficient in explaining the different biologies that exist in different tumour regions. Multiple pharmacokinetic models are required for a complete description of the tissue. This variability is testified by the model selection procedure which showed how, in a single slice of one tumour, multiple models perform better. This result confirms the study of Bagher-Ebadian where they implemented a selection method based on nested models [3]. They found that in the necrotic core of the tumour, models describing vascular filling with no microvascular leakage (similar to the TM) and leakage without vascular reabsorption were selected because of the lack of blood flow. They also hypothesised that the model describing leakage with reabsorption (similar to the ETM) would be selected in the fast growing rims of the lesion. Our results show that there are a number of regions in the tumour where the CA exudation is prevented by the BBB and where the concentration of CA is so low that the evidence of perfusion is missing. In this case, the use of the NEM is recommended as the use of different models could result only in overfitting the data. In fact, our results showed that no leakage of the CA into the interstitium and the lack of flow of the CA through the tissue made the NEM the model of choice for the majority of regions, particularly in the non-enhancing lesions (37.5% of voxels). The result is very close to the ETM (32%), which was the model of choice in the enhancing lesions (54.8%). This suggests that, in areas where there is enhancement, a model with three parameters performs better and that the choice is dependent on the underlying state of the tissue. In fact, both the ETM and SSM are fitted by three parameters but the third parameter is very different between the two models (v p for the ETM describing a vascular component in the tissue, and τ i for the SSM describing the transcytolemmal water exchange). Furthermore, with the implementation of the ESSM, we saw that the transendothelial water exchange did not have any impact on the signal (compared with parameters derived by the fitting of simpler models). It is necessary to consider that the ESSM required nine parameters to be fitted and that the cost of fitting extra parameters is often contrary to the principle of parsimony. In fact, in fitting data to a noiselimited dataset, the estimation could be very poor and dependent on the optimization procedure itself (the initial conditions, for example) [3]. We compared the AIC values from the different fitting procedures to check whether a model with more parameters is more appropriate than a simple one. The ESSM was selected as the model of choice by < 0.1% of voxels, indicating that a model with three parameters performed better in the description of brain tumours and further confirming the poor quality of fit observed for the ESSM model. Our outcome agreed with the results of Duan et al. [39]. Using representative in silico and clinical (cervical cancer) DCE-MRI data, they demonstrated the sensitivity of complicated models (parameters > 3) to noise and their decreasing probability of being selected in low signal-to-noise data [39]. The reliability of DCE-MRI data is not only based on the goodness of fit of the chosen pharmacokinetic model, but also on the robustness of the extracted parameters. For this reason, we assessed, for each lesion and for each parameter, the coefficient of variation. We worked in a simulation environment where we added Gaussian noise to our signal and we fitted the noisy curves 500 times. This procedure resulted with a heterogeneous distribution of CVs that was not linked to contrast enhancement. In fact, Fig. 6 shows the plot of four different tissue activity curves together with the K trans value and its CV in an enhancing and non-enhancing tumour. The curves in Fig. 6b, c belong to the same enhancing lesion and while they correspond to similar K trans values of 0.81 and 0.52 [1/min], they varied by 12 and 128%, respectively. On the other hand, the curves in Fig. 6e, f belong to the same non-enhancing lesion and show regions with both a low (4%) and a high (97%) variability.
We set three different thresholds for the CV to evaluate the variability of the K trans . Table 1 shows the percentage of voxels and their relative mean K trans value at different CVs thresholds (10, 20 and 50%), for each lesion and for each patient. Higher grade glioma tend to have more voxels with a lower CVand also a more stable value of K trans while, for some of the other patients, the mean value of K trans is highly affected by the portion of voxels taken into consideration (P04, P07). This result confirmed the improper use of one average value in statistical comparisons of brain tumours. Not only because of the heterogeneity of the tissue under investigation but also, and more importantly, because it is affected by the reliability of fit within voxels.
Finally, Fig. 6 gives a graphical representation of this effect showing K trans values under the 20% threshold of CV covering only 25% of voxels (an average percentage value evaluated among all patients). This result suggested that only this selection of voxels represents robust values, which can be used in the following statistical analyses, as, more importantly, in clinical evaluations. The selection of the threshold that makes DCE-MRI robust is, however, dependent on the effect size that is being measured and hence will vary across studies.
The main limitation of this study is the small size of the dataset. Furthermore, the sensitivity of DCE-MRI data to water exchange effect was reduced by the 26°flip angle acquisition (exchange-minimized approach) [40]. As a consequence, the precision of the τ i parameter extracted might be low.
In conclusion, DCE-MRI methods hold great promise for quantitative in vivo evaluation of permeability and vascular properties under different pathophysiological conditions. It allows us to identify, and quantitatively measure, smaller changes in permeability for pathological conditions effecting the BBB, than would be observed through visuassessment of post-contrast T 1 -weighted images. Different models yield different pharmacokinetic parameters and, for this reason, a model selection is critical for the appropriate analysis of DCE-MRI time courses based on the regional tissue biology, specifically permeability and vasculature. Future work needs to assess the physiological basis for each model in the reliable selection of DCE-MRI data. The applicability of each model depends on the physiology, anatomy and heterogeneity of the tumour and the tumour microenvironment. In addition, due to the noisy nature of DCE-MRI data, a model selection procedure alone is not enough: pharmacokinetic parameters need to be validated with a stability test in order to give only robust results for statistical analyses and clinical evaluation.
Funding This work was funded in part by the National Institute for Health Research (NIHR) Imperial Biomedical Research Centre (BRC), the Brain Tumour Charity and the Brain Tumour Research Campaign.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Ethical approval The study was reviewed and approved by the London-Fulham Research Ethics Committee.
Informed consent Written and informed consent was obtained from all participants before recruitment to the study and all data was anonymised in accordance with the EU General Data Protection Regulation.