Assessing Treatment Response Through Generalized Pharmacokinetic Modeling of DCE-MRI Data.

Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) enables the quantification of contrast leakage from the vascular tissue by using pharmacokinetic (PK) models. Such quantitative analysis of DCE-MRI data provides physiological parameters that are able to provide information of tumor pathophysiology and therapeutic outcome. Several assumptive PK models have been proposed to characterize microcirculation in the tumoral tissue. In this paper, we present a comparative study between the well-known extended Tofts model (ETM) and the more recent gamma capillary transit time (GCTT) model, with the latter showing initial promising results in the literature. To enhance the GCTT imaging biomarkers, we introduce a novel method for segmenting the tumor area into subregions according to their vascular heterogeneity characteristics. A cohort of 11 patients diagnosed with glioblastoma multiforme with known therapeutic outcome was used to assess the predictive value of both models in terms of correctly classifying responders and nonresponders based on only one DCE-MRI examination. The results indicate that GCTT model's PK parameters perform better than those of ETM, while the segmentation of the tumor regions of interest based on vascular heterogeneity further enhances the discriminatory power of the GCTT model.

Introduction

Magnetic resonance imaging (MRI) is the most commonly used medical imaging technique to demonstrate tumor morphology and the relationships between malignant lesions and the neighboring structures. MRI also provides relevant clinical information for clinical management and surgical planning. In recent years, sequential sets of MRI data and the development of small molecular weight paramagnetic contrast agents have had a major impact in the assessment and monitoring of tumor treatment response.1

The blood–brain barrier (BBB)2 is a membrane that separates the parenchyma of the central nervous system from the blood. It consists of endothelial cells interconnected by tight junctions that cause a selective permeability based on molecular weight and osmotic characteristics. In pathologies such tumors,3 multiple sclerosis,4 and acute ischemic strokes,5 the BBB is disrupted by various mechanisms. The BBB disruption is reflected by, for example, MRI contrast enhancement in pathological areas.

Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is a useful imaging technique for assessing the BBB leakage; it is based on the extravasation of the contrast agent (CA) from arteries to the parenchyma, which leads to decreased longitudinal relaxation time and therefore increased MRI signal intensity. It relies on fast T1-weighted MRI sequences obtained before, during, and after the intravenous administration of a gadolinium (Gd)-based CA.

There are several models for the quantification of DCE-MRI data that can be applied in many anatomical areas such as brain, prostate, and breast tissue, aiming to extract tissue physiological parameters. These models provide independent markers of angiogenic activity, and therefore act as diagnostic and prognostic indicators in a broad range of tumor types. A number of physiological parameters of vascularization coming from pharmacokinetic (PK) models play an important role in quantitatively assessing tumor angiogenesis within a region. Neovascularity and leaky vessels are associated with higher k (volume transfer constant between plasma and EES) values, while response to therapy has been correlated with a drop in ktrans values between imaging sessions.6 The results can be combined with other quantitative measures such as T2* perfusion7 in brain and H215 O-PET8 in body imaging for the evaluation of treatment response. In this paper, we study the applicability of PK models for the early characterization of therapy response in glioblastoma multiforme (GBM) patients and also present a novel framework for enhancing the imaging biomarkers by pre-segmenting the tumor region of interest (ROI) into subregions according to vascular heterogeneity metrics.

Methods

MRI acquisition protocol

All MRI acquisitions were performed on a 3-T Philips scanner. Prior to the CA injection, two acquisitions were made for the variable flip angle (VFA) data using flip angles (FA) of 5° and 15°, repetition time (TR) of 10 ms, and echo time (TE) of 2 ms. DCE-MRI data were acquired using fast three-dimensional spoiled gradient echo (SPGR) with TR 3.6 ms, TE 1.75 ms, and FA 6°. The dimensions of the image was 192 × 192 voxels with 36 slices and 50 time points for DCE-MRI protocol, temporal resolution 6 s, and the size of one voxel 1.15 × 1.15 × 2.99 mm. The type of CA that was used was Dotarem, and the quantity was 0.1 mmol/kg of body weight.

Preprocessing

