Next Article in Journal
BTLA/HVEM Axis Induces NK Cell Immunosuppression and Poor Outcome in Chronic Lymphocytic Leukemia
Next Article in Special Issue
Prediction of Human Papillomavirus (HPV) Association of Oropharyngeal Cancer (OPC) Using Radiomics: The Impact of the Variation of CT Scanner
Previous Article in Journal
A Novel Two-Lipid Signature Is a Strong and Independent Prognostic Factor in Ovarian Cancer
Previous Article in Special Issue
MR Image Changes of Normal-Appearing Brain Tissue after Radiotherapy
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Towards an Image-Informed Mathematical Model of In Vivo Response to Fractionated Radiation Therapy

by
David A. Hormuth II
1,2,*,
Angela M. Jarrett
1,2,
Tessa Davis
3 and
Thomas E. Yankeelov
1,2,3,4,5,6
1
Oden Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX 78712, USA
2
Livestrong Cancer Institutes, The University of Texas at Austin, Austin, TX 78712, USA
3
Department of Biomedical Engineering, The University of Texas at Austin, Austin, TX 78712, USA
4
Department of Diagnostic Medicine, The University of Texas at Austin, Austin, TX 78712, USA
5
Department of Oncology, The University of Texas at Austin, Austin, TX 78712, USA
6
Department of Imaging Physics, The University of Texas MD Anderson Cancer Center, Houston, TX 77030, USA
*
Author to whom correspondence should be addressed.
Cancers 2021, 13(8), 1765; https://doi.org/10.3390/cancers13081765
Submission received: 23 February 2021 / Revised: 1 April 2021 / Accepted: 4 April 2021 / Published: 7 April 2021
(This article belongs to the Special Issue Transformational Role of Medical Imaging in Oncology)

Abstract

:

Simple Summary

Using medical imaging data and computational models, we develop a modeling framework to provide personalized treatment response forecasts to fractionated radiation therapy for individual tumors. We evaluate this approach in an animal model of brain cancer and forecast changes in tumor cellularity and vasculature.

Abstract

Fractionated radiation therapy is central to the treatment of numerous malignancies, including high-grade gliomas where complete surgical resection is often impractical due to its highly invasive nature. Development of approaches to forecast response to fractionated radiation therapy may provide the ability to optimize or adapt treatment plans for radiotherapy. Towards this end, we have developed a family of 18 biologically-based mathematical models describing the response of both tumor and vasculature to fractionated radiation therapy. Importantly, these models can be personalized for individual tumors via quantitative imaging measurements. To evaluate this family of models, rats (n = 7) with U-87 glioblastomas were imaged with magnetic resonance imaging (MRI) before, during, and after treatment with fractionated radiotherapy (with doses of either 2 Gy/day or 4 Gy/day for up to 10 days). Estimates of tumor and blood volume fractions, provided by diffusion-weighted MRI and dynamic contrast-enhanced MRI, respectively, were used to calibrate tumor-specific model parameters. The Akaike Information Criterion was employed to select the most parsimonious model and determine an ensemble averaged model, and the resulting forecasts were evaluated at the global and local level. At the global level, the selected model’s forecast resulted in less than 16.2% error in tumor volume estimates. At the local (voxel) level, the median Pearson correlation coefficient across all prediction time points ranged from 0.57 to 0.87 for all animals. While the ensemble average forecast resulted in increased error (ranging from 4.0% to 1063%) in tumor volume predictions over the selected model, it increased the voxel wise correlation (by greater than 12.3%) for three of the animals. This study demonstrates the feasibility of calibrating a model of response by serial quantitative MRI data collected during fractionated radiotherapy to predict response at the conclusion of treatment.

1. Introduction

With over a century of development, radiotherapy (RT) remains a principal component of oncological care for many disease sites including brain cancer. In the particular case of high-grade gliomas, RT is administered alongside chemotherapy to treat residual and unresected disease [1]. The delivery of highly conformal RT doses has been enabled by the development of intensity-modulated RT and, more recently, magnetic resonance imaging (MRI) guided linear accelerators [2]. These conformal RT techniques can deliver the prescribed dose to the tumor while significantly reducing the dose delivered to healthy appearing tissue [3]. However, the full potential of RT techniques have not yet been fully realized as current RT treatment plans are based predominately on anatomical or structural imaging measures that report on the extent of the disease, rather than the underlying biology which may provide more fundamental information on treatment response [4]. Several quantitative imaging techniques have emerged as potential imaging biomarkers [4] to adapt or boost RT plans [5]; however, they still only provide a temporal snapshot of the tumor biology. Imaged-based mathematical modeling techniques [6] may provide a means to forecast tumor dynamics in response to RT, and thereby enable the personalization of RT plans for each patient’s unique tumor biology.
Cell survival following RT is commonly described by the linear quadratic (LQ) model, which reports the fraction of cells that will survive and proliferate to form a colony of cells given a specific radiation dose and cell-specific radiosensitivity parameters [7]. One strength of the LQ model is its simplicity for evaluating the efficacy of different RT regimens for tumor and healthy appearing tissue, thus enabling the selection of fractionated RT doses that have the most therapeutic benefit. One shortcoming, however, is that it does not capture the temporal dynamics of tumor response. That is, following a single fraction of RT the LQ model cannot estimate when a tumor will shrink or when tumor cells will die. Response to RT is not binary (i.e., cells do not either die or emerge unscathed from exposure to radiation), rather there is a spectrum of response ranging from cells that receive sufficient damage to undergo apoptosis, to cells that continue proliferating for a few generations before becoming senescent, or cells that proliferate indefinitely due to successful DNA repair [8]. In recent years, we [9] and others [10,11,12,13] have attempted to connect models capturing the dynamics of tumor growth with the efficacy of RT. Beyond characterizing response to RT, these mechanism-based and machine learning techniques can also be used to predict response [9,14], identify optimal treatment strategies [15,16,17,18], and evaluate alternative schedules based off different growth assumptions [19]. In the works of Rockne et al. [13] and Prokopiou et al. [12], the effects of RT were assumed to result in the immediate loss of tumor cells (or reduction in tumor volume) as estimated from the LQ model. In a leave-one-out cross validation, Rockne et al. observed that the predicted post-RT tumor volume was well within the inter-observer tumor volume uncertainty [13] demonstrating a potential approach to evaluating the efficacy of RT regimens. Prokopiou et al. [12] proposed the use of a proliferation saturation index to identify patients that might benefit from non-standard fractionation regimens; this modeling framework was later shown to provide strong predictive accuracy of future response [20]. Alternatively, non-instantaneous cell death following RT was considered by Brüningk et al. [11] and Hormuth et al. [9] in their models of response in the in vitro and in vivo pre-clinical setting, respectively. Both Brüningk et al. [11] and Hormuth et al. [9] observed that models with an instantaneous cell death greatly increased the error in model fits and/or predictions compared to models with a delayed death or altered proliferation, respectively. However, it should be noted that models of instantaneous death may be more appropriate for longer time scales (i.e., days to weeks) when sufficient time has passed for the tumor to shrink. While these efforts have shown promising results capturing the temporal dynamics of response to RT, techniques capable of also predicting the spatial response to fractionated RT are still lacking. Towards this goal, we have extended our image-informed modeling approach [6,9,21] to consider the response of tumor and vasculature to fractionated RT in a pre-clinical model of high-grade glioma.
Our previous modeling studies considered single, large dose, RT [9], which does not mirror the standard-of-care for high-grade glioma patients which consists of 2 Gy per weekday for 6 weeks [22]. To address this limitation, we revised our experimental approach to delivery 2 Gy or 4 Gy per weekday for up to 2 weeks. Additionally, animals were imaged with quantitative magnetic resonance imaging (MRI) techniques sensitive to tumor cellularity and tumor vasculature [23,24] before, during, and following the completion of RT. By using quantitative MRI, we are able to non-invasively observe and quantify tumor physiology in 3D throughout the study, which is then used for calibration of animal specific tumor growth and response parameters as well as validation of predicted tumor response to treatment. From the observed dynamics of tumor growth and response, we developed three models of response to RT. The first model assumed immediate cell death and removal (similar to the approaches used by [12,13]). The second model assumed immediate cell death and a reduction of tumor proliferation rate (due to cell cycle arrest), while the third model assumed a reduction of tumor proliferation rate and an increase in tumor cell death rate (due to delayed mitotic cell death). The three models of response to RT were incorporated in a family of 18 models that included three approaches to spatially-vary the efficacy of RT and two approaches to parameterizing the tumor cell proliferation rate. In this preliminary study, we calibrated the entire family of 18 models to each individual animal and evaluated the descriptive (i.e., how well the model fits the data) and predictive (i.e., how well the model predicts future response) accuracy of these models. We then used the Akaike Information Criterion (AIC) to select the best or most parsimonious model and calculate model weights. We then quantified the descriptive and predictive error for both the selected model and the ensemble averaged model.

2. Materials and Methods

2.1. Theory

2.1.1. Mathematical Description of Tumor and Vasculature Growth

