Tracer kinetic model selection for dynamic contrast-enhanced magnetic resonance imaging of locally advanced cervical cancer.

Abstract Background. Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) offers a unique capability to probe tumour microvasculature. Different analysis of the acquired data will possibly lead to different conclusions. Therefore, the objective of this study was to investigate under which conditions the Tofts (TM), extended Tofts (ETM), compartmental tissue uptake model (C-TU) and 2-compartment exchange model (2CXM) were the optimal tracer kinetic models (TKMs) for the analysis of DCE-MRI in patients with cervical cancer. Material and methods. Ten patients with locally advanced cervical cancer (FIGO: IIA/IIB/IIIB/IVA – 1/5/3/1) underwent DCE-MRI prior to radiotherapy. From the two-parameter TM it was possible to extract the forward volume transfer constant (Ktrans) and the extracellular-extravascular volume fraction (ve). From the three-parameter ETM, additionally the plasma volume fraction (vp) could be extracted. From the three-parameter C-TU it was possible to extract information about the blood flow (Fp), permeability-surface area product (PS) and vp. Finally, the four-parameter 2CXM extended the C-TU to include ve. For each voxel, corrected Akaike information criterion (AICc) values were calculated, taking into account both the goodness-of-fit and the number of model parameters. The optimal model was defined as the model with the lowest AICc. Results. All four TKMs were the optimal model in different contiguous regions of the cervical tumours. For the 24 999 analysed voxels, the TM was optimal in 17.0%, the ETM was optimal in 2.2%, the C-TU in 23.4% and the 2CXM was optimal in 57.3%. Throughout the tumour, a high correlation was found between Ktrans(TM) and Fp(2CXM), ρ = 0.91. Conclusion. The 2CXM was most often optimal in describing the contrast agent enhancement of pre-treatment cervical cancers, although this model broke down in a subset of the tumour voxels where overfitting resulted in non-physiological parameter estimates. Due to the possible overfitting of the 2CXM, the C-TU was found more robust and when 2CXM was excluded from comparison the C-TU was the preferred model.

the shape of the AiF and the tissue enhancement curve it is possible to extract information about these tissue-specific parameters through tracer kinetic modelling (tKM) [1].
Most types of cancers exhibit a chaotic vasculature once they grow beyond a diameter of 2 mm [2] making DCE imaging a unique tool for investigating cancer progression. A further consequence of tumour growth is a heterogeneous supply of blood, which indirectly may cause hypoxic tumour sub-regions [3,4] resulting in radio resistance [5,6] and more aggressive tumour phenotypes [7,8]. Clinical and preclinical evidence linking DCE-Mri and hypoxia has been inferred for xenografts [9] and patients with advanced cervical cancer using a hypoxic gene profile [10]. An earlier study using the Eppendorf probe to measure partial oxygen pressure reported similarly that low enhancement was related to low oxygenation [3,11]. other studies have also showed the value of DCE-Mri to predict local tumour control and overall patient survival [12][13][14].
the analysis of DCE data has changed gradually during the last decades towards increasing complexity along with improvements in Mr equipment. However, tKMs are often applied disregarding the conditions for which they have been developed. For example, the early brix model assumes constant infusion of the contrast agent [15], and the tofts model (TM) and extended tofts model (EtM) assumes that the capillary mean transit time (Mtt) of the contrast agent is negligible [16,17]. similarly, the compartmental tissue uptake model (C-tU) assumes that the concentration of contrast agent is much higher in plasma compared to the in extracellularextravascular space and the two-compartment exchange model (2CXM) presupposes the existence of two distinct compartments. these are all idealised scenarios and are likely not reflecting the true biology. the question arises as to which model is the optimal for quantification of DCE data. Here the term 'optimal' is used to mean the model that best characterised the data. Concerning cervical cancers, a DCE model comparison study has been carried out comparing the tM, EtM, C-tU and the 2CXM by Donaldson et al. [18]. they concluded that the 2CXM and C-tU were better suited for cervical tumour DCE-analysis than the tM and the EtM. their model comparison study was, however, performed on tumour-averaged enhancement curves although it was realised that heterogeneous enhancement was present within the tumour volume. to our knowledge, no studies have been undertaken to account for tumour heterogeneity when comparing different models for quantification of DCE for cervical tumours.
Accordingly, the aim of this study was to investigate the conditions under which the tM, EtM, C-tU and 2CXM are the optimal tKMs for the voxel level analysis of DCE-Mri in patients with cervical cancer.