All datasets (both DCE-MRI and VFA) were co-registered to the arterial phase [maximum signal-to-noise ratio (SNR)] of the DCE-MRI, and a temporal smoothing algorithm was subsequently used to remove the intrinsic noise of the MRI signal. Due to the large discrepancies in the values of arterial concentrations among different subjects, we used a theoretical arterial input function (AIF) for all subjects. To this end, biexponential decay was assumed to describe the plasma concentration using values from earlier work.9

Estimation of contrast agent concentration

To proceed to the quantification of CA kinetics from signal intensities (SIs), the concentration of CA should be determined at each time point of the dynamic scan; possibly the most crucial step. A first approach is to assume that the CA concentration is proportional to the change in SI. However, when CA concentrations are high, the relationship between concentration and SI becomes nonlinear10 and the previous approximation may lead to significant errors.

There are several techniques for measuring changes in T1 due to CA presence, such as inversion recovery,11–13 Look-Locker,14,15 and VFA.3,16–18 The VFA method was chosen in our work because of the good SNR and reduced acquisition times.

The CA concentration is related to the change in relaxation time via the following formula: where r1 is the longitudinal relaxivity of the CA, T10 is the longitudinal relaxation time without CA, and T1(t) is the longitudinal relaxation at time t after the injection.

The MRI signal S1 acquired with TE ≪ T2 is given by the following equation: where S0 is the relaxed signal for a 90° pulse when TR ≫ T1, α is the flip angle, and TR is the repetition time. Substituting in Equation (2) the measurements acquired from the VFAs data and solving the nonlinear problem per voxel, the vector [S0, T10] is calculated.

Substituting in Equation (2)S1 with S(t) (the SI in DCE-MRI data) and α with a (flip angle that was used in DCE-MRI protocol), and solving for T1, the time course of the longitudinal relaxation time can be calculated from Equation (3).

The overall procedure is depicted in Figure 1, where using the VFAs data, vector [S0, T10] is estimated using Equation (2). Afterward, using DCE-MRI data, the time course of the longitudinal relaxation time (T1(t)) is calculated by Equation (3). Substituting the pre-contrast relaxation time (T10) and the time course of the longitudinal relaxation time (T1(t)) in Equation (1); the CA concentration is calculated.

Figure 1

Schematic diagram for converting signal intensity into concentrations. The rectangle shows the start point, the trapeziums show the parameters, and the circle includes the equations to be solved.

Impulse response function

Using systems theory and assuming that biological tissues are linear, time-invariant systems, tracer kinetics can be described by the impulse response function (IRF),19 which incorporates the physiological parameters at each voxel. This way, the tissue concentration is given by Equation (4): where IRF(t) = IRF(t) + IRF(t), with IRF the vascular impulse response function and IRF the parenchymal impulse response function, F is the flow, C is the concentration in the tissue, C is the concentration in the arterial blood (AIF), and ⊗ represents convolution.

Tofts and extended Tofts model

The most commonly used model in literature is the Tofts model (TM),20 which is a single-compartment model in which CA diffuses from an external vascular space into a well-mixed tissue compartment. Tofts et al assumed that when CA is injected to the bloodstream, it will pass the disrupted blood vessel endothelium and move to the extravascular extracellular space (EES) with a rate proportional to the difference of CA concentration between the plasma (C) and the EES space (C(t)):

Using the convolution theorem, the solution of Equation (5) is given by the following formula: with IRF = 0 and IRF . In the above equations, k represents the volume transfer constant from the plasma space to EES, v is the volume of EES, and k=k/v is the transfer constant from EES to the plasma space. The negligible plasma volume assumption of TM is invalid for several tissue types, especially for brain tumors, which might lead to significant errors.

Tofts et al extended the original model by introducing the vascular term as an external compartment. The result was to separate the enhancement caused by contrast leakage from that caused by intravascular contrast. The extended Tofts model (ETM)21 is described by the following equation: where IRF = v·δ(t)/F and IRF , and v is the volume of vascular space. Given that C and C are known by converting the tissue and the artery SIs, respectively, and using Equation (7), the vector [k] is estimated per voxel. The weak point of these models is that k can be interpreted either as plasma flow in flow-limited cases or as tissue permeability in permeability-limited cases, but does not allow separate estimation of these two independent parameters. Moreover, TM can provide accurate PK parameters if and only if tissue is weakly vascularized, while ETM is also accurate in highly perfused tissues.22