Tumor and vasculature growth are described using a coupled set of 3D partial differential equations built upon the reaction-diffusion model that has been extensively studied at the pre-clinical and clinical settings in glioma [9,25,26,27] and other tumors [28,29]. In our previous efforts, we have applied this model in the setting of untreated tumor growth [21] and single fraction radiotherapy [9]. Here, it is used as the base model upon which we build a set of models of response to fractionated radiotherapy (described in the Section 2.2). We now present the salient features of the model which captures tumor and vasculature growth in the absence of treatment and list the model parameters and variables in Table 1. (The extension of the model to account for the effects of radiation are described below in Section 2.1.2) While a more detailed description of the development and model assumptions are found elsewhere [21], Equation (1) is the main operational equation and describes the spatial and temporal evolution of the tumor volume fraction:
ϕ T ( x ¯ , t ) t = · ( D T ( x ¯ , t ) · ( ( ( 1 ϕ V ( x ¯ , t ) θ T , V ( x ¯ , t ) ) ϕ T ( x ¯ , t ) θ T , V ( x ¯ , t ) ) + ( ϕ T ( x ¯ , t ) θ T , V ( x ¯ , t ) ϕ V ( x ¯ , t ) θ T , V ( x ¯ , t ) ) ) ) D i f f u s i o n + k p , T ϕ T ( x ¯ , t ) ( 1 ϕ T ( x ¯ , t ) / θ T ( x ¯ , t ) ) L o g i s t i c   G r o w t h
where ϕ T ( x ¯ , t ) is the tumor cell volume fraction at three-dimensional position x ¯ and time t, D T ( x ¯ , t ) is the tumor cell diffusion coefficient coupled to local mechanical properties, ϕ V ( x ¯ , t ) is the blood volume fraction, θ T , V ( x ¯ , t ) is the summation of tumor and the blood volume fraction carrying capacities, kp,T is the tumor cell proliferation rate, and θ T ( x ¯ , t ) is the tumor cell carrying capacity (i.e., the maximum packing fraction that a voxel can functionally support). We assume tumor cell diffusion (the first term on the right hand side of Equation (1)) is influenced both by the space occupied by tumor associated vasculature (i.e., ϕ V terms) as well as local mechanical stress that alter D T ( x ¯ , t ) according to:
D T ( x ¯ , t ) = D T , 0 exp ( λ 1 · σ v m ( x ¯ , t ) )
In Equation (2), DT,0 is the unrestricted (or maximal) tumor cell diffusion coefficient in the absence of stress, λ 1 is the stress-tumor cell diffusion coupling constant, and σ v m ( x ¯ , t ) is the von Mises stress. The von Mises stress is used as it reflects the total stress experienced for a given section of tissue, and is determined by solving for tissue displacement, u , using the linear elastic, isotropic equilibrium equation defined as:
· G u + G 1 2 v ( · u ) λ 2 ϕ T ( x ¯ , t ) = 0
where G is the shear modulus, ν is the Poisson’s ratio, and λ 2 is the second coupling constant. Literature values [30] are used to assign tissue specific G and ν.
Our previous modeling study [21] identified that the optimal coupling of tumor vasculature to tumor cell growth was via the carrying capacity θ T ( x ¯ , t ) . Thus, θ T ( x ¯ , t ) is evolved spatially and temporally in response to changes in the local tumor vasculature according to:
θ T ( x ¯ , t ) = { θ max ϕ V ( x ¯ , t ) ϕ V , t h r e s h θ min + ϕ V ( x ¯ , t ) ( θ max θ min ϕ V , t h r e s h ) ϕ V ( x ¯ , t ) < ϕ V , t h r e s h
where the range of expected carrying capacities is from θ min to θ max , and ϕ V , t h r e s h represents a critical value for tumor vasculature that would begin to change the number of cells a voxel can support. The parameters θ min and θ max are assigned as the lowest and highest volume fractions, respectively, observed during the image visits used for calibration. In a similar fashion, we describe the spatial-temporal evolution of tumor vasculature using:
ϕ V ( x ¯ , t ) t = · ( D V ( x ¯ , t ) · ( ( ( 1 ϕ T ( x ¯ , t ) θ T , V ( x ¯ , t ) ) ϕ V ( x ¯ , t ) θ T , V ( x ¯ , t ) ) + ( ϕ V ( x ¯ , t ) θ T , V ( x ¯ , t ) ϕ T ( x ¯ , t ) θ T , V ( x ¯ , t ) ) ) ) D i f f u s i o n + k p , V ϕ V ( x ¯ , t ) ( 1 ϕ V ( x ¯ , t ) / θ V ) d L o g i s t i c   G r o w t h / A n g i o g e n e s i s k d , V ϕ V ( x ¯ , t ) ( 1 d ) D e a t h
where D V ( x ¯ , t ) is the mechanically coupled vascular diffusion coefficient (described in a similar fashion as Equation (2)), kp,V is the tumor vasculature growth rate, θ V is the blood volume fraction carrying capacity (assigned as the maximum observed blood volume fraction), d is a normalized parameter describing the distance to the periphery of the tumor (1 for voxels at the periphery, 0 for voxels furthest from the periphery), and kd,V is the vascular death rate. We note that occupancy by tumor and blood volume fractions are considered in the cross-diffusion terms, but not explicitly in the logistic growth terms in Equations (1) and (5). There are two main reasons for this choice. First, the blood and tumor volume fractions have different maximum carrying capacities. By not considering the different maximum carrying capacities it is possible to simulate voxels that are completely vascularized (which is not observed in these tumors) and would, therefore, inhibit further growth of tumor cells within that voxel. Second, we assume the blood volume fraction to dynamically update the carrying capacity of the tumor cells (via Equation (4)). In general, the carrying capacities should be interpreted as the limitations in the number of cells that a voxel can support based on the available resources (e.g., oxygen, glucose) and space. While a decrease in vascularization would allow more space for tumor cells to grow, there would be insufficient resources to support that growth thereby resulting in a decrease of the carrying capacity.
Equations (1)–(5) provide a coupled set of partial different equations that describe the spatio-temporal development of tumor cells and vasculature in the absence of RT. The following sections describe how we generate a family of models that incorporate the effects of RT in different ways.

2.1.2. Mathematical Descriptions of Response to Fractionated Therapy

We developed a set of three models (RTM1 through RTM3) that are used to describe the response of tumor and vasculature to fractionated radiation therapy. The first model (RTM1) assumes there is a surviving fraction (SF) of tumor cells that remains following RT described by:
ϕ T , p o s t R T = ϕ T , p r e R T ( C i · S F )
where ϕT,post-RT is the tumor cell volume fraction immediately after RT, ϕT,pre-RT is the tumor cell volume fraction immediately before RT, SF is the surviving fraction of tumor cells (between 0 and 1), and Ci is one of three coupling approaches discussed below. RTM1 is similar to how the LQ model has been used in other modeling approaches [13]; however, here we do not assume radiosensitivity parameters (e.g., α or β) a priori and simply fit (for each animal) the effect of RT as the SF. During simulation time steps that correspond to RT treatment times, Equation (6) is evaluated prior to solving Equations (1)–(5). The post-RT tumor volume fraction, ϕT,post-RT, is then used as the current estimate of the tumor distribution used in the finite difference approximation of Equations (1)–(5). RTM2 includes the death term from RTM1, while also reducing kp,T, in Equation (1), with each additional fraction of RT:
k p , T = k p , T , 0 ( C i · S F n )
where kp,T,0 is the pre-treatment proliferation rate, and SFn represents the fraction of tumor cells that are able to proliferate following n fractions of RT. RTM3 removes the assumption of instantaneous death (i.e., Equation (6)), and instead assumes tumor cells die at a rate kd,T:
k d , T = k d , T , 0 ( 1 C i · S F n )
where kd,T,0 is the death rate due to RT, and ( 1 C i · S F n ) is the fraction of dying cells after n fractions of RT. RTM3 amends Equation (1), by adding an exponential death term shown below:
ϕ T ( x ¯ , t ) t = · ( D T ( x ¯ , t ) · ( ( ( 1 ϕ V ( x ¯ , t ) θ T , V ( x ¯ , t ) ) ϕ T ( x ¯ , t ) θ T , V ( x ¯ , t ) ) + ( ϕ T ( x ¯ , t ) θ T , V ( x ¯ , t ) ϕ V ( x ¯ , t ) θ T , V ( x ¯ , t ) ) ) ) D i f f u s i o n + k p , T ϕ T ( x ¯ , t ) ( 1 ϕ T ( x ¯ , t ) / θ T ( x ¯ , t ) ) L o g i s t i c   G r o w t h k d , T ϕ T ( x ¯ , t ) D e a t h   d u e   t o   R T    
Identical formulations of Equations (6)–(8) are used to model the effect of RT on the vasculature. In addition to these three models (i.e., Equations (6)–(8)) of radiation induced cell death, we also developed three ways to spatially vary the efficacy of RT through the coupling coefficient Ci [9]. The first coupling approach (C1) maximizes the SF as ϕT approaches θT due to an assumed slower tumor cell proliferation (thereby making the cells less susceptible to RT):
C 1 ( x ¯ , t ) = ϕ T ( x ¯ , t ) / θ T ( x ¯ , t )
The second coupling approach relates the efficacy of RT to ϕV via:
C 2 = exp ( α 1 ( ϕ V ( x ¯ , t ) / θ V ) )
where α1 is a coupling constant. Thus, regions with high ϕV experience an increased treatment efficacy (due to assumed increased oxygenation). For the third coupling approach (C3), we assume C3 is equal to 1 everywhere. Combining the three models of response to RT (Equations (6)–(8)) and the three approaches to spatially vary response (i.e., C1C3) we have a total of nine models of response to radiation therapy.

2.2. Experimental Methods

2.2.1. Animal Model and Radiation Therapy Protocol

All experimental procedures were approved by our Institutional Animal Care and Use Committee. Seven female athymic Hsd:RH-Foxn1RNU nude rats (weighing from 186 to 220 g) were purchased from Envigo (Indianapolis, IN, USA). The human glioma cell line U-87 MG was purchased from ATCC (ATCC HTB-14, Manassas, VA, USA) and were grown as per the packaging. U87-MG cells were cultured to 80% confluence (1–2 weeks) in Eagle’s Minimal Essential Media (ATCC 30-2003) supplemented with 10% FBS (A3160502, Thermo Fisher, Waltham, MA, USA), 100U/mL penicillin-streptomycin (15-140-122, Fisher Scientific, Houston, TX, USA), 500 ng/mL amphotericin B (BS721, BioBasic, Amherst, NY, USA), 500 ng/mL plasmocin mycoplasma prophylactic (ant-mpp, Invivogen, San Diego, CA, USA). Cells were then trypsinized, washed 2×, and resuspended in 1 × PBS at a concentration of 8 × 107/mL. For tumor cell injections [31], the animals were anesthetized (2% isoflurane in 100% oxygen) and the needle was positioned 1 mm posterior and 2 mm right lateral to the Anterior-Posterior (AP) and Medial-Lateral (ML) coordinates of bregma, respectively, then advanced 4 mm ventral (deep) to skull surface. The needle was then lowered into the skull slowly over 2 min, then allowed to remain in situ another 2 min prior to injection. Injection was carried out at a rate of 0.4 µL/min over approximately 6 min. After injection, pressure was allowed to equalize with the needle in situ for 10 min prior to retraction. A total of 2 × 105 U-87 MG glioma cells in a 2.5 µL volume. 48 h prior to their first MRI study, a permanent jugular catheter was placed in each rat.
Imaging studies typically began when the tumor volume reached 30 to 50 µm3. Three days later, animals were randomly assigned to either the 2 Gy or 4 Gy per fraction treatment group. Fractionated RT was delivered five days per week for two consecutive weeks. Whole brain RT was delivered using the MultiRad350 (350 kVp/4.7 mA, Precision X-ray Irradiation, North Branford, CT, USA). During the irradiation procedure, the animals were anesthetized with a mixture of 2% isoflurane in 100% oxygen, and then positioned within the irradiation chamber in a prone position with the animal’s head placed at the center of the irradiation platform. Lead blocks were placed to shield the animal’s torso. RT was delivered at a dose rate of 1.0 Gy/min in the dorsal to ventral direction, and total dose delivered was controlled by setting the exposure time to either 2 min or 4 min. The RT beam was filtered using a Sn-Al-Cu filter (Precision X-ray Irradiation). The dose rate and planned total dose for both the 2 min and 4 min exposure was verified daily using the MultiRad 350’s integrated dosimeter prior to animal irradiation. Imaging continued up to three times per week during and after the completion of RT. We do note, however, that the imaging and RT experimental time line varied between animals due to animal health concerns (e.g., excessive weight loss or dehydration, tumor growth impacting mobility) that required them to be euthanized prior to the end of the study. Table 2 reports the exact imaging and RT time line for each animal.

2.2.2. Imaging Procedure

Multiparametric MRI was acquired using a 7.0T horizontal-bore magnet (Bruker Biospec, Billerica, MA, USA) with a 60 mm diameter volume coil over a 32 × 32 × 16 mm3 field of view. Multiparametric MRI consisted of inversion recovery data to construct a T1-map, T2-weighted MRI for anatomical images, diffusion-weighted (DW-) MRI data to compute the apparent diffusion coefficient (ADC), and dynamic contrast-enhanced (DCE-) MRI data to compute the blood volume fraction. All images were acquired with a 128 × 128 matrix and 16 slices.
Data for the pre-contrast T1 map were acquired using a segmented FLASH (segFLASH) inversion recovery sequence with: α = 15°, TE = 3.2 ms, segment size = 12. The longitudinal relaxation curve was sampled at 30 inversion times (TI) ranging from 125 to 3025 with a 100 ms spacing. Voxel-wise T1 values were estimated by fitting Equation (12) to each voxel’s relaxation curve [32]:
S ( T I ) = A B exp ( T I / T 1 * ) A = S 0 T 1 * / T 1 B = S 0 ( 1 + T 1 * / T 1 ) T 1 = T 1 * ( B A 1 )
where T 1 * is the effective T1, and S0 is the inherent signal intensity. A, B and T 1 * were estimated using a non-linear least square optimization (lsqcurvefit in MATLAB, Mathworks, Natick, MA, USA) and then used to calculate T1.
T2-weighted MRI data was acquired using a fast spin-echo or rapid imaging with refocused echoes (RARE) sequence with the following pulse sequence parameters: TR = 3500 ms, TE = 14 ms, RARE factor of 8, and number of excitations = 10.
DW-MRI data was acquired using a pulsed gradient echo sequence with three b-values (150, 350, and 800 s/mm2) and gradients applied simultaneously along the three orthogonal directions with the following pulse sequence parameters: TR = 2500 ms, TE = 28.7 ms, number of excitations = 2, Δ = 23 ms, and δ = 3 ms. The DW-MRI data was modeled by a standard mono-exponential decay to estimate the apparent diffusion coefficient (ADC) voxel-wise within the tumor. The tumor cell volume fraction, ϕ T ( x ¯ , t ) , was then estimated directly from the ADC map using Equation (13):
ϕ T ( x ¯ , t ) = ( A D C w A D C ( x ¯ , t ) A D C w A D C min )
where ADCw is the ADC of free water at 37 °C [33], A D C ( x ¯ , t ) is the ADC value at position x ¯ and time t, and ADCmin is the minimum ADC observed within the tumor regions-of-interest (ROIs) across all animals. We note that while we have extensively used the ADC to estimate tumor cellularity [9,21,25], we acknowledge that there are other factors that may influence the measured ADC (e.g., cell size and permeability). This point is discussed further in [21].
DCE-MRI data was collected using a T1-weighted spoiled gradient echo sequence with TR = 101 ms, TE = 1.9 ms, and a flip angle of 22°. A 200 µL bolus (0.05 mmol kg−1) of Gado-DTPATM (BioPhysics Assay Lab, Worcester, MA, USA) was injected after 25 image volumes were acquired. The relative blood volume fraction, ϕ V ( x ¯ , t ) , was calculated by computing the ratio of the area under the curve for the concentration of the contrast agent time course for each tissue voxel to the area under the arterial input function [34] over the first 60 s [9,21].
Following the first imaging session, a mutual information based rigid registration algorithm was used at the beginning of each subsequent imaging session to register the current animal placement to the first imaging session placement [25]. The resulting spatial offsets and rotations were applied when selecting the field of view on the console to minimize the need for post-acquisition registration.
After all imaging studies were completed for an individual rat, the same registration algorithm was employed (as needed) to maximally align the imaging data across time. The rigid registration transformation was estimated using imregtform in MATLAB R2020a using the multimodal configuration which employs a mutual information-based cost function. Registration performance was visually assessed by overlaying the registered and target (i.e., initial image) over the central eight slices. If the initial automatic registration performance was inadequate, a manual transformation was applied as an initial transformation for imregtform. Tumor regions of interest were segmented using a semi-automated approach by an imaging scientist with over ten years of experience in segmenting contrast-enhancing lesions in murine models of glioma. The segmented tumor consisted of the enhancing lesion on a post-contrast T1-weighted MRI (from the DCE-MRI dataset). First, a region of interest is manually drawn around the contrast-enhancing lesion. Second, a k-means clustering in MATLAB R2020a is used to identify voxels that are enhancing and non-enhancing within the region of interest. Third, imfill in MATLAB R2020a is used to fill in holes within the regions identified as enhancing tissue. Finally, the k-means segmented tumor is visually inspected before proceeding to modeling. The robustness of this approach was evaluated in an in silico study where noise was added from a normal distribution (equivalent to an SNR of 20) to each animal’s day 0 post-contrast T1-weighted MRI to generate 100 unique imaging volumes which were then segmented using the semi-automated approach just described. We then calculated the variability in volume estimates and the degree of spatial overlap using the standard error and Dice correlation coefficient, respectively. We observed that this semi-automated approach is robust to the noise level observed in the image used for segmentation resulting in a standard error of less than 0.26 mm3 and Dice correlation coefficients greater than 0.91 for all animals. Complete results are reported in the Supplemental Material and Table S1. The first panel in Figure 1 provides a schematic of our image processing and modeling approach.

2.3. Numerical and Computational Methods

2.3.1. Finite Difference Approximation

While complete numerical details are described elsewhere [30], here we present the salient details. A finite difference approximation implemented in MATLAB R2020a was used to determine the spatial-temporal evolution of ϕ T ( x ¯ , t ) and ϕ V ( x ¯ , t ) using a fully explicit in time differentiation (time step = 0.01 days) and three dimensions in space (Δx = 250 µm, Δy = 250 µm, Δz = 1000 µm) central difference spatial differentiation. No flux boundary conditions were applied at the skull boundary for ϕ T ( x ¯ , t ) and ϕ V ( x ¯ , t ) at the skull boundary. The boundary condition for u 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).

