Image-based personalization of computational models for predicting response of high-grade glioma to chemoradiation

High-grade gliomas are an aggressive and invasive malignancy which are susceptible to treatment resistance due to heterogeneity in intratumoral properties such as cell proliferation and density and perfusion. Non-invasive imaging approaches can measure these properties, which can then be used to calibrate patient-specific mathematical models of tumor growth and response. We employed multiparametric magnetic resonance imaging (MRI) to identify tumor extent (via contrast-enhanced T1-weighted, and T2-FLAIR) and capture intratumoral heterogeneity in cell density (via diffusion-weighted imaging) to calibrate a family of mathematical models of chemoradiation response in nine patients with unresected or partially resected disease. The calibrated model parameters were used to forecast spatially-mapped individual tumor response at future imaging visits. We then employed the Akaike information criteria to select the most parsimonious member from the family, a novel two-species model describing the enhancing and non-enhancing components of the tumor. Using this model, we achieved low error in predictions of the enhancing volume (median: − 2.5%, interquartile range: 10.0%) and a strong correlation in total cell count (Kendall correlation coefficient 0.79) at 3-months post-treatment. These preliminary results demonstrate the plausibility of using multiparametric MRI data to inform spatially-informative, biologically-based predictive models of tumor response in the setting of clinical high-grade gliomas.

registration algorithm (imregister in MATLAB (Mathworks, Natick, MA)) was used to provide intra-and inter-visit registration. For intra-visit registration, each image and parameter map were registered to the T2-FLAIR image collected at that visit. The intra-visit registration consisted of resampling of images to match the image resolution of the T2-FLAIR image. For the inter-visit registration, each image and parameter map acquired at post-baseline image sessions were registered to the T2-FLAIR image acquired at baseline. Panel A in Figure 1 shows an example of the registration process for a given patient.
For each patient, the gross tumor volume (GTV) was defined as the enhancing tumor on the postcontrast T1-weighted images and clinical tumor volume (CTV) defined as non-enhancing, T2-hyperintense region on the T2-FLAIR images were segmented using a semi-automated approach using thresholding methods in combination with manual adjustments by a radiation oncologist and secondary quality review by a second senior radiation oncologist. A k-means clustering of signal intensity was used to segment the white matter, gray matter, and cerebrospinal fluid from T2-FLAIR images [1]. The brain-skull interface was manually segmented from the T2-FLAIR image. Panel B in Figure 1 shows an example of this segmentation process.
The ADC calculated from DWI data was used to estimate the tumor cell volume fraction at each imaging visit. ADC was estimated voxel-wise from an echo planar imaging based DWI sequence using Eq. (1): , (1) where b1 and b2 are the b-values of 0 and 1000 s/mm 2 , respectively, and and are the signal intensities corresponding to b1 and b2, respectively. The tumor volume fraction at 3D position and time t , , was then estimated on a voxel-specific basis as described in [2][3][4][5], using the ADC of free water (ADCw) [6], the minimum ADC measured (ADCmin), and Eq. (2): We have used this approximation previously to provide non-invasive estimates of tumor cellularity [3][4][5]7,8]; however, we note that this approximation is a simplification of all the biological aspects that contribute to changes in ADC. This point is discussed further in [3,9]. Within the GTV, we assumed that the primary cellular contribution is from tumor cells, therefore was calculated using Eq. (2).
However, within the CTV the cell density or relationship to imaging features is less clear, thus we used a fixed value of 0.16 [10] everywhere within that region. An alternative approach for assigning celluarity in the non-enhancing regions (used in [10,11] and elsewhere) is to assume a spatially varying value of cellularity decreasing from the value observed at the interface of the enhancing region to a fixed value at the periphery of the non-enhancing region. For the two-species model, the tumor volume fraction within the enhancing or GTV region was calculated using Eq. (2) and set to zero outside the GTV, while in the non-enhancing or CTV region, it was set to a fixed value of 0.16.

S.2. Numerical implementation
The spatial-temporal evolution of was determined using a 3D finite difference approximation implemented in MATLAB R2019b (Mathworks, Natick, MA) using a fully explicit in time differentiation and central difference spatial differentiation. The numerical time step was selected to maintain numerical stability based on the grid spacing and diffusion coefficient values. No flux (Neumann) boundary conditions were used for at the skull boundary. The boundary condition for was assumed to be zero displacement in the normal direction, while it was assumed that the tissue in the tangential directions was free to move (i.e., slip condition). Identical numerical approaches were also used for and . For a complete description of the numerical implementation of these techinques, the interested reader is referred to [12].