Gamma capillary transit time

The gamma capillary transit time (GCTT) model23 unifies four models: TM,20 ETM,21 the two-compartment exchange (2CX) model,24 and the adiabatic tissue homogeneity (ATH) model.25 A major drawback of the aforementioned models is that every voxel is treated as a single capillary tissue unit with a single capillary transit time. The distributed capillary adiabatic tissue homogeneity (DCATH) model26 overcame this drawback by assuming a statistical distribution (normal, corrected normal, and skewed) of the transit times in the parenchyma and vascular IRFs. However, the DCATH model failed in the sense that certain results did not correspond to realistic values (eg, negative transit times) and the model could not provide closed-form solutions.23

The GCTT model overcame the limitations of the DCATH model by assuming that capillary transit times are governed by the gamma distribution. This way, each voxel is assumed to have different characteristics that are described by the parameter α−1, which is defined as the width of the distribution of the capillary transit times within a tissue voxel, and actually expresses the heterogeneity of tissue microcirculation and microvasculature.19,23 Depending on this para meter, IRF is adapted to the specific properties of the tissue voxel. Mathematically, the parameter α−1 is expressed by the following equation: where τ is the scale parameter of the gamma distribution, and t is the capillary transit time.

The vascular and parenchymal components of the IRF in the GCTT model are given by the following equations: where D(u) is the gamma distribution of capillary transit times, E is the extraction fraction, which indicates the fraction of CA that is extracted from v into v in a single capillary time, γ is the upper incomplete gamma function, and k is the CA transfer rate from EES to the vascular space. Replacing Equations (9) and (10) in Equation (4), the formula for the GCTT model can be derived as:

As in ETM, C(t) and C(t) are computed by converting the tissue and the artery SIs, respectively, and then estimating using Equation (11) the vector [F, E, k−1] per voxel. As was mentioned earlier, the GCTT model converges to preexisting models in certain limits of its parameters, as shown in Table 1.

Table 1

GCTT convergence to other models.

CONDITION	CONVERGENCE
α−1→0	IRFATH→ IRFGCTT
α−1→1	IRF2cx→IRFGCTT
α−1→1, τ→ε	IRFETK→IRFGCTT
α−1→1, τ→0	IRFTK →IRFGCTT

Heterogeneity tumor region segmentation based on 0−1

In order to enhance the application of the GCTT model in real clinical data, a preprocessing step is proposed by exploiting the α−1 parameter. For this purpose, the MR image of the tumor was segmented based on each voxel’s α−1 value, in order to separate tumor into subregions of similar vascular heterogeneity characteristics. After extensive experimentation and observation of the α−1 histogram characteristics in the tumor ROIs, a four-subregion segmentation scheme was proposed, as shown in Figure 2. It was observed that, in most of the cases, the histogram distributions of the α−1 parameter include three characteristic peaks and a plateau in the same value intervals. The α−1 boundary values were determined empirically based on the average histogram profile. The first subregion consists of the first peak, the second subregion consists of the linear histogram part, the third subregion consists of the second peak, and the fourth subregion consists of the third peak. Provided that the average profile is similar in future studies, it is possible to use the same values and test by other researchers.

Figure 2

Proposed methodology for separating the tumor image area into four subregions based on the vascular heterogeneity characteristics (α−1). This way, in each subregion as well as in the whole tumor region, the value of all PK GCTT parameters can be assessed separately.

The different subregions are defined as follows:

“Full Homogeneous” subregion: when α−1 ∈ (0,0.2] – in tables of results referred to as 1;

“Homogeneous” subregion: when α−1 ∈ (0.2,0.5] – in tables of results referred to as 2;

“Heterogeneous” subregion: when α−1 ∈ (0.5,0.85] – in tables of results referred to as 3;

“Full Heterogeneous” subregion: when α−1 ∈ (0.85,1) – in tables of results referred to as 4.