2.3.2. Parameter Calibration (Scenario 1) and Tumor Response Forecasting (Scenario 2)

As outlined in Section 2.1, we developed a family of 18 models consisting of nine approaches accounting for the effects of RT (three models of RT response (RTM1-RTM3) and three approaches (C1C3) to spatially vary the effect of RT), and two approaches to calibrate kp,T. The model parameters listed in Table 1 were calibrated (second panel in Figure 1) for all 18 models and each animal using a hybrid simulated-annealing Levenberg-Marquardt algorithm [21]. All parameters were considered to be global (uniform throughout the domain) with the exception of kp,T, which was calibrated either as a global variable or varied spatially throughout the domain. When kp,T was defined as a field, parameter values were only calibrated within a subset of points within the tumor and then interpolated elsewhere to reduce the number of individual parameters needed to be calibrated. More specifically, for a given 3 × 3 region of voxels 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. This calibration approach also regularizes (or smooths) the parameter field spatially. (Bounds used for model calibration are reported in Table S2). We performed two calibration scenarios to determine how well the model describes the data and how well it forecasts future response. For scenario 1, we calibrated each model to all imaging days for an individual animal to evaluate how well each model describes the data. Once we have calibrated the model parameters, we then determined the parameters 95% confidence intervals (third panel in Figure 1) using nlparci in MATLAB. We then sampled the parameter confidence intervals to generate 100 sets of model parameters that were then used in 100 additional finite difference forward evaluations to estimate tumor growth and response. For each of the 100 tumor growth calculations, we assessed the error as described below to generate confidence intervals in our model output. For scenario 2, we seek to determine the predictive strength of these models. To do so, we first calibrated each model to the first half of the available data for each animal. We then (similarly to scenario 1) sampled the confidence intervals of the calibrated parameters to generate a set of 100 model parameters that were then used in 100 additional finite difference forecast of tumor growth and treatment response that was directly compared to the last half of the available data for each animal. For each of these 100 forecasts, we also determined the error as described below to generate confidence intervals in our predictions.

2.3.3. Model Selection and Ensemble Average