S.3 Model parameter calibration
A total of 40 models ( Figure 1C) were developed from two base models (i.e., one-species or twospecies), 10 therapy coupling combinations (i.e., C1 to C4), and two proliferation parameterization approaches (described below). First, we assumed kp,T and kp,E were uniform throughout the domain.
Second, we assumed that kp,T and kp,E varied as a field within the tumor regions of interest. When the proliferation rates were calibrated as a field, we only calibrated for a subset of points within the tumor and then interpolated elsewhere. For example, for a given 3 × 3 voxel sub-region within the tumor, the parameter values were calibrated at the corner and center positions while the remaining four points were interpolated from the nearest calibrated values [9]. This parameter calibration approach significantly reduces the number of individual parameters and results in a spatially smooth parameter field. To reduce the number of model parameters for the two-species model, we assumed that the DN,w and DN,g are proportional to DE,w and DE,g, respectively, by a single calibrated scalar factor. Likewise, we assumed kp,N is proportional to kp,E. bNE and bEN were empirically determined and assigned to values of 4 and 1, respectively. The remaining calibrated model parameters (in Table 2) were fit as a global variable. Model parameters are bounded as shown in Supplemental Table 2 based on their physical definition, numerical stability, or within an order of magnitude of values reported in literature.

Supplemental
where Dt is the simulation time step, Dx is the image resolution in the x-direction, Dy is the image resolution in the y-direction, and Dz is the image resolution in the z-direction.
Panel D in Figure 1 depicts the calibration approach used in this study. For each patient, we considered three different calibration/prediction scenarios. For the first scenario, we calibrated each model to all of the available data to see how well the models describe that data. For the second and third scenarios, we calibrated each model to a subset of the available data (baseline and 1-month for scenario 2; baseline, 1-month, and 3-month for scenario 3) and then those calibrated parameters are used to run the model forward in time to predict the tumor response at that patient's remaining imaging visits (i.e., the 3-month and 5-month visits).
While the salient details of the calibration approach are highlighted here, a more complete description of the calibration algorithm can be found in the referenced publication [12]. Figure 2 provides a schematic of our model calibration approach. Briefly, we used the Levenberg-Marquardt [12,13] 0.25 > D stability Δt algorithm to minimize the sum of the squared errors between the measured and simulated tumor growth where and are the model estimated tumor volume fractions for a given set of parameters b for the enhancing and non-enhancing region, respectively, and and are the measured tumor volume fractions for the enhancing and non-enhancing regions, respectively. We note that in Eq. (5) the residuals between both species are normalized by their carrying capacity to account for magnitude differences in the two volume fractions.

S.4 Sensitivity of model calibration to measurement noise
We evaluated the model calibration approach using an in silico tumor (i.e., a tumor whose growth is governed by the models above with user-selected parameter values). The purpose of this in silico study was to evaluate the ability of our parameter calibration approach to accurately determine the true model parameter governing the in silico tumor's growth and response. For this in silico study we assessed the robustness of the parameter calibration in the presence of 10% measurement noise. For the model with the most parameters (i.e., the two-species model with a calibrated proliferation field), 250 noisy datasets were generated by adding random noise from a Normal distribution with a standard deviation of 10%. (In our study, the largest variation we observed for a given patient in voxel-wise values of ADC for normal appearing white matter was a mean difference of -2.9% with a standard deviation of 9.7% compared to the baseline image. Thus, the 10% value is appropriate.) Model parameters were then estimated from the noisy datasets and then compared to the known parameters used to "grow" the in silico tumor. Initial guesses for model parameters for each calibration scenario were selected randomly from a uniform distribution with lower and upper bounds shown in Supplemental Table 2 to ensure parameter guess were not near the true values. We observed that less than 5.6% error was observed when two time points were used for calibration or less than 3.33% when three or more time points were used for calibration.
Supplemental Table 3 reports the full results for this exercise.

S.5 Model selection
The Akaike Information criterion (AIC, [14]) was used to select the model that optimally balances model complexity and model-data agreement. The AIC was calculated using: , (6) where k is the number of parameters calibrated for a given model, n is the number of data points used to calibrate the model, and RSS is the residual sum squares between the measured and model estimated tumor growth. We calculated the AIC for each model over the timepoints used for model calibration. We then selected the model with the lowest average AIC across all patients as the most parsimonious model. We then calculated the Akaike weights for each model defined as: where wi is the weight for the i-th model, is equal to AICi -AICmin, and AICmin is the minimum observed AIC. The model with the lowest average AIC was the two species model, with a locally varying proliferation rate and radiation and chemotherapy both coupled to approach C2 (i.e., coupled to ER). This model will be used in all of the model calibrations and predictions reported in the results. Supplemental Table 4 reports the ten model combinations with the lowest AIC. Of these ten models, six were the twospecies model, four were the single-species model, and seven had a spatially varying proliferation rate.

S.6 Individual patient results for the prediction phase
Supplemental Figure 1 reports prediction results for the selected model to an assumed "no-growth" or static model. That is, for the no-growth model we assume that the tumor growth at the 3-month and 5-  These are the parameter values when the selected model is calibrated to the entire tumor growth time course. The model that was selected via the Akaike Information Criterion was the two species reaction-diffusion model, with a field proliferation rate, and coupling combination 8. For the two-species model the following parameters are not calibrated: kp,T,, qT, DT,w, DT,g,. *kp,E is reported as an average of the parameters over that field.