This segmentation comprises a new framework which in essence divides each tumor image region into four subregions of differing homogeneity as given by the GCTT model. For each subregion, the PK parameters were collected, and the classification results for the two groups of patients were reported. For completeness, in the comparative study the results in the entire tumor region (“All regions” – in tables of results referred to as 0) were also computed in order to test whether the proposed heterogeneity segmentation framework can add value to the strength of the PK GCTT candidate imaging biomarkers.

Implementation of the models

Figure 3 presents the workflow for the conversion of concentrations to PK parameters. The PK parameters were calculated per voxel via nonlinear least squares problems†, by solving Equation (7) for ETM and Equation (11) for the GCTT model.

Figure 3

Workflow for calculation of PK parameters. The rectangles show the start and end points, the trapeziums show the extracted PK parameters, and the ellipsoids include the equations to be solved.

In ETM, all parameters were assumed positive, and the initial values of k, k, and v were 0.001 (min−1), 0.009 (min−1), and 0.01 (none), respectively.

In GCTT, all parameters were assumed positive and the initial values of F, E, k, τ, and α−1 were 0.002 (min−1), 0.5 (none), 0.009 (min−1), 1.5 (s), and 0.5 (none), respectively. In the last step, the additional parameters (k, PS, t, and v) were calculated by the following equations: where PS is the permeability surface area product per unit mass of tissue and Hct is the hematocrit (reference value 0.45).

Clinical question for assessing and comparing the PK models

We retrospectively analyzed data from 11 patients (anonymized) diagnosed with GBM from University Hospital of Leuven at a time point before definite diagnosis regarding the therapy outcome. As shown in Table 2, each patient belonged in one of the following three categories: response, immune reaction, and progression, depending on the actual final outcome. Regarding treatment outcome, patients in the first two categories exhibited similar characteristics; thus for testing the discriminatory power of the PK models in this work, we divided patients into two groups: “Response” and “Non-Response”, as shown in Table 2. This decision mainly reflects the most important aspect of this work, ie, to use PK imaging biomarkers for assisting the clinician to early assess the treatment outcome and better manage the therapeutic alternatives.

Table 2

Clinical cases.

PATIENT	OUTCOME	GROUP
P1	Immune reaction	Response
P2	Immune reaction	Response
P3	Response	Response
P4	Response	Response
P5	Progression	Non-Response
P6	Progression	Non-Response
P7	Progression	Non-Response
P8	Progression	Non-Response
P9	Progression	Non-Response
P10	Progression	Non-Response
P11	Progression	Non-Response

The data in this study are part of a larger longitudinal imaging study in which patients with GBM treated with immune therapy receive monthly an extensive MRI exam. The institutional review board of the University Hospitals Leuven approved this study and informed consent of every patient was obtained before participation. This was done in order to characterize the temporal pattern on imaging in patients treated with a novel therapy and, secondly, in case an immune reaction was expected, to be able to characterize this with advanced MRI. Clinical outcome of the therapy is thus based on imaging follow-up and clinical examination. The standard imaging guidelines were used to define progression, response, and immune reaction, as defined in the Response Assessment in Neuro-Oncology (RANO) criteria.27 These criteria are currently defined as the standard clinical guidelines and count as good clinical practice. Performing histopathological examination obtained by biopsy or craniotomy at every time point in every patient is unethical. For this reason, no histological information was available in our study.

For each subject, an expert radiologist annotated the ROI of the tumor. The PK parameters were calculated per voxel, and regions with parameters out of bounds or with insufficient goodness of fitness (R2) were excluded. The focus of our work was to define which imaging biomarkers could better predict outcome before definite diagnosis in the two groups of patients with only one imaging study available.

Statistical analysis