The Akaike Information Criterion (AIC) [35] was used to select the model that optimally balanced model complexity and agreement with the data. The AIC is defined as:
A I C = 2 k + n ln ( R S S n ) + 2 k ( k + 1 n k 1 )
where k is equal to the number of calibrated parameters for each model, n is the number of data points used to calibrate the model, and RSS is the residual sum squares between the model and measured tumor and blood volume fractions calculated over the total data used for calibration. For each animal in the first calibration scenario (i.e., the scenario where all models are calibrated to the entire time-course of imaging data obtained for each animal), we calculated the AIC for each model. We then report results for the model with the lowest AIC as well as the ensemble average model. The AIC for each model was used to calculate the ensemble average weights defined as:
w i = exp ( δ i 2 ) j = 1 18 exp ( δ j 2 )
where wi is the weight for the i-th model, δ i is equal to A I C i A I C min , and AICmin is the minimum observed AIC. The ensemble averaged tumor volume fraction was calculated using:
ϕ T , e n s = j = 1 100 i = 1 18 w i ϕ T , i , j
where ϕ T , e n s is the ensemble average, and ϕ T , i , j is tumor volume fraction for the i-th model and the j-th set of parameters sampled from the parameter confidence intervals (described in Section 2.3.2). The ensemble averaged blood volume fraction was calculated in a similar way as Equation (16). For the second calibration scenario (i.e., the prediction scenario), we report results using the model with the lowest AIC (calculated for each animal during the second calibration scenario) as well as the ensemble average of the model forecasts.

2.4. Error Quantification

The model calculated tumor and blood volume fractions were compared directly to the measured values (see the fourth panel of Figure 1). At the global level, we calculated the percent error in predicted tumor volume and the Dice coefficient (describing the degree of overlap of the measured and model calculated volumes). At the local level we calculated the Pearson correlation coefficient (PCC) and the concordance correlation coefficient (CCC), which characterize the degree of correlation and agreement, respectively, between the model and the measurement at each voxel location. For each model calculated tumor and blood volume fractions derived from each set of model parameters (sampled from the parameter confidence interval), we assessed the global and local level errors at each time point. Each of the above-mentioned error metrics were then averaged for each set of model parameters across time. That is, for each model and animal we have 100 estimates of PCC, CCC, Dice, and percent error in tumor volume. The results for each animal are reported as box and whisker plots.

3. Results

3.1. Model Calibration Scenario

Figure 2 and Figure 3 report the global and local-level error analysis, respectively, for the first calibration scenario. The top panel in Figure 2 shows the model weights for each individual animal as well as the average weight across animals. In general, models with a voxel-specific kp,T (models 1–9) were weighted more heavily (contributing 80.7% of the ensemble weight) compared to the global kp,T. Additionally, the models that incorporated a death rate (rather than instantaneous death) were weighted the most (models 7–9) and contributed 56.9% of the ensemble weight.
At the global level (bottom two panels of Figure 2), the selected model (orange box plots) had median percent errors in tumor volume predictions less than 8.7%, while the ensemble averaged model had percent errors in tumor volume predictions less than 18.1%. A high level of spatial overlap was observed across all animals resulting in Dice values greater than 0.68. Six out of seven animals had statistically significant differences (p < 0.05) between the selected model and the ensemble average for both the percent error in tumor volume and Dice values.
Figure 3 reports the local-level error analysis for both tumor (left column) and blood volume estimates (right column). At the local level for both the selected and ensemble average model (Figure 3), we observed a strong degree of correlation (PCCs > 0.74) for all animals between the observed and model estimated values of tumor volume fraction. Similarly, a high degree of agreement (CCCs > 0.76) was observed for animals 2–4, while the remaining animals, comparatively, had a lower degree of agreement (CCCs > 0.53). All seven animals had statistically significant differences (p < 0.05) between the selected and ensemble average model.
For blood volume estimates, a high degree of correlation (PCCs > 0.70) was observed for the selected model for all animals, while lower correlation was observed for the ensemble average for animals 5 and 6 (PCCs < 0.56). CCCs greater than 0.50 were observed for the selected model estimate of blood volume fraction for all but animal 5. The ensemble average generally resulted in lower agreement (CCCs > 0.14) compared to the selected model. Statistically significant differences (p < 0.05) were observed between the selected and ensemble model for six out of seven animals. Calibrated model parameters for each animal are reported in Table 3 for the selected model.

3.2. Model Prediction Scenario

Figure 4 and Figure 5 report the global and local error analysis, respectively, for the prediction scenario for animals 1, 3, 5, 6, and 7. The top panel in Figure 4 shows the model weights for each individual animal as well as the average weight across animals. Similar to the first calibration scenario, models with a voxel-specific kp,T (models 1–9) were weighted more heavily (contributing 77.5% of the ensemble weight) compared to the global kp,T. Likewise, models that incorporated a death rate (rather than instantaneous death) were weighted the most (models 7–9) and contributed to 54.9% of the ensemble weight.
In general, a low error (less than 16.2%) was observed in tumor volume predictions across all animals for the selected model. The ensemble average model resulted in significantly higher error (p < 0.05) for all animals. Additionally, the Dice correlation coefficient was greater than 0.60 for the selected model. The ensemble average model resulted in significantly lower (p < 0.05) Dice correlation coefficients for four of the animals.
Predictions of the tumor volume fraction at the local level (left column in Figure 5) resulted in strong correlation for animals 1 and 3 (PCCs > 0.80) and lower correlation for the remaining animals (PCCs > 0.57) for the selected model. The ensemble average model resulted in PCCs greater than 0.69 for all animals. Reduced agreement compared to the calibration scenario was observed resulting in CCCs greater than 0.58 for animals 1, 3, and 7 for both the selected and ensemble average model. The ensemble average model resulted in lower agreement compared to the selected model for three out of five animals.
For blood volume fraction predictions at the local level (right column in Figure 5), the selected and ensemble average models both resulted in strong correlation (PCCs > 0.65). Lower agreement was observed across the cohort resulting in a median CCC greater than 0.49 for the selected model and 0.36 for the ensemble average model.
Figure 6 and Figure 7 display the results for the tumor growth and response predictions. Figure 6 and Figure 7 display the central slice predictions of tumor and blood volume fraction, respectively, for all five animals. Figure 6 shows the results for animals 1 and 3, while Figure 7 shows the results for animals 5 through 7. Sagittal views of Figure 6 and Figure 7 are shown in Figures S1 and S2. Figure 8 details the tumor volume predictions overtime for both the selected model and ensemble average for each animal presented in Figure 6 and Figure 7. Time-resolved errors in tumor volume and the Dice correlation coefficient are shown in Figures S3 and S4. For animal 1, a high level of voxel-wise agreement was observed for tumor volume predictions (PCC = 0.86 and CCC = 0.80) compared to the blood volume predictions (PCC = 0.70 and CCC = 0.42) for the selected model. Notably, both models were able to predict the area of low cell density in the first four prediction time points. For animal 3, the model tends to underestimate the tumor area in the central slice; however, a high level of agreement is observed at the voxel level (PCC > 0.84 and CCC > 0.72) for both the selected and ensemble model. For animal 5, we observed a median error less than 8.7% for both models in tumor volume predictions, and strong (PCC = 0.93 and CCC = 0.74) voxel level agreement for the ensemble average model. Animal 6, had predictions over 11 days which resulted in a median error of 16.0% in tumor volume predictions for the model with the lowest AIC, while the ensemble average model had greater than 100% error in tumor predictions. For animal 7, we predicted tumor growth from days 10 to 17, which resulted in a median error less than 21.5% for both models in tumor volume predictions. Additionally, both models tended to overestimate blood volume predictions at day 14 and 17 resulting in CCCs less than 0.26.

4. Discussion