Material and methods
Patients this prospective study was approved by the local medical ethics research board, and all patients gave written informed consent. in brief, 10 cervical cancer patients with advanced stages (Figo: iiA/iib/iiib/ iVA -1/5/3/1) underwent DCE-Mri prior to radiotherapy. Median (range) age was 54  years and the median (range) tumour size was 90.6 (40.2-202.7) cm 3 .

Imaging protocol
DCE-Mri was performed using a 3t Philips Achieva scanner and a three-dimensional (3D) saturation recovery spoiled gradient echo technique with 20-24 slices having a 5 mm slice thickness, tr/tE of 2.9/1.4 ms, t sat of 25ms, flip angle (FA) of 10°, in-plane resolution 2.3 mm  2.3 mm and 2.1 s time resolution. the bolus injected was 0.1 mmol/kg Dotarem at 4 ml/s, followed by a 50 ml saline flush. A total of 120 dynamic scans were obtained of which 18 time-points were scanned before the bolus arrived at the iliac arteries. A t 1 relaxation map was constructed before contrast injection using a 3D gradient recalled echo sequence with five different FA scans (5°, 10°, 15°, 20°, 25°) with the same orientation and field of view as the dynamic scan.
Image processing the t1 map was produced following the work by Fram et al. [19]. this was afterwards used to convert the dynamic scan from signal intensity to contrast agent concentration as described earlier [20]. to increase the signal to noise ratio, we spatially convolved all volumes with a 3D block kernel of 0.19 cm 3 .