An exploratory histogram analysis, using the derived PK parameter maps from the ROIs, was performed, yielding several metrics including the minimum, maximum, mean, standard deviation, median, skewness, kurtosis, variance, as well as the 5%, 30%, 70%, and 95% percentiles for each parameter. A Kruskal–Wallis test28 was applied to a total of 588 parameters (49 PK parameters multiplied by 12 histogram metrics) to investigate the differences between responders and nonresponders in terms of the histogram metrics. The estimated P-values were declared statistically significant only at the 1% significance level. Because of the small sample size, the estimation of the predictive power of each parameter relied on the leave-one-out cross validation (LOO-CV). For every single parameter in the dataset, a LOO-CV took place in which each sample acted as a validation set and the remaining samples were used for training. The aggregated predictive probabilities of each parameter over the entire dataset composed a new set for further analysis. The receiver operator characteristic (ROC) curves were then computed and the corresponding areas under the receiving operating curves (AUCs) were calculated to assess the performance of each parameter in predicting the responsiveness of the treatment accurately. Sensitivity, specificity, negative and positive predictive values, accuracy, the confusion matrix, and the optimal cutoff value of all ROC curves were calculated. The confusion matrix reports the number of true positive (TP), true negative (TN), false positive (FP), and false negative (FN).

All possible linear regression model combinations were generated from the parameters having a P-value <1%, and a model selection framework was intended to select the “best” model that differentiated the two groups. This framework relied on a wrapper approach equipped with parameter(s) selection criteria based on the small-sample corrected Akaike information criterion (AICc).29 AICc was selected instead of the basic AIC30 because the former is strongly recommended in case of having a small cohort for analysis.31 Using wrapper techniques, all possible univariate and multivariate models were first generated, fitted by a regression analysis function, and, based on their AICc value, the optimal model (ie, yielding the lowest AICc) was finally selected.

For the purposes of the analysis, in-house software was written in R32 to fit a linear model of the form “y ~ x1 + x2,…, + xn”, where y is the dependent variable (ie, responders vs nonresponders), and x1, x2,…,x are the explanatory variables or the parameters derived from the histogram analysis otherwise. The packages “pROC”33 and “glmulti”34 were used for computing the ROC curves, their relative statistical measurements, and the optimal model, respectively. For illustrative purposes, a graphical summary of the analysis was also obtained, containing information about the AICc profile of each model and the estimated importance of each parameter computed as the sum of the relative evidence weights of all models in which the specific parameter appears.

Alternatively, a Lasso regularized generalized linear model, using R package “glmnet”,35 was applied to the entire dataset comprised of 588 parameters. A LOO-CV was performed first to calculate the optimum tuning parameter lambda (λ) referring to the lowest CV error. The optimum lambda was used for fitting the model to the data and the resulting fitted model was then added to R function “predict”, to return the coefficients of the best model (nonzero coefficients).

Results

Parameters with P-values <1% are summarized in Table 3. The name of each parameter is composed by four parts: 1) the applied PK model (ie, GCTT or ETM as described previously), 2) the computed PK parameter (ie, k), 3) the examined tumor image subregion according to our method presented earlier (ie, 0, 1, 2, 3, and 4 for “All region”, “Full Homogeneous”, “Homogeneous”, “Heterogeneous”, and “Full Heterogeneous”, respectively), and 4) the estimated histogram metric (ie, 30 for the 30th percentile, Std for the standard deviation, etc). Table 4 shows the ROC analysis performed in the most statistically signifi-cant (by the Kruskal–Wallis test) parameters. All parameters achieved complete separation of the two groups in Table 2.

Table 3

Statistically significant parameters at the 1% significance level of the Kruskal–Wallis test.

PARAMETERS	P1	P2	P3	P4	P5	P6	P7	P8	P9	P10	P11	P-VALUE
GCTT_E_3_70	0.3511	0.3987	0.3739	0.3240	0.4140	0.5045	0.4023	0.4319	0.4172	0.4508	0.4507	<0.001
ETM_Ktrans_0_70	0.0694	0.1080	0.0842	0.0897	0.1226	0.1557	0.1232	0.1198	0.1146	0.1157	0.1169	<0.001
GCTT_Kep_4_95	0.5459	0.5284	0.5194	0.5068	0.6092	1.083	0.5582	0.8082	0.7289	0.8537	0.8836	<0.001
GCTT_PS_0_Variance	0.0008	0.0027	0.0022	0.0013	0.0033	0.0079	0.0030	0.0041	0.0028	0.0065	0.0035	<0.001
GCTT_PS_0_Std	0.0280	0.0523	0.0469	0.0356	0.0578	0.0890	0.0549	0.0643	0.0530	0.0803	0.0594	<0.001
GCTT_PS_4_5	0.0100	0.0086	0.0085	0.0102	0.0125	0.0113	0.0103	0.0122	0.0102	0.0124	0.0139	<0.001
GCTT_tc_0_95	20.47	21.11	21.18	19.73	22.21	25.93	26.44	25.84	25.90	21.56	26.86	<0.001
GCTT_Vp_1_30	0.0199	0.0243	0.0216	0.0249	0.0312	0.0381	0.0347	0.0255	0.0346	0.0416	0.0285	<0.001
GCTT_Vp_1_70	0.0297	0.0431	0.0468	0.0426	0.0501	0.0544	0.0621	0.0493	0.0535	0.0593	0.0510	<0.001