We have developed and applied an experimental-computational framework to predict the response of murine tumors to fractionated radiation on an individual subject basis. More specifically, we parameterized 18 biologically-based models of tumor and vasculature response to fractionated RT via quantitative MRI data obtained in seven animals, and then used the calibrated models to predict the spatio-temporal development of key tumor characteristics. Non-invasive imaging data from DW- and DCE-MRI enabled estimates of tumor and blood volume fractions that served as ground truth for model calibration, selection, and evaluation of prediction accuracy. We evaluated two calibration scenarios to assess how well these models: (1) describe the entire tumor growth time course, as well as how well they (2) predict the remaining imaging visits when half of the data is used for calibration. For the first scenario, when calibrated to all available imaging time points, the model with the lowest AIC resulted in less than 8.7% error in tumor volume estimates across all animals. While the ensemble average model resulted in less than 18.1% error in tumor volume estimates across all animals for the first scenario. Notably, the ensemble average model resulted in a decrease in tumor volume error for three out of seven animals, and an increase in CCCs for all seven animals. For the second scenario, when we calibrate over the first half of the time points and predict the remaining time points, a median error of less than 16.2% in tumor volume estimates were observed across all animals for the model with lowest AIC. Unlike the first calibration scenario, the ensemble average prediction resulted in increased tumor volume estimates compared to the model predictions with the selected model. However, increased correlation (PCCs) and agreement (CCCs) were observed for three out of the five and two out of the five animals, respectively, for the ensemble average prediction. We note, that we did not observe significantly different response between treatment groups. Notably animals 1, 6, 7 were all imaged at least 5 days post RT despite receiving different treatment doses. We hypothesize the varied response could be due to variations in tumor growth properties or sensitivity to RT itself. For example, animal 1 presents with the lowest growth rate (Table 3) of the animals, had the longest overall survival, and received only 2 Gy/day. This result may provide further motivation for personalizing RT regimens to adapt to individual tumor properties. This preliminary study indicates a promising approach for personalizing mathematical models of response to fractionated RT.
Our approach extends our previously developed image-driven modeling framework applied in the presence [21,25] and absence of radiotherapy [9] in the C6 glioma line. In the present effort, we refined our experimental approach in three main ways to improve upon previous studies. First, we applied our experimental-computational framework to the U-87 cell line. Second, animals received smaller radiation fractions (2 or 4 Gy) instead of single large dose (20 or 40 Gy) to more closely mimic how RT is delivered clinically. While the delivery of RT is still substantially different (whole brain versus highly conformal RT beams), this approach facilitated the study of temporal response to RT that is not possible through single fraction of RT. Third, we constructed for each animal an ensemble average model weighted using the AIC. In general, the ensemble average model did not outperform the selected model, but by weighting model outcomes we have a forecast that considers all of the possible tumor response patterns. In addition, the model weights themselves provide some insight into the subject-specific growth and response characteristics. We investigated ensemble averages (and individual model simulations) as they could deliver a powerful tool for clinicians providing a forecast of the average, best, and worst response scenarios for an individual subject at the beginning of therapy. Similar to weather forecasting [36], as additional data is observed model weights could be adjusted based on forecast performance to provide an updated ensemble forecast.
Accurate characterizations and forecasts of the temporal response of tumors to radiation therapy are critical to the development of patient optimized treatment regimens. Optimized radiotherapy regimens may be able to address variations in tumor properties (e.g., hypoxia) that alter radiosensitivity and influence response to radiotherapy [37]. Several promising computational oncology studies have investigated applying mathematical modeling to systematically evaluate alternative regimens [17,38,39]. These approaches focus primarily on adapting or optimizing regimens based on changes in tumor geometry. However, we posit that these approaches could be refined further to adapt or optimize RT regimens based on both tumor geometry and intratumor radiosensitivity. By reducing tumor biology to simply the shape and location of the tumor, we are blind to spatial variations in tumor response that might ultimately lead to disease progression. Our coupled model of tumor growth and angiogenesis forecasts a spatial map of response that could be targeted by intensity modulated RT. The results in this preliminary cohort indicate strong predictive accuracy in tumor geometry (low error in tumor volume predictions) while additional development is needed to improve local level predictions. Our current model of response to RT is relatively simple compared to the complex biophysical mechanisms of response to radiotherapy, which are undoubtedly varying spatially and temporal during fractionated RT. However, additional experimental data is required to properly initialize and constrain a more complete description of tumor radiobiology. In the clinical setting, one potential application for this modeling framework is to predict long term response (or time to progression) to guide alternative fractionation schemes with the end goal of improving patient outcomes [17,28]. While short-term predictions performed well, further development may be needed to improve long term predictions. Once a model is established that can accurately reproduce the spatiotemporal development and response of the tumor it can then be used to identify alterative treatment regimens or fractionation schemes that may outperform standard methods.
With regards to calibrated parameter values, we observed that the calibrated kp,T,0 and DT,0 are within ranges reported in literature for image-based estimates of proliferation and diffusion [25,40] except for animal 6. For animal 6, kp,T,0 and DT,0 were at least 2.3 and 1.5 times higher, respectively, than the other animals. The larger growth rates for animal 6 might contribute to the higher errors observed in both the calibration and prediction scenarios. The poorer prediction may be due to the rapid or aggressive tumor early on that results in higher proliferation and diffusion coefficients. Individual tumor forecasts could be compared to a population averaged forecast to identify subjects whose tumor growth or response forecasts deviate significantly from the population. Additionally, we observed a dose-dependency on SF where animals that received 2 Gy was higher (SF > 0.95) than those animals that received 4 Gy (SF < 0.91).
There are several opportunities for further investigation and development of our experimental-computational approach. First, we assume that the measured ADC can provide reasonable estimates of tumor volume fraction. As we have previously discussed [21], the ADC is influenced by a combination of both cellular and tissue properties including cell density, cell size, and cell permeability. However, in our formulation, we assume the changes observed in the ADC over time are influenced predominately by cell density. Advanced diffusion techniques [41] may be able to probe these properties further and provide a more complete description of tumor tissue and cell properties. Second, while the present data set does provide evidence that personalizing mathematical models of response to fractionated RT is a promising avenue to investigate, further experimental studies are needed to increase the cohort size and provide a more complete characterization of the predictive accuracy of these models. Third, future studies should consider including additional tumor cell lines that are more infiltrative and better recapitulate human glioblastoma characteristics, thereby testing the generalizability of the experimental-computational approach. While the U-87 line, is a human-derived glioblastoma line, it may fail to capture some of the key characteristics of the human disease [42]. Notably, U-87 tumors tend not to be diffuse or infiltrative and remain well circumscribed and vascularized. As such, the U-87 line enables reliable estimation of tumor burden (and assignment of ground truth) making it suitable for model development and refinement. However, alternative cell lines that better recapitulate human glioblastoma should be considered for further development of this technique. Other human derived cell lines which have a more invasive pathology (e.g., the U-251) should be considered to evaluate the generalizability of this approach to different physio-pathological conditions. One challenge with any of these human derived cells is the use of immunocompromised animals which may not accurately recapitulate the interactions between the tumor and host, or the immune response to RT [43]. Fourth, genetic and phenotypic diversity [44] is a significant factor in tumor growth and response of tumors to systemic therapy and radiotherapy. Although this diversity is not explicitly described in our mathematical model, we hypothesize that these variations may be implicitly captured through the individualized model parameterization. To some extent, phenotypic diversity may be accounted for through calibrating a proliferation field rather than assuming homogenous tumor growth properties. To test this hypothesis, tumor growth predictions should be compared to other approaches that explicitly integrate genetic and phenotypic diversity with mechanism-based models [45]. Fifth, future studies, in larger cohorts, are needed to investigate whether it is important to perform subject-specific model selection or use an ensemble average model whose weights are determined from a training set. Sixth, future efforts should consider more complete descriptions of tumor-induced angiogenesis and regression. Our model of tumor-induced angiogenesis and regression is an over-simplification of the complex mechanisms of vasculature creation and destruction observed in vivo [46]. As such, it may fail to accurately capture angiogenesis in all tumors. Finally, in this modeling approach we assumed the delivered dose to be uniform throughout the tumor, which may not be the case for a single angle delivery of RT. Variability in the delivered dose should be considered to more accurately capture response dynamics spatially. In the clinical setting, RT treatment plans could be used to provide a spatial map of the planned RT throughout the brain.

5. Conclusions

We have developed and applied a novel image-driven and biologically-based modeling framework for characterizing and forecasting response of both tumor and vasculature tissue to fractionated radiation therapy. Preliminary results indicate low error in characterizing and predicting tumor volume and local cell number. Thus, further investigation is warranted, and future efforts should apply this approach to a larger cohort and a broader range of fractionation schemes.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/cancers13081765/s1, Figure S1: Model predictions (scenario 2) for animal 1 and 3 (Sagittal plane), Figure S2: Model predictions (scenario 2) for animal 5, 6, and 7 (Sagittal plane), Figure S3: Error for tumor growth predictions for animals 1 and 3. Figure S4: Error for tumor growth predictions for animals 5, 6, and 7. Table S1: Standard error in segmented volume, Table S2: Parameter ranges used for model calibration.

Author Contributions

Conceptualization, D.A.H.II and T.E.Y.; methodology, D.A.H.II, A.M.J., T.D., and T.E.Y.; software, D.A.H.II; validation, D.A.H.II, A.M.J.; formal analysis, D.A.H.II, A.M.J., and T.E.Y.; investigation, D.A.H.II; resources, D.A.H.II; data curation, D.A.H.II; writing—original draft preparation, D.A.H.II; writing—review and editing, A.M.J., T.D., and T.E.Y.; visualization, D.A.H.II; supervision, T.E.Y.; project administration, D.A.H.II, T.E.Y.; funding acquisition, D.A.H.II, T.E.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Cancer Institute, grant numbers R01CA235800, U24CA226110, U01CA174706, and U01CA253540, the Cancer Prevention and Research Institute of Texas grant number RR160005, and the AAPM via the 2018-19 Research Seed Funding.

Institutional Review Board Statement

All experimental procedures were approved by The University of Texas at Austin’s Institutional Animal Care and use Committee (protocol No. AUP-2018-00338 on 11 February 2019).

Informed Consent Statement

Not applicable.

Data Availability Statement

The datasets generated or analyzed, and the code used during the current study are available from the corresponding author on reasonable request.

Acknowledgments