TKM fitting
Four of the most often used kinetic models developed for fast infusion protocols in patients with cervical cancer were investigated: tM, EtM, C-tU and 2CXM [16,21]. From these tKMs, it was possible to extract information about the forward volume transfer constant (K trans ), the extracellular-extravascular volume fraction (v e ), the plasma volume fraction (v p ), the blood flow (F p ) and permeability-surface area product (Ps) (supplementary Appendix, available online at http://informahealthcare.com/doi/abs/ 10.3109/0284186X.2014.937879). All model parameters were calculated using a non-linear least squares technique that minimises the sum squared difference between model fit and data. K trans , F p and Ps were constrained to be positive while v p and v e were constrained between 0 and 1.Voxels where the solution space meet the limits set by the constraints (e.g. v e  1) were excluded. Furthermore, voxels where the fit of the bi-exponential 2CXM was reduced to a mono-exponential model were excluded from the comparison of parameters from different models (supplementary Appendix available online at http://informahealthcare.com/doi/abs/10.3109/ 0284186X.2014.937879).
in each patient, an AiF was derived from the non-convolved dataset by averaging the AiF measured in the left femoral arteries where the b 1 field was consistently most homogenous. Due to inflow artefacts in the femoral arteries, the pre-contrast longitudinal relaxation (t 1,0 ) determination is connected with some uncertainty, and a literature value for t 1,0 was chosen instead: t 1,0 (blood)  1660 ms [22]. to correct for differences in large and small vessel haematocrit, the AiF was multiplied by a factor 1.18, based on an assumed haematocrit of 0.38 and an assumed small to large vessel ratio of 0.7 [23].
When more complex models are applied, data with better signal to noise ratios (sNr) are needed [24,25]. similar to luypaert et al. [25], the contrast to noise ratio (CNr) was defined as the ratio of the peak contrast enhancement to the standard deviation of signal intensities at the baseline before contrast injection.
the region of interest (roi) was chosen to be the clinical gross tumour volume (gtV) delineated by an experienced oncologist on a transversal t2-weighted Mri following the recommendations of the gEC-Estro working group [26].

Model comparison
For the comparison of the four tKMs, we used the same approach described by brix et al. [27] and Korporaal et al. [23]. the four tKMs described above form a set of candidate models to quantify DCE data. We consider the optimal model as the model with the smallest possible number of parameters for adequate representation of the data. therefore, for each voxel, Akaike information criterion (AiC) was calculated and corrected for small sample sizes (AiCc) [28]. based on the AiCc values, for each voxel, the optimal model was chosen by selecting the minimum AiCc (AiCc min ). to assess the impact of a particular choice, we investigated the interchangeability of the K trans values between the different models. if one uses the 2CXM model for analysis, a K trans value can be derived from the parameters F p and Ps or a combination thereof. similarly, when using one of the tofts models, the parameter K trans reflects F p , Ps, or a combination thereof [17]. this depends on the conditions under which the contrast agent leaks out of the vessels into the extravascular-extracellular space. Under flowlimited conditions (Ps F p ) K trans equals F p . Under permeability-limited conditions (Ps F p ) K trans equals the Ps, and in the mixed flow-and permeability condition in between, K trans can be derived: E x F p , where E is the extraction fraction (Appendix available online at http://informahealthcare.com/ doi/abs/10.3109/0284186X.2014.937879). We chose to analyse all three conditions for the subset of voxels in which tM, EtM and 2CXM performed best. in this way, we determined how a K trans value should be derived from the tKM parameters of the 2CXM model and how K trans values from the tofts models should be interpreted. the tM, EtM and C-tU are all special case solutions of the more general 2CXM. the interpretation of the hemodynamic parameter of the tM is valid only if the underlying tissue is weakly perfused (F p ) or weakly vascularised (v p ). similarly, for EtM the assumptions are highly perfused (F p ) or weakly vascularised (v p ) [17]. the validity of these assumptions may be investigated by calculating the plasma Mtt which may be derived from the central volume theorem. For the 2CXM, the plasma Mtt is the ratio of the plasma volume (v p ) over the sum of all influxes (F p  Ps) [18].

Statistical analysis
Pearson's correlation and Wilcoxon's rank sum test were applied for statistical inference. the correlations were compared using Hotelling's t 2 -test. the statistical significance level was chosen as p  0.05. Figure 1 shows an overview of all parameter maps that were analysed for a single transverse slice in a single patient (corresponding relative uncertainty for the parameters maps are shown in supplementary Figure 1, available online at http://informahealth care.com/doi/abs/10.3109/0284186X.2014.937879). Comparing the individual voxels in the model respective AiCc maps, an AiCc min map could be generated reflecting the model that yielded the lowest AiCc value. Figure 2 shows that all four tKMs are the optimal models in different contiguous regions of the cervical tumour. For individual patients, the tM was the optimal in 6.4-34.9%, the EtM in 0.8-7.9%, the C-tU in 9.8-39.3% and the 2CXM in 40.3-82.8%. For all the 24,999 voxels analysed, the tM was optimal in 17.0%, the EtM was optimal in 2.2%, the C-tU was optimal in 23.4% and the 2CXM was optimal in 57.3%. When excluding the 2CXM from the comparison, the distribution was as follows: the tM was optimal in 17.5%, the EtM was optimal in 12.7% and the C-tU was optimal in 69.8% (supplementary Figure 2, available online at http://informahealthcare.com/doi/abs/10.3109/ 0284186X.2014.937879).

Results
the percentage of voxels that were discarded from the analysis was 29.3%. Due to overfitting of the 2CXM, 12.5% of the voxels were removed while the remaining were discarded because either v e (tM), v e (EtM) or v p (C-tU) were constrained at upper limits set [i.e. v e (tM)  1 or v e (EtM)  1 or v p (C-tU)  1]. the overfitting of the 2CXM occurred predominantly in the regions where tM was optimal (9.5%) and less frequently in regions C-tU was optimal (3.0%).
Examples of tissue enhancement curves under which the different models were optimal are shown in supplementary Figure 3 (available online at http://informahealthcare.com/doi/abs/10.3109/ 0284186X.2014.937879). in the case, where tM was optimal, there was very little difference in the fits of the four tKMs and, therefore, the one with the lowest number of free parameters will have the lowest AiCc. similarly, for where the EtM was optimal, the EtM and 2CXM models show similar fits, and the model with the lowest number of free parameters was again chosen as optimal. in this case the tM and C-tU were in unable to model the steep upslope. Where the C-tU model was optimal, the 2CXM shows a similar fit while tM and EtM have discrepancies at the initial peak and during the wash-out phase. in the final case, where the 2CXM was optimal, only this model was able to adequately describe the shape of the tissue enhancement curve.
there was a significant difference in the median values of the CNr between the regions where the four tKMs were optimal (table i). the lowest median CNr was found for optimal tM followed by optimal EtM and optimal C-tU, while the highest median CNr was found where 2CXM was optimal. table i also showed that the Mtt varied little across the different optimal regions for all models and was 5-6 times larger than the temporal resolution of the dynamic scan. From table i, the ratio F p /Ps was found to be greater than 1 in 86% of the voxels and the quantitative estimate of v e (2CXM) was smaller than that of the tM and EtM by 28%. the interchangeability of the K trans was assessed in Figure 3, showing the different distributions of the calculated and derived K trans values based on the different flow and permeability conditions in the cervix. the calculated K trans was directly estimated from either tM or EtM while the derived K trans was estimated using the 2CXM. in the range where tM was optimal, a high correlation was found between K trans (tM) and F p , r  0.91. similarly, where the 2CXM model was optimal, there was again a strong correlation between K trans (tM) and F p , r  0.90. For EtM, K trans was also correlated better with F p compared with the other regimes  r  0.86. All three correlations under flow-limited regime were significantly higher as compared with the mixed flow and permeability regime.
both K trans (tM) to K trans (EtM), and v e (tM) to v e (EtM), were highly correlated with Pearson's r  0.96 and r  0.99, respectively (Figure 4a-d). A poorer correlation was observed between the estimates of v e (2CXM) to v e (tM) and v e (EtM), with r  0.75 and r  0.78, respectively (Figure 4e,f). No correlation was observed between the plasma volume estimate of the EtM and 2CXM (Figure 4g,h). However, correlation between v p (2CXM) to v e (EtM) and v e (tM) was observed, with r  0.66 and r  0.69, respectively (Figure 4i,j). When comparing v p v e (2CXM) with v e (tM) and v e (EtM), a high degree of correlation was found (Figure 4k,l). Comparing parameter estimates from C-tU and 2CXM showed good agreement between the estimates of flow (F p ) and significant better correlation between v p estimated by C-tU and 2CXM compared to that estimated by EtM and 2CXM. No correlation was identified between Ps estimated by C-tU and 2CXM (supplementary Figure 4, available online at http://informahealthcare.com/doi/abs/10.3109/ 0284186X.2014.937879).

Discussion
this study addressed the question: under what conditions do the tM, EtM, C-tU and 2CXM perform optimally on a voxel level. None of the models were generally optimal across all voxels, but were optimal in individually distinct subregions. in regions with high CNr, the 2CXM outperformed the tM, EtM and C-tU. this result was supported by earlier simulations studies [24,25] investigating uncertainties in parameters estimated at different levels of CNr. the higher requirements of CNr for successful use of the 2CXM supports the use of scanners with higher field strengths and, therefore, a similar analysis performed at lower field strengths will likely skew the optimal model selection distribution towards the simpler models. AiCc min in all patients reflected a heterogeneous pattern, with contiguous regions in which one of the models outperformed the others (Figure 1). the fact that the regions were contiguous suggests that this was not caused by noise, but rather reflected different underlying properties of the tissue. the 2CXM was most often the optimal model for describing the DCE-Mri data for cervical cancers with 57.3%, followed by C-tU with 23.4%, tM with 17.0% and the EtM was only optimal in 2.2%. this is somewhat contradictory to common practice, since the extended tofts model has been avidly applied to cervix data, justified by the fact that the influence of the intravascular compartment was non-negligible [29]. the results of the 2CXM supported the notion that the intravascular compartment cannot be neglected. Figure 3 shows that all regions had high correlation between F p and K trans , as estimated by either of the tofts models. this suggested that K trans mainly reflects flow. theoretically, K trans will reflect flow when Ps is infinite. However, table i showed that Ps F p in practically all regions of the tumour, suggesting that K trans should mainly reflect Ps. these contradictory observations suggest that the assumptions of the tMs were not supported by the underlying physiology. the assumption leading to the special case solution with dynamic signal profile similar to tM for non-zero flow is that v p is equal to zero. given this assumption, K trans becomes F p if Ps is infinite. However, a different special case solution of a similar shape as tM can be deduced for non-zero flow assuming: v e ≠0≠v p and Ps approaching zero [17]. in this case, the impulse response function is equal to F p * exp(-F p /v p ) where v p and v e have been exchanged compared to tM. this  it should be noted that although the 2CXM seems to predominantly be the optimal model, it should be used with caution. the double exponential decay of the 2CXMs impulse response function allows for 'overfitting' of the data that is not necessarily supported by the underlying physiology. the 12.5% of voxels that were removed from the analysis were removed because the data only supported a mono-exponential decay and the fit parameters of the 2CXM, thus could not be reliably determined. may also explain why v p (2CXM) is correlated to v e (EtM) and v e (tM). Finally, plasma Mtt was considerably longer than that of the temporal resolution of the acquired DCE dataset and was therefore not negligible as assumed by tM and EtM. this finding was also supported by the work done previously by Donaldson et al. [18].
the lower estimate of v e of the 2CXM compared to the tM or EtM (table i) was also reported by Donaldson et al. [18]. the significant increase in correlation between v e (tM) and v e v p (2CXM) compared with v e (tM) and v e (2CXM), suggests that v e  where 2CXM was optimal (d-f). Comparison between v p (2CXM) and v p (EtM) where EtM and 2CXM were optimal (g-h). Comparison between v p (2CXM) and v e (EtM) where EtM and 2CXM were optimal (i-j). Correlation between v p (2CXM)+v e (2CXM) to v e (EtM) and v e (tM) where 2CXM was optimal (k-l). therefore, from a practical perspective it may be more feasible to use C-tU model for cervical cancer patient, since fever parameters, limits the possible solution space, thus making it more robust.
Although the local tumour control is high in patients with cervical cancer is [30], the use of functional imaging like DCE-Mri could further improve on local control, particular in tumours responding poorly to radio-chemotherapy. recommendations regarding the analysis of the DCE-Mri data and what parameters should be reported for outcome evaluation have previously been given by leach et al. [31,32]. they recommended reporting primarily K trans or the integrated area under the concentration time curve (iAUgC). since we showed that K trans primarily reflects flow (F p ) for cervical cancer it seems plausible to report this parameter if the 2CXM or C-tU is chosen as the preferred model for analysis. generally, it is difficult to recommend one model rather than another, since there are many factors that may influence this choice, including CNr and the tissue physiology. it does, however, seem prudent to use the C-tU instead of the EtM, and if tM is used, one should be aware that the volume extracted is likely not the extravascular-extracellular volume for cervical cancer patients, but rather the entire distribution volume.
in conclusion, all four models provided reasonable fits to DCE-Mri data from the cervix. For most voxels, the 2CXM was best in describing the contrast agent enhancement of pre-treatment cervical cancers. Due to the redundant nature of the 2CXM, it may, however be more useful to apply a more stable model like the C-tU. the assumptions for correct interpretation of tM and EtM seems to be violated for cervical cancer. Finally, the data shown here supported that K trans as estimated by tM or EtM throughout the tumour mainly reflect flow.