Table 4

ROC curve analysis.

PARAMETERS	ACCURACY	SENSITIVITY	SPECIFICITY	PPV	NPV	TN	TP	FN	FP	ROC CURVE CUTOFF
GCTT_E_3_70	90.9%	100%	85.7%	80%	100%	6	4	0	1	0.4005
ETM_Ktrans_0_70	90.9%	100%	85.7%	80%	100%	6	4	0	1	0.1113
GCTT_Kep_4_95	90.9%	100%	85.7%	80%	100%	6	4	0	1	0.5520
GCTT_PS_0_Variance	90.9%	100%	85.7%	80%	100%	6	4	0	1	0.0028
GCTT_PS_0_Std	90.9%	100%	85.7%	80%	100%	6	4	0	1	0.0527
GCTT_PS_4_5	100%	100%	100%	100%	100%	7	4	0	0	0.0102
GCTT_tc_0_95	100%	100%	100%	100%	100%	7	4	0	0	21.37
GCTT_Vp_1_30	90.9%	100%	85.7%	80%	100%	6	4	0	1	0.0252
GCTT_Vp_1_70	100%	100%	100%	100%	100%	7	4	0	0	0.0480

Abbreviations: PPV, positive predictive value; NPV, negative predictive value; TN, true negative; TP, true positive; FP, false positive; FN, false negative.

Following an extensive screening of all possible model combinations in the identification of the optimal model for predicting the therapeutic outcome, the information criterion profile of all models is graphically depicted in Figure 4. According to this profile, the model with the lowest AICc is composed of the parameters “GCTT_tc_0_95”, “GCTT_PS_4_5”, and “GCTT_Vp_1_30” with AICc = 7.4395. To facilitate comparison among the different parameters through the wrapper model selection, Figure 5 highlights the estimated importance of each parameter. The nonzero coefficients from the Lasso model, related to the parameters from the best fitted model, resulted in parameters “GCTT_tc_0_95”, “GCTT_PS_4_5”, “GCTT_Vp_1_30” (same parameters from the selected model using AICc), “GCTT_Vp_1_70” (fourth top-ranked parameter as depicted in Fig. 5), “GCTT_ Ktrans_3_Std”, and “GCTT_Ktrans_4_5”. The parameters “GCTT_Ktrans_3_Std” and “GCTT_Ktrans_4_5” were rejected at the preprocessing phase of the analysis using AICc because the corresponding P-values were higher than 1%.

Figure 4

Information criterion profile of all possible model combinations through the wrapper evaluation.

Figure 5

Model-averaged importance of the statistically significant parameters.

According to Figure 5, the three top-ranked parameters are the 95th percentile of the transit times in the whole region of tumor, the 5th percentile of the permeability surface area product in the “Full Heterogeneous” subregion, and 30th percentile of the plasma volume in the “Full Homogeneous” subregion, which are also the same three explanatory variables of the optimal multivariate linear regression model. Capillary transit times are indicative of hypoxia and provide important information about tumor pathophysiology. Therefore, this biomarker can be potentially related with oxygenation and patient outcome. Permeability surface area product in the “Full Heterogeneous” subregion is also a potentially meaningful biomarker since increased permeability can be attributed to a highly heterogeneous vascular bed, one of the hallmarks of malignant tumors that results in contrast leakage within the tumor tissue. It is also important to mention that ETM provided only one statistically significant parameter, the 70th percentile of k, in the whole tumor area and is related to the wash-in rate of the CA in the tissue which is the most indicative characteristic of a disrupted BBB. As can be noted in Figure 5, highly ranked parameters are also the 30th and 70th percentile of the blood plasma volume in the “Full Homogeneous” subregion. The high accuracy of this parameter in the classification might be due to the expected alteration in the vasculature of a responder, but there is no published data yet to support this. Concerning the rest of the parameters shown in Figure 5, they were found to be less significant in our study.