The authors acknowledge the Texas Advanced Computing Center for providing high-performance computing resources. T.E.Y. is a CPRIT Scholar of Cancer Research.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Omuro, A.; DeAngelis, L.M. Glioblastoma and other malignant gliomas: A clinical review. JAMA 2013, 310, 1842–1850. [Google Scholar] [CrossRef]
  2. Lagendijk, J.J.W.; Raaymakers, B.W.; Van Vulpen, M. The Magnetic Resonance Imaging–Linac System. Semin. Radiat. Oncol. 2014, 24, 207–209. [Google Scholar] [CrossRef] [PubMed]
  3. Jaffray, D.A. Image-guided radiotherapy: From current concept to future perspectives. Nat. Rev. Clin. Oncol. 2012, 9, 688–699. [Google Scholar] [CrossRef]
  4. Jones, K.M.; Michel, K.A.; Bankson, J.A.; Fuller, C.D.; Klopp, A.H.; Venkatesan, A.M. Emerging Magnetic Resonance Imaging Technologies for Radiation Therapy Planning and Response Assessment. Int. J. Radiat. Oncol. Biol. Phys. 2018, 101, 1046–1056. [Google Scholar] [CrossRef] [Green Version]
  5. Troost, E.G.C.; Thorwarth, D.; Oyen, W.J.G. Imaging-Based Treatment Adaptation in Radiation Oncology. J. Nucl. Med. 2015, 56, 1922–1929. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Yankeelov, T.E.; Atuegwu, N.; Hormuth, D.A.; Weis, J.A.; Barnes, S.L.; Miga, M.I.; Rericha, E.C.; Quaranta, V. Clinically Relevant Modeling of Tumor Growth and Treatment Response. Sci. Transl. Med. 2013, 5, 187ps9. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. McMahon, S.J. The linear quadratic model: Usage, interpretation and challenges. Phys. Med. Biol. 2019, 64, 01TR01. [Google Scholar] [CrossRef]
  8. Wouters, B.G. Cell death after irradiation: How, when and why cells die. In Basic Clinical Radiobiology; Joiner, M.C., Van der Kogel, A.J., Eds.; CRC Press: Boca Raton, FL, USA, 2009; pp. 27–40. [Google Scholar]
  9. Hormuth, D.A.; Jarrett, A.M.; Yankeelov, T.E. Forecasting tumor and vasculature response dynamics to radiation therapy via image based mathematical modeling. Radiat. Oncol. 2020, 15, 4–14. [Google Scholar] [CrossRef] [PubMed]
  10. Enderling, H.; Chaplain, M.A.J.; Hahnfeldt, P. Quantitative modeling of tumor dynamics and radiotherapy. Acta Biotheor. 2010, 58, 341–353. [Google Scholar] [CrossRef] [PubMed]
  11. Brüningk, S.; Powathil, G.; Ziegenhein, P.; Ijaz, J.; Rivens, I.; Nill, S.; Chaplain, M.; Oelfke, U.; Ter Haar, G. Combining radiation with hyperthermia: A multiscale model informed by in vitro experiments. J. R. Soc. Interface 2018, 15, 20170681. [Google Scholar] [CrossRef] [Green Version]
  12. Prokopiou, S.; Moros, E.G.; Poleszczuk, J.; Caudell, J.; Torres-Roca, J.F.; Latifi, K.; Lee, J.K.; Myerson, R.; Harrison, L.B.; Enderling, H. A proliferation saturation index to predict radiation response and personalize radiotherapy fractionation. Radiat. Oncol. 2015, 10, 159. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Rockne, R.; Rockhill, J.K.; Mrugala, M.; Spence, A.M.; Kalet, I.; Hendrickson, K.; Lai, A.; Cloughesy, T.; Alvord, E.C.; Swanson, K.R. Predicting the efficacy of radiotherapy in individual glioblastoma patients in vivo: A mathematical modeling approach. Phys. Med. Biol. 2010, 55, 3271–3285. [Google Scholar] [CrossRef] [Green Version]
  14. Elazab, A.; Anter, A.M.; Bai, H.; Hu, Q.; Hussain, Z.; Ni, D.; Wang, T.; Lei, B. An optimized generic cerebral tumor growth modeling framework by coupling biomechanical and diffusive models with treatment effects. Appl. Soft Comput. J. 2019, 80, 617–627. [Google Scholar] [CrossRef]
  15. Mizutani, T.; Magome, T.; Igaki, H.; Haga, A.; Nawa, K.; Sekiya, N.; Nakagawa, K. Optimization of treatment strategy by using a machine learning model to predict survival time of patients with malignant glioma after radiotherapy. J. Radiat. Res. 2019, 60, 818–824. [Google Scholar] [CrossRef] [PubMed]
  16. Holdsworth, C.; Corwin, D.; Stewart, R.; Rockne, R.; Trister, A.; Swanson, K.; Philips, M. Adaptive IMRT using a multiobjective evolutionary algorithm integrated with a diffusion–invasion model of glioblastoma. Phys. Med. Biol. 2012, 57, 8271–8283. [Google Scholar] [CrossRef]
  17. Corwin, D.; Holdsworth, C.; Rockne, R.C.; Trister, A.D.; Mrugala, M.M.; Rockhill, J.K.; Stewart, R.D.; Phillips, M.; Swanson, K.R. Toward Patient-Specific, Biologically Optimized Radiation Therapy Plans for the Treatment of Glioblastoma. PLoS ONE 2013, 8, e079115. [Google Scholar] [CrossRef]
  18. Kim, M.; Kotas, J.; Rockhill, J.; Phillips, M. A Feasibility Study of Personalized Prescription Schemes for Glioblastoma Patients Using a Proliferation and Invasion Glioma Model. Cancers 2017, 9, 51. [Google Scholar] [CrossRef] [Green Version]
  19. Budia, I.; Alvarez-Arenas, A.; Woolley, T.E.; Calvo, G.F.; Belmonte-Beitia, J. Radiation protraction schedules for low-grade gliomas: A comparison between different mathematical models. J. R. Soc. Interface 2019, 16. [Google Scholar] [CrossRef] [Green Version]
  20. Sunassee, E.D.; Tan, D.; Ji, N.; Brady, R.; Moros, E.G.; Caudell, J.J.; Yartsev, S.; Enderling, H. Proliferation saturation index in an adaptive Bayesian approach to predict patient-specific radiotherapy responses. Int. J. Radiat. Biol. 2019, 95, 1421–1426. [Google Scholar] [CrossRef] [Green Version]
  21. Hormuth, D.A.; Jarrett, A.M.; Feng, X.; Yankeelov, T.E. Calibrating a Predictive Model of Tumor Growth and Angiogenesis with Quantitative MRI. Ann. Biomed. Eng. 2019, 47, 1539–1551. [Google Scholar] [CrossRef]
  22. Stupp, R.; Mason, W.P.; Van den Bent, M.J.; Weller, M.; Fisher, B.; Taphoorn, M.J.B.; Belanger, K.; Brandes, A.A.; Marosi, C.; Bogdahn, U.; et al. Radiotherapy plus Concomitant and Adjuvant Temozolomide for Glioblastoma. N. Engl. J. Med. 2005, 352, 987–996. [Google Scholar] [CrossRef] [PubMed]
  23. Padhani, A.R.; Liu, G.; Mu-Koh, D.; Chenevert, T.L.; Thoeny, H.C.; Takahara, T.; Dzik-Jurasz, A.; Ross, B.D.; Van Cauteren, M.; Collins, D.; et al. Diffusion-Weighted Magnetic Resonance Imaging as a Cancer Biomarker: Consensus and Recommendations. Neoplasia 2009, 11, 102–125. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. O’Connor, J.P.B.; Tofts, P.S.; Miles, K.A.; Parkes, L.M.; Thompson, G.; Jackson, A. Dynamic contrast-enhanced imaging techniques: CT and MRI. Br. J. Radiol. 2011, 84, S112–S120. [Google Scholar] [CrossRef] [Green Version]
  25. Hormuth, D.A., II; Weis, J.A.; Barnes, S.L.; Miga, M.I.; Rericha, E.C.; Quaranta, V.; Yankeelov, T.E. Predicting in vivo glioma growth with the reaction diffusion equation constrained by quantitative magnetic resonance imaging data. Phys. Biol. 2015, 12, 46006. [Google Scholar] [CrossRef] [Green Version]
  26. Alfonso, J.C.L.; Talkenberger, K.; Seifert, M.; Klink, B.; Hawkins-Daarud, A.; Swanson, K.R.; Hatzikirou, H.; Deutsch, A. The biology and mathematical modelling of glioma invasion: A review. J. R. Soc. Interface 2017, 14, 20170490. [Google Scholar] [CrossRef] [PubMed]
  27. Rockne, R.; Alvord, E.C.; Rockhill, J.K.; Swanson, K.R. A mathematical model for brain tumor response to radiation therapy. J. Math. Biol. 2009, 58, 561–578. [Google Scholar] [CrossRef] [Green Version]
  28. Jarrett, A.M.; Hormuth, D.A.; Wu, C.; Syed, A.K.; Virostko, J.; Sorace, A.; DiCarlo, J.C.; Patt, D.; Goodgame, B.; Avery, S.; et al. Evaluating patient-specific neoadjuvant regimens for breast cancer via a mathematical model constrained by quantitative magnetic resonance imaging data. Neoplasia 2020, 22, 820–830. [Google Scholar] [CrossRef]
  29. Chen, X.; Summers, R.; Jianhua Yao, J. FEM-Based 3-D Tumor Growth Prediction for Kidney Tumor. IEEE Trans. Biomed. Eng. 2011, 58, 463–467. [Google Scholar] [CrossRef]
  30. Hormuth, D., II; Eldridge, S.B.; Weis, J.; Miga, M.I.; Yankeelov, T.E. Mechanically Coupled Reaction-Diffusion Model to Predict Glioma Growth: Methodological Details. In Springer Methods and Protocols: Cancer Systems Biology; Von Stechow, L., Ed.; Springer: New York, NY, USA, 2018; pp. 225–241. ISBN 978-1-4939-7493-1. [Google Scholar]
  31. Jove Rodent Stereotaxic Surgery. Available online: https://www.jove.com/t/5205/rodent-stereotaxic-surgery (accessed on 23 March 2021).
  32. Steinhoff, S.; Zaitsev, M.; Zilles, K.; Shah, N.J. Fast T(1) mapping with volume coverage. Magn. Reson. Med. 2001, 46, 131–140. [Google Scholar] [CrossRef]
  33. Whisenant, J.G.; Ayers, G.D.; Loveless, M.E.; Barnes, S.L.; Colvin, D.C.; Yankeelov, T.E. Assessing reproducibility of diffusion-weighted magnetic resonance imaging studies in a murine model of HER2+ breast cancer. Magn. Reson. Imaging 2014, 32, 245–249. [Google Scholar] [CrossRef] [Green Version]
  34. Hormuth, D.A., II; Skinner, J.T.; Does, M.D.; Yankeelov, T.E. A comparison of individual and population-derived vascular input functions for quantitative DCE-MRI in rats. Magn. Reson. Imaging 2014, 32, 397–401. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Akaike, H. A new look at the statistical model identification. IEEE Trans. Autom. Control 1974, 19, 716–723. [Google Scholar] [CrossRef]
  36. Yankeelov, T.E.; Quaranta, V.; Evans, K.J.; Rericha, E.C. Toward a Science of Tumor Forecasting for Clinical Oncology. Cancer Res. 2015, 75, 918–923. [Google Scholar] [CrossRef] [Green Version]
  37. Kallman, R.F.; Dorie, M.J. Tumor oxygenation and reoxygenation during radiation theraphy: Their importance in predicting tumor response. Int. J. Radiat. Oncol. 1986, 12, 681–685. [Google Scholar] [CrossRef]
  38. Le, M.; Delingette, H.; Kalpathy-Cramer, J.; Gerstner, E.R.; Batchelor, T.; Unkelbach, J.; Ayache, N. Personalized Radiotherapy Planning Based on a Computational Tumor Growth Model. IEEE Trans. Med. Imaging 2016, 36, 815–825. [Google Scholar] [CrossRef]
  39. López Alfonso, J.C.; Parsai, S.; Joshi, N.; Godley, A.; Shah, C.; Koyfman, S.A.; Caudell, J.J.; Fuller, C.D.; Enderling, H.; Scott, J.G. Temporally feathered intensity-modulated radiation therapy: A planning technique to reduce normal tissue toxicity. Med. Phys. 2018, 45, 3466–3474. [Google Scholar] [CrossRef]
  40. Rutter, E.M.; Stepien, T.L.; Anderies, B.J.; Plasencia, J.D.; Woolf, E.C.; Scheck, A.C.; Turner, G.H.; Liu, Q.; Frakes, D.; Kodibagkar, V.; et al. Mathematical Analysis of Glioma Growth in a Murine Model. Sci. Rep. 2017, 7, 2508. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  41. Jiang, X.; Li, H.; Xie, J.; McKinley, E.T.; Zhao, P.; Gore, J.C.; Xu, J. In vivo imaging of cancer cell size and cellularity using temporal diffusion spectroscopy. Magn. Reson. Med. 2017, 78, 156–164. [Google Scholar] [CrossRef]
  42. Jacobs, V.L.; Valdes, P.A.; Hickey, W.F.; De Leo, J.A. Current review of in vivo GBM rodent models: Emphasis on the CNS-1 tumour model. ASN Neuro 2011, 3, e00063. [Google Scholar] [CrossRef] [Green Version]
  43. Willey, C.D.; Gilbert, A.N.; Anderson, J.C.; Gillespie, G.Y. Patient-Derived Xenografts as a Model System for Radiation Research. Semin. Radiat. Oncol. 2015, 25, 273–280. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Schulz, A.; Meyer, F.; Dubrovska, A.; Borgmann, K. Cancer stem cells and radioresistance: DNA repair and beyond. Cancers 2019, 11, 862. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Johnson, K.; Howard, G.R.; Morgan, D.; Brenner, E.A.; Gardner, A.L.; Durrett, R.E.; Mo, W.; Al’Khafaji, A.; Sontag, E.D.; Jarrett, A.M.; et al. Integrating multimodal data sets into a mathematical framework to describe and predict therapeutic resistance in cancer. Phys. Biol. 2021, 18, 016001. [Google Scholar] [CrossRef] [PubMed]
  46. Hardee, M.E.; Zagzag, D. Mechanisms of glioma-associated neovascularization. Am. J. Pathol. 2012, 181, 1126–1141. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Figure 1. Overview of the image processing and computational approach. Pre-processing: Images are first aligned to the first imaging session using a rigid registration. We then analyze DW-MRI and DCE-MRI to estimate the ADC and blood volume, respectively. Contrast-enhanced images are then used for the segmentation of tumor burden. Calibrate model family for each animal: Each model within the model family is calibrated for each individual animal. Tumor and blood volume fractions at baseline (t0) and an initial guess of model parameters, (a), are used to initialize our 3D model of tumor growth and response. The mathematical model is evaluated at subsequent imaging days t1 to tn, (b), and the error between the model and the measured tumor and blood volume fractions are assessed and are used to refine the estimates of model parameters (c). Post-calibration processing: Confidence intervals for each set of model parameters are determined for each model and animal. We then sample the parameter confidence intervals and evaluate the model with each set of sampled parameters to assess the variations in model outcomes. The Akaike Information Criterion is then used to select the most parsimonious model and determine model weights. The ensemble averaged model is then constructed by weighting each model output. Error analysis: Error between the measured and model estimated tumor and blood volume fraction is then analyzed at the global and local levels.