The statistical analysis was extended in the context of identifying any potential discrimination between subpopulations “Immune Reaction” and “Response” within the “Res ponse” group. To this end, the statistical analysis framework using AICc as the criterion for selecting the parameters with the optimum discriminatory power was applied to a cohort of four samples/patients (two patients in each group). Despite the fact that more than 50 from the 588 parameters achieved accuracy of 100%, their estimated P-values were declared as statistically insignificant (P-value >0.12). Consequently, all parameters failing to meet the criterion of the P-value threshold at the preprocessing phase (<1%) were excluded from the linear model combination process. Also, a Lasso regression model was applied to the same cohort of 588 parameters, using LOO-CV for tuning the model and calculating its nonzero coefficients. All the returned coefficients had a zero value; thus the analysis failed to define a model for predicting accurately the treatment response outcome of the examined cohort. In conclusion, the analysis showed that both models have no predictive power in classifying correctly the subpopulations “response” and “immune reaction” for the specific dataset that was used in the analysis.

Discussion

Biomarkers extracted from DCE-MRI have been used in several studies for the prediction of therapeutic outcome to chemoradiation36,37 and therapies that target vasculature1,6,38–40 in several tumor types,1,41–43 as well as for the grading of gliomas through histogram analysis.44 All these studies indicate that the relative changes of the DCE biomarkers have the potential to predict or assess patient outcome. However, pertinent studies are needed to define the actual meaning of these parameters and to confirm their robustness in describing tumor vascular heterogeneity.

In our work, statistical analysis was performed on the results derived from two PK models, the well-established ETM and the more recent GCTT; the latter has shown initial promising results in the literature.45,46 The statistical analysis was able to differentiate responders from nonresponders, although the cohort of patients was relatively small (N = 11) and all subjects were diagnosed with GBM. The statistical analysis indicates that the GCCT model outperforms ETM. This is in line with results from previous research,23 where authors argue that GCTT biomarkers are expected to be more robust. The previous assumption is reliable, given that the IRF of the model is able to be adapted to the specific tissue (per voxel) based on the heterogeneity of the transit times by adjusting the α−1 parameter. Moreover, results derived from three different types of brain tumors23 (GBM, pleomorphic xanthoastrocytoma, and anaplastic meningioma) have shown remarkable discrepancies in DCE biomarker values when different PK models were used. This inconsistency could be attributed, according to previous work,22 to the fact that TM is accurate only in weakly vascularized tissues while ETM is also accurate in highly perfused tissues, but considering that tumor vasculature is highly heterogeneous, a model such as TM or ETM that assumes a single capillary transit time is expected to diverge from reality in certain cases.

In essence, the GCCT model generalizes well-established models that have been extensively used so far and offers a new variable (α−1) that describes vascular heterogeneity. In this work we proposed an additional preprocessing step by exploiting the α−1 parameter. The tumor region of interest was segmented into four sub-regions based on the aforementioned parameter. Although these bounds were determined empirically, after several experiments we found that the classification outcome was not affected by small perturbations in the proposed boundary values for each subregion. However, when only three subregions were used instead of four, the classification failed to give satisfactory results. Finally, considering the limited number of patients and the restriction to only GBM tumor, the proposed classification based on α−1 boundaries should be verified in a more extensive patient group. A thorough sensitivity analysis should be performed as well in order to define consistent subregion boundaries that could be used prospectively in clinical decision support.

In our study, the GCTT model’s PK parameters performed better than Tofts’, while the segmentation of the tumor ROIs based on vascular heterogeneity further enhanced the discriminatory power of the GCTT model. The results indicate that the GCTT model may be more efficient in characterizing the chaotic vascular structure of tumor and subsequent pathophysiological characteristics such as the delayed extraction of the CA in the EES. Also, in search for the best model according to the optimal multivariate linear regression presented, the GCCT PK parameters outperformed the Tofts ones. k computed from ETM performed well in our study, and in particular the 70% percentile. Figure 6 shows the histogram profiles for k for all the patients. The ROC results showed, however, that the overall best PK parameter is the extraction fraction E (70% percentile) of the GCCT model but only when computed in the “Heterogeneous” subregion (in tables of results referred to as region 3). This is further highlighted in Figure 7, where it can be observed that the histograms of this parameter between responders and nonresponders become morphologically more separable in the “Heterogeneous” subregion than in the whole tumor ROI as annotated in the MRI data. If this result is confirmed in future studies, it has the potential to enhance the robustness of PK imaging biomarkers from DCE-MRI and widen their clinical adoption for aiding the therapy monitoring process. However, being aware of the limited cohort group of this study and the restriction to only GBM tumor, we suggest more extensive studies in a broader range of patients and tumor types in order to establish which model is better for the early prediction of response in cancer patients.

Figure 6

k PK parameter computed in the whole tumor region of interest also offers a good separability while the 70% percentile was the best metric.

Figure 7

Histograms for the extraction fraction GCCT model parameter E for the whole tumor ROI (left) and the “Heterogeneous” subregion (right). It can be observed that there is a morphological enhancement in terms of separability between responders and nonresponders in the “Heterogeneous” subregion in comparison to the whole tumor ROI as annotated in the MRI data.

Table A1

Abbreviation and symbols.

ADC	Apparent diffusion coefficient
AIF	Arterial input function
AIC	Akaike information criterion
AICc	Small-sample corrected akaike information criterion
ATH	Adiabatic tissue homogeneity
AUCs	Areas under the roc-curve
BBB	Blood brain barrier
CA	Contrast agent
DCE-MRI	Dynamic contrast-enhanced magnetic resonance imaging
EES	Extracellular extravascular space
ETM	Extended tofts model
FN	False negative
FP	False positive
GBM	Glioblastoma multiforme
GCTT	Gamma capillary transit time
Gd	Gadolinium
IRF	Impulse response function
MRI	Magnetic resonance imaging
MTT	Mean transit time
PK	Pharmacokinetic
rCBF	Regional cerebral blood flow
rCBV	Regional cerebral blood volume
ROC	Receiver operator characteristics
ROI	Region of interest
SNR	Signal to noise ratio
SPGR	Spoiled gradient echo
TM	Tofts model
TN	True negative
TP	True positive
VOI	Volume of interest
2CX	Two compartment exchange

Table A2

PK parameters derived from DCE-MRI data.6

PARAMETER	DEFINITION	UNIT
Models parameters
F	Blood flow (normalized by tissue volume)	mL/mL/min−1
E	Extraction fraction	none
ktrans	Volume transfer constant between plasma and EES	min−1
kep	Rate constant between EES and plasma	min−1
τ	Scale parameter of gamma distribution	sec
α−1	Distribution of capillary transit times	none
vp	Blood plasma volume per unit volume of tissue	none
ve	EES volume per unit volume of tissue	none
Additional parameters
PS	Permeability surface area product per unit mass of tissue	mL/mL/min−1
tc	Capillary transit time	sec
Biophysical parameters
T10	Pre-contrast longitudinal relaxation time	sec
T1(t)	Post-contrast longitudinal relaxation time	sec
S0	Equilibrium magnetization	none
Ct(t)	Concentration in tissue	mM
Ca(t)	Concentration in artery blood (AIF)	mM
IRFv(t)	Vascular impulse response function	none
IRFp(t)	Parenchyma impulse response function	none
Hct	Hematocrit	none
Imaging parameters
TS	Temporal sampling interval	sec
TR	Repetition time	sec
S1	Signal intensity	none
α	Flip Angle (FA)	degrees
αf	Variable Flip Angle (VFA)	degrees
r1	Relaxation constant of the contrast agent	sec−1·mM−1