Figure 1. Overview of the image processing and computational approach. Pre-processing: Images are first aligned to the first imaging session using a rigid registration. We then analyze DW-MRI and DCE-MRI to estimate the ADC and blood volume, respectively. Contrast-enhanced images are then used for the segmentation of tumor burden. Calibrate model family for each animal: Each model within the model family is calibrated for each individual animal. Tumor and blood volume fractions at baseline (t0) and an initial guess of model parameters, (a), are used to initialize our 3D model of tumor growth and response. The mathematical model is evaluated at subsequent imaging days t1 to tn, (b), and the error between the model and the measured tumor and blood volume fractions are assessed and are used to refine the estimates of model parameters (c). Post-calibration processing: Confidence intervals for each set of model parameters are determined for each model and animal. We then sample the parameter confidence intervals and evaluate the model with each set of sampled parameters to assess the variations in model outcomes. The Akaike Information Criterion is then used to select the most parsimonious model and determine model weights. The ensemble averaged model is then constructed by weighting each model output. Error analysis: Error between the measured and model estimated tumor and blood volume fraction is then analyzed at the global and local levels.
Cancers 13 01765 g001
Figure 2. Global error analysis from the calibration scenario. The top panel shows model weights determined from the Akaike Information Criterion for all seven animals. C1 through C3 correspond to the different approaches to introduce spatial variation in tumor response, while RTM1 through RTM3 correspond to the different RT response models. Models with a local kp,T and RTM3 accounted for 56.9% of the ensemble average model. In the bottom panel, the average percent error in tumor volume and the Dice correlation coefficient are reported for the selected model (orange) and the ensemble averaged model (blue) for all seven animals. The median error for the selected model resulted in less than 9% error for all animals. The median error for the ensemble average was also less than 6% error, with the exception of animal 6 which had greater than 19% error. A strong level of spatial overlap (Dice > 0.68) was observed for all animals. Statistically significant differences (p < 0.05) between model estimates are indicated by the ‘+’ symbol. Model 7 was weighted the highest with an average weight of 0.21.
Figure 2. Global error analysis from the calibration scenario. The top panel shows model weights determined from the Akaike Information Criterion for all seven animals. C1 through C3 correspond to the different approaches to introduce spatial variation in tumor response, while RTM1 through RTM3 correspond to the different RT response models. Models with a local kp,T and RTM3 accounted for 56.9% of the ensemble average model. In the bottom panel, the average percent error in tumor volume and the Dice correlation coefficient are reported for the selected model (orange) and the ensemble averaged model (blue) for all seven animals. The median error for the selected model resulted in less than 9% error for all animals. The median error for the ensemble average was also less than 6% error, with the exception of animal 6 which had greater than 19% error. A strong level of spatial overlap (Dice > 0.68) was observed for all animals. Statistically significant differences (p < 0.05) between model estimates are indicated by the ‘+’ symbol. Model 7 was weighted the highest with an average weight of 0.21.
Cancers 13 01765 g002
Figure 3. Local error analysis from the calibration scenario. The PCC (top row) and CCC (bottom row) are reported for the tumor (left column) and blood volume (right column) fractions for the calibration scenario. In general, a high level of agreement and correlation was observed for the tumor volume fraction (PCCs > 0.74 and CCCs > 0.53) for the selected and ensemble average model estimates. Similarly, a high level of correlation was observed for voxel-wise estimates of blood volume fraction (PCCs > 0.70). CCCs were greater than 0.50 for all but animal five. Statistically significant differences (p < 0.05) between model estimates are indicated by the ‘+’ symbol.
Figure 3. Local error analysis from the calibration scenario. The PCC (top row) and CCC (bottom row) are reported for the tumor (left column) and blood volume (right column) fractions for the calibration scenario. In general, a high level of agreement and correlation was observed for the tumor volume fraction (PCCs > 0.74 and CCCs > 0.53) for the selected and ensemble average model estimates. Similarly, a high level of correlation was observed for voxel-wise estimates of blood volume fraction (PCCs > 0.70). CCCs were greater than 0.50 for all but animal five. Statistically significant differences (p < 0.05) between model estimates are indicated by the ‘+’ symbol.
Cancers 13 01765 g003
Figure 4. Global error analysis from the prediction scenario. The top panel shows model weights determined from the Akaike Information Criterion for the five animals from the prediction scenario. C1 through C3 correspond to the different approaches to introduce spatial variation in tumor response, while RTM1 through RTM3 correspond to the different RT response models. Models with a local kp,T and RTM3 accounted for 54.9% of the ensemble average model. In the bottom panel, the average percent error in tumor volume and the Dice correlation coefficient are reported for the selected model (orange) and the ensemble averaged model (blue) for five animals. (Animals 2 and 4 had insufficient imaging visits for the prediction scenario). The median error for the selected model resulted in less than 16.2% error for all animals. For animal 1 and 6, the ensemble average greatly overestimated the tumor volume with a median error of greater than 100%. A high level of spatial overlap (Dice values greater than 0.60) was observed for the selected model. Statistically significant differences (p < 0.05) between model estimates are indicated by the ‘+’ symbol.
Figure 4. Global error analysis from the prediction scenario. The top panel shows model weights determined from the Akaike Information Criterion for the five animals from the prediction scenario. C1 through C3 correspond to the different approaches to introduce spatial variation in tumor response, while RTM1 through RTM3 correspond to the different RT response models. Models with a local kp,T and RTM3 accounted for 54.9% of the ensemble average model. In the bottom panel, the average percent error in tumor volume and the Dice correlation coefficient are reported for the selected model (orange) and the ensemble averaged model (blue) for five animals. (Animals 2 and 4 had insufficient imaging visits for the prediction scenario). The median error for the selected model resulted in less than 16.2% error for all animals. For animal 1 and 6, the ensemble average greatly overestimated the tumor volume with a median error of greater than 100%. A high level of spatial overlap (Dice values greater than 0.60) was observed for the selected model. Statistically significant differences (p < 0.05) between model estimates are indicated by the ‘+’ symbol.
Cancers 13 01765 g004
Figure 5. Local error analysis from the prediction scenario. The Pearson correlation coefficient (top row) and concordance correlation coefficient (bottom row) are reported for the tumor (left column) and blood volume (right column) fractions for the prediction scenario for five animals. (Animals 2 and 4 had insufficient imaging visits for the prediction scenario). Lower PCCs were observed for animal 5 and 6 for the selected model compared to the rest of the cohort. Elevated PCCs were observed for three out of five animals. Similarly, voxel level agreement was greater than 0.58 for animals 1, 3, and 7 for the selected model. Blood volume predictions at the voxel level had a high level of correlation (PCCs greater than 0.65) for both the selected and ensemble averaged model. Lower agreement was observed between the predicted and measured blood volume resulting in CCCs greater than 0.49 for the selected model and 0.36 for the ensemble model.
Figure 5. Local error analysis from the prediction scenario. The Pearson correlation coefficient (top row) and concordance correlation coefficient (bottom row) are reported for the tumor (left column) and blood volume (right column) fractions for the prediction scenario for five animals. (Animals 2 and 4 had insufficient imaging visits for the prediction scenario). Lower PCCs were observed for animal 5 and 6 for the selected model compared to the rest of the cohort. Elevated PCCs were observed for three out of five animals. Similarly, voxel level agreement was greater than 0.58 for animals 1, 3, and 7 for the selected model. Blood volume predictions at the voxel level had a high level of correlation (PCCs greater than 0.65) for both the selected and ensemble averaged model. Lower agreement was observed between the predicted and measured blood volume resulting in CCCs greater than 0.49 for the selected model and 0.36 for the ensemble model.
Cancers 13 01765 g005
Figure 6. Model predictions (scenario 2) for animal 1 and 3. Model predictions of tumor and blood volume fractions at the central slice are shown for animal 1 and animal 3. For both animals, the model predictions from the selected model and the ensemble average are shown. The bottom two rows show the T2-weighted and post-contrast T1-weighted MRI. For animal 1, both the selected and ensemble average model predict the areas of intratumor heterogeneity in the tumor volume predictions. Blood volume predictions tended to predict increased blood volume in the interior of the tumor relative to the periphery as observed. For animal 3, the observed tumor is more homogeneous (compared to animal 1) and the model is able to predict this distribution. Additional, both models predict a higher blood volume fraction towards the interior of the tumor at the final time point which is not present in the measured blood volume fraction. For animal 1, tumor growth was predicted from day 10 to 21 and resulted in −1.0% and 261% error in tumor volume predictions for the selected and ensemble models, respectively. For animal 3, both the selected and ensemble average model resulted in <4.6% error in tumor volume predictions. A sagittal view of this figure is shown in Figure S1.
Figure 6. Model predictions (scenario 2) for animal 1 and 3. Model predictions of tumor and blood volume fractions at the central slice are shown for animal 1 and animal 3. For both animals, the model predictions from the selected model and the ensemble average are shown. The bottom two rows show the T2-weighted and post-contrast T1-weighted MRI. For animal 1, both the selected and ensemble average model predict the areas of intratumor heterogeneity in the tumor volume predictions. Blood volume predictions tended to predict increased blood volume in the interior of the tumor relative to the periphery as observed. For animal 3, the observed tumor is more homogeneous (compared to animal 1) and the model is able to predict this distribution. Additional, both models predict a higher blood volume fraction towards the interior of the tumor at the final time point which is not present in the measured blood volume fraction. For animal 1, tumor growth was predicted from day 10 to 21 and resulted in −1.0% and 261% error in tumor volume predictions for the selected and ensemble models, respectively. For animal 3, both the selected and ensemble average model resulted in <4.6% error in tumor volume predictions. A sagittal view of this figure is shown in Figure S1.
Cancers 13 01765 g006
Figure 7. Model predictions (scenario 2) for animal 5, 6, and 7. Model predictions of tumor and blood volume fractions at the central slice are shown for animal 5, 6, and 7. For all animals, the model predictions from the selected model and the ensemble average are shown. The bottom two rows show the T2-weighted and post-contrast T1-weighted MRI. For animal 5, only one prediction time point was available and both the selected and ensemble average model resulted in <8.7% error in tumor volume predictions. For animal 6, tumor growth was predicted from day 8 to 19 and resulted in 16.0% and 521% error in tumor volume predictions for the selected and ensemble models, respectively. For animal 7, four prediction time points were available and both the selected and ensemble average model resulted in <21.5% in tumor volume predictions across all time points. A sagittal view of this figure is shown in Figure S2.
Figure 7. Model predictions (scenario 2) for animal 5, 6, and 7. Model predictions of tumor and blood volume fractions at the central slice are shown for animal 5, 6, and 7. For all animals, the model predictions from the selected model and the ensemble average are shown. The bottom two rows show the T2-weighted and post-contrast T1-weighted MRI. For animal 5, only one prediction time point was available and both the selected and ensemble average model resulted in <8.7% error in tumor volume predictions. For animal 6, tumor growth was predicted from day 8 to 19 and resulted in 16.0% and 521% error in tumor volume predictions for the selected and ensemble models, respectively. For animal 7, four prediction time points were available and both the selected and ensemble average model resulted in <21.5% in tumor volume predictions across all time points. A sagittal view of this figure is shown in Figure S2.
Cancers 13 01765 g007
Figure 8. Tumor volume prediction time courses (scenario 2). For each animal the measured (blue dots), model mean (black line), model 95% confidence interval (shaded red), and dates of radiotherapy (red dots) are shown for both the selected model (top panel) and ensemble average (bottom panel). For the selected model, larger confidence intervals were observed for animal 1 relative to the other models. For the ensemble average, tumor growth was generally overestimated.
Figure 8. Tumor volume prediction time courses (scenario 2). For each animal the measured (blue dots), model mean (black line), model 95% confidence interval (shaded red), and dates of radiotherapy (red dots) are shown for both the selected model (top panel) and ensemble average (bottom panel). For the selected model, larger confidence intervals were observed for animal 1 relative to the other models. For the ensemble average, tumor growth was generally overestimated.
Cancers 13 01765 g008
Table 1. Model parameters and variables.
Table 1. Model parameters and variables.
Parameter or
Variable Group
Parameter or VariableInterpretationSource
Measured ϕ T ( x ¯ , t ) Tumor cell volume fractionDW-MRI
ϕ V ( x ¯ , t ) Blood volume fractionDCE-MRI
θminMinimum value for carrying capacityDW-MRI
θVMaximum blood volumeDCE-MRI
ϕ V , p r e t r e a t m e n t ¯ Average   pre - treatment   ϕ V ( x ¯ , t ) DCE-MRI
Calibratedkp,T,0Tumor cell proliferation rateCalibrated globally or locally
kd,T,0Tumor cell death rateCalibrated globally
DT,0, DV,0ϕT and ϕV diffusion coefficients in the
absence of mechanically coupling
Calibrated globally
ϕV,threshThreshold on ϕV for carrying capacity
to decrease
Calibrated globally
θmaxMax carrying capacityCalibrated globally
kp,VVasculature proliferation rateCalibrated globally
kd,VVasculature death rateCalibrated globally
SFSurviving fractionCalibrated globally
α1Coupling constantCalibrated globally
AssigneddDistance to the periphery
of the tumor
Calculated
θ T , V ( x ¯ , t ) Combined ϕT and ϕV
carrying capacity
Calculated
λ1Coupling ConstantSet to 0.25
λ2Coupling ConstantSet to 1
G, νShear modulus and Poisson’s ratio
assigned for white and gray matter
Literature [30]
Table 2. Experimental timeline.
Table 2. Experimental timeline.
Animal
Number
Dose/DayImage
Days
Treatment
Dates
Calibration
Dates
Prediction
Dates
812 Gy0,3,5,7,10,
12,17,19,21
3–7,10–140–710–21
24 Gy0,3,53–50–5NA
32 Gy0,3,5,7,10,123–7,10–120–57–12
42 Gy0,3,53–60–5NA
52 Gy0,3,5,73–60–57
64 Gy0,2,5,7,8,
12,14,16,19
0–2,6–9,12,130–78–19
74 Gy0,3,5,7,10,
12,14,17
3–7,10–140–710–17
Table 3. Calibrated model parameters for selected model.
Table 3. Calibrated model parameters for selected model.
Parameter or
Variable
Animal Number (Dose per Day)
1 (2 Gy)2 (4 Gy)3 (2 Gy)4 (2 Gy)5 (2 Gy)6 (4 Gy)7 (4 Gy)
kp,T,0 (day−1)0.84 ± 0.251.20 ± 0.351.40 ± 0.481.10 ± 0.161.82 ± 0.214.65 ± 5.452.04 ± 0.43
kd,T,0 (day−1)0.22 ± 0.040.40 ± 0.070.11± 0.010.03 ± 0.060.5 ± 0.100.35 ± 0.04NA
DT,0
(104 m2 day−1)
2.32 ± 0.074.50 ± 0.072.58 ± 0.071.37 ± 0.023.44 ± 0.056.93 ± 0.303.38 ± 0.02
DV,0
(104 m2 day−1)
2.72 ± 0.063.00 ± 0.093.10 ± 0.093.00 ± 0.193.11 ± 0.053.25 ± 0.123.00 ± 0.04
ϕV,thresh0.03 ± 0.120.05 ± 0.060.03 ± 0.050.03 ± 0.020.02 ± 0.040.02 ± 0.020.02 ± 0.08
θmax0.90 ± 0.010.90 ± 0.020.94 ± 0.020.87 ± 0.030.92 ± 0.010.91 ± 0.000.92 ± 0.01
kp,V (day−1)0.65 ± 0.031.10 ± 0.060.82 ± 0.041.25 ± 0.171.49 ± 0.130.17 ± 0.471.76 ± 0.02
kd,V (day−1)0.18 ± 0.020.32 ± 0.060.29 ± 0.040.54 ± 0.100.97 ± 0.170.07 ± 0.040.07 ± 0.04
SF0.95 ± 0.010.80 ± 0.080.97 ± 0.010.99 ± 0.010.98 ± 0.020.90 ± 0.010.91 ± 0.01
α1NANA4.00 ± 0.17NANANANA
The mean 95% confidence interval are reported for each parameter. NA indicated that the parameter was not estimated for the selected model.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Hormuth, D.A., II; Jarrett, A.M.; Davis, T.; Yankeelov, T.E. Towards an Image-Informed Mathematical Model of In Vivo Response to Fractionated Radiation Therapy. Cancers 2021, 13, 1765. https://doi.org/10.3390/cancers13081765

AMA Style

Hormuth DA II, Jarrett AM, Davis T, Yankeelov TE. Towards an Image-Informed Mathematical Model of In Vivo Response to Fractionated Radiation Therapy. Cancers. 2021; 13(8):1765. https://doi.org/10.3390/cancers13081765

Chicago/Turabian Style

Hormuth, David A., II, Angela M. Jarrett, Tessa Davis, and Thomas E. Yankeelov. 2021. "Towards an Image-Informed Mathematical Model of In Vivo Response to Fractionated Radiation Therapy" Cancers 13, no. 8: 1765. https://doi.org/10.3390/cancers13081765

APA Style

Hormuth, D. A., II, Jarrett, A. M., Davis, T., & Yankeelov, T. E. (2021). Towards an Image-Informed Mathematical Model of In Vivo Response to Fractionated Radiation Therapy. Cancers, 13(8), 1765. https://doi.org/10.3390/cancers13081765

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop