PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Phys Biol. Author manuscript; available in PMC 2016 July 1.
Published in final edited form as:
PMCID: PMC4486062
NIHMSID: NIHMS699087

Predicting in vivo glioma growth with the reaction diffusion equation constrained by quantitative magnetic resonance imaging data

Abstract

Reaction-diffusion models have been widely used to model glioma growth. However, it has not been shown how accurate this model can predict future tumor status using model parameters (i.e., tumor cell diffusion and proliferation) estimated from quantitative in vivo imaging data. Towards this end, we used in silico studies to develop the methods needed to accurately estimate tumor specific reaction-diffusion model parameters, and then tested the accuracy with which these parameters can predict future growth. The analogous study was then performed in a murine model of glioma growth. The parameter estimation approach was tested using an in silico tumor “grown” for ten days as dictated by the reaction-diffusion equation. Parameters were estimated from early time points and used to predict subsequent growth. Prediction accuracy was assessed at global (total volume and Dice value) and local (concordance correlation coefficient, CCC) levels. Guided by the in silico study, rats (n=9) with C6 gliomas, imaged with diffusion weighted magnetic resonance imaging (DW-MRI), were used to evaluate the model’s accuracy for predicting in vivo tumor growth. The in silico study resulted in low global (tumor volume error < 8.8 %, Dice > 0.92) and local (CCC values > 0.80) level errors for predictions up to six days into the future. The in vivo study showed higher global (tumor volume error > 11.7%, Dice < 0.81) and higher local (CCC < 0.33) level errors over the same time period. The in silico study shows that model parameters can be accurately estimated and used to accurately predict future tumor growth at both the global and local scale. However, the poor predictive accuracy in the experimental study suggests the reaction-diffusion equation is an incomplete description of in vivo C6 glioma biology and may require further modeling of intra-tumor interactions including segmentation of (for example) proliferative and necrotic regions.

1. Introduction

Mathematical models have been constructed to describe tumor growth and invasion over a large range of spatial scales (nm to cm) and temporal scales (ns to years). Substantial discussions have focused on translating these models to clinical care with the long term goal of providing clinicians with patient-specific predictions of future tumor growth and therapy response in order to optimally select and guide patient therapy [13]. Approaches for patient-specific predictions may focus on changes in a single property such as tumor volume, or changes in tumor growth as a function of several related properties (e.g., cellularity, vascularity, nutrient distribution). Models that focus on the change in a single tumor property can be parameterized readily with experimental data [4,5], but may fail to capture spatial and temporal tumor heterogeneity of, for example, cellularity, vasculature density, proliferation rates, and the level of response (or lack thereof) of cells to treatment that is observed within tumors [6,7]. Patient-specific models that capture a tumor’s spatial and temporal heterogeneity could be used to more accurately describe the delivery of treatment and subsequent response [811]. Unfortunately, modeling these characteristics frequently requires knowledge of parameters that can only be measured by highly invasive methods or within idealized (in vitro) settings [1215]. The reliance of the existing modeling literature on parameters that are either extraordinarily difficult or impossible to measure non-invasively fundamentally limits their clinical application. Recasting these models in terms of parameters measured via non-invasive imaging measurements would dramatically improve the clinical relevance of patient-specific tumor growth predictions [1].

Magnetic resonance imaging (MRI) and positron emission tomography (PET) can be used to provide an array of non-invasive, quantitative, and functional measurements in 3D and at multiple time points of tumor growth. More specifically, MRI and PET can provide measurements of cellularity [16], blood volume [17,18], blood flow [17,18], hypoxia [19], oxygen saturation [20], and metabolism [21]. Additionally, the ability to make repeatable, non-invasive, spatially discretized, quantitative measurements of tumor growth supports the development, testing, and refinement of mathematical descriptions of in vivo tumor growth. Several groups [1,5,2229] have incorporated imaging measurements from MRI, PET, and x-ray computed tomography into mathematical models of tumor growth. Preliminary efforts in both breast [23] and pancreatic [24] cancers, have shown that patient specific imaging data can potentially accurately predict future tumor growth. This, however, has not been demonstrated for gliomas.

One common model for glioma growth is the reaction-diffusion model, whereby the spatio-temporal change in tumor cellularity is due to proliferation and invasion (described by random diffusion) of tumor cells. The proliferation and invasion of cells are typically characterized with a proliferation rate and a diffusion coefficient, respectively. The reaction-diffusion model of glioma growth described by Swanson et al [29], uses proliferation and diffusion coefficients of tumor cells estimated from T2-weighted and post-contrast T1-weighted MRI data obtained at two time points. The estimated tumor cell proliferation and diffusion values can then be used to simulate tumor growth following surgical resection [29] or simulate a virtual control to assess patient response to radiotherapy [25]. Jbadbi et al [28] extended this approach by allowing anisotropic diffusion of tumor cells by replacing the diffusion coefficient with a diffusion tensor measured using diffusion tensor imaging (DTI). The authors showed that simulated anisotropic tumor growth better matched the shape of glioma growth observed in patients. Spatially varying estimates of diffusion and proliferation were included in the work of Ellingson et al [27]. In this work, serial diffusion weighted MR images were used to develop a voxel-wise analytical solution (when certain assumptions are satisfied) to a reaction-diffusion model of glioma growth. The proliferation and diffusion values were compared to MR spectroscopy measurements, but these values were not used to simulate tumor growth. Another extension of the reaction-diffusion model is the incorporation of mechanical properties of healthy and tumor tissue into a description of tumor growth [23,24,26,30]. The work of Hogea et al [26] showed the benefit of incorporating mechanical deformations caused by the invading tumor growth into a reaction-diffusion model. Their effort also demonstrated the means to invert their model system to estimate parameters from imaging data.

In this work we use a reaction-diffusion model of glioma growth with proliferation and diffusion values estimated from quantitative in vivo imaging data to predict future tumor growth and then validate (or refute) that prediction by direct comparison to future in vivo measurements. Using an in silico tumor we first developed the means to accurately estimate model parameters and assessed the accuracy of tumor growth predictions. We then performed the analogous in vivo study, where model parameters were estimated from serial diffusion-weighted MRI data in a murine model of glioma to predict future tumor status which could then be directly compared to experimental outcome. The in silico experiments show that model parameters can be accurately estimated from tumor growth datasets and then used to predict future tumor growth with low global and local errors. However, when the approach is applied to in vivo glioma measurements, it is shown that the reaction-diffusion model provides poor predictive ability of future tumor growth.

2. Materials and Methods

2.1. Modeling Approach

The reaction-diffusion equation describes the spatio-temporal rate of change in tumor cell number and distribution due to the random movement of tumor cells (diffusion; the first term on the right hand side of Eq. (1)), and proliferation (reaction; the second term on the right hand side of Eq. (1)):

N(x¯,t)t=(D(x¯)N(x¯,t))+k(x¯)N(x¯,t)(1-N(x¯,t)θ),
(1)

where N ([x with macron],t) is the number of tumor cells at three-dimensional position [x with macron] and time t, D([x with macron]) is the tumor cell diffusion coefficient at position [x with macron], k ([x with macron]) is the net tumor cell proliferation at position [x with macron], and θ is the tumor cell carrying capacity. Note that the proliferation term varies temporally as a function of cell density, N ([x with macron],t), although it is assumed that the proliferation rate, k ([x with macron]), is temporally constant. As described below, quantitative DW-MRI provides estimates of N ([x with macron],t). These data are obtained at multiple time points, early in the tumor’s life cycle, and used to solve an inverse problem using Eq. (1) to return estimates of k ([x with macron]) for each voxel within the tumor, and two D([x with macron]) values: one for white matter (Dwm) and one for gray matter (Dgm). The forward evaluation of Eq. (1) is solved using a three dimension in space, fully explicit finite difference (FD) in time simulation written in Matlab (Mathworks, Natick, MA). The simulation domain has no diffusive flux of tumor cells at brain tissue boundaries (i.e., at the skull) and the grid spacing matches the spatial resolution of the MRI data (Δx = 250 μm, Δy = 250 μm, Δz = 1000 μm). The simulation time step was set at 0.01 days.

2.2. In Silico Experiments

Figures 1 and and22 show the approach for the in silico experiments. An initial distribution of tumor cells, N ([x with macron],t0), was seeded within a rat brain domain. A spatial map of k was determined from DW-MRI estimates of N ([x with macron],t) from a C6 glioma bearing rat using Eq. (2):

k(x¯)=-log(N(x¯,t1)/N(x¯,t0))/(t1-t0),
(2)

where N ([x with macron],t0) and N ([x with macron], t1) represent the distribution of cells at time t0 and t1, respectively, while k ([x with macron]) remained constant in time. The initial model parameters, or Ptrue (P = (k(1), k(2),… k(n), Dgm, Dwm), where n is the number of voxels within the tumor), were selected by iteratively scaling k ([x with macron]) until a 12 fold increase in total tumor volume was observed over 10 days, matching the average observed tumor volume increase observed in the in vivo study. An FD simulation of Eq. (1) with an initial distribution of tumor cells N ([x with macron],t0) and parameters Ptrue was used to grow an in silico tumor for 10 days (1000 time steps/iterations). Tumor cell distributions, N ([x with macron],t), were then sampled at days 0, 2, and 4 thru 10. Panel (a) in Figure 1 shows a central anatomical axial slice through a rat head and cropped images of the in silico tumor cell distribution seeded at day 0 and “grown” to days 2 and 4. Three different combinations of these three time points were then used to estimate three sets of model parameters (P1, P2, and P3). P1 was estimated using tumor cell measurements from days 0 and 4 (N ([x with macron],t0) and N ([x with macron],t4), respectively). P2 was estimated using tumor cell measurements from days 2 and 4 (N ([x with macron],t2) and N ([x with macron],t4), respectively), while P3 was estimated using tumor cell measurements from days 0, 2, and 4 (N ([x with macron],t0), N ([x with macron],t2), and N ([x with macron],t4), respectively). Panel (b) in Figure 1 shows the parameter optimization approach for these model parameters (P1, P2, and P3). For each set of estimated parameters (P1, P2, and P3) a FD simulation of Eq. (1) was initialized with N ([x with macron],t0) (for P1 and P3) or N ([x with macron],t2) (for P2) and used to grow a tumor to day 4 resulting in a model estimate of N for each parameter set (Nest,1, Nest,2, and Nest,3). The error between N and Nest,1, Nest,2, and Nest,3, respectively, was calculated. If this error is minimized, the optimal values of P1, P2, and P3 have been determined, otherwise P1, P2, and P3 are updated with new values. The optimized P values are then used to predict future tumor growth.

Figure 1
Parameter optimization approach
Figure 2
Predicted tumor growth modeling approach

Figure 2 shows the modeling approach for predicting future tumor cell distributions. For each set of optimized parameters (P1, P2, and P3) an FD simulation of Eq. (1) was used to “grow” the tumor from day 4 to day 10 (6 days; 600 iterations) resulting in predicted tumor cell distributions for each parameter set (Npred,1 ([x with macron],t5–10), N pred,2 ([x with macron],t5–10), and N pred,3 ([x with macron],t5–10), respectively). Error between the true N and Npred was calculated at both the global and local levels by calculating the percent error in tumor volume, the Dice similarity coefficient, the normalized root mean square error (nRMS error), and the concordance correlation coefficient (CCC).

2.3. In Vivo Experiments

All experimental procedures were approved by Vanderbilt University’s Institutional Animal Care and Use Committee. Female Wistar rats (n = 9, 236–263 g) were anesthetized, given analgesics, and inoculated with C6 glioma cells (1 × 105) via stereotaxic injection. During each MRI procedure body temperature was maintained near 37° C by a flow of warm air directed over the animal and respiration was monitored using a pneumatic pillow. Each rat was anesthetized using 2% isoflurane in 98% oxygen for all surgical and imaging procedures. Rats were imaged beginning 10 days post-surgery (defined as day 0). Rats were imaged up to 10 days after the first imaging time point. The first three imaging measurements for all rats occurred on days 0, 2, and 4. Rats 1–3 were then imaged on days 5, 8, and 10. Rats 4–5 were imaged on days 5, 6, and 9. Rat 6 was imaged on days 5, 6, 8 and 10, while rats 7–8 were imaged only on days 5 and 6. Rat 9 was only imaged at one additional time on day 5.

MRI was performed on a 9.4 T horizontal-bore magnet (Agilent, Santa Clara, CA, USA). The animal’s head was positioned in a 38 mm diameter Litz quadrature coil (Doty Scientific, Columbia, SC, USA) and was secured by a bite bar. All MR images were sampled with a 128 × 128 × 16 matrix acquired over a 32 × 32 × 16 mm3 field of view. In order to facilitate the modeling, the imaging volumes obtained at time points two through the end of the experiment were registered to the first time point via a mutual information based rigid registration algorithm performed at the scanner [31]; this ensures that the image volumes obtained at each time point are very nearly identical (See Supplementary Figures 1–3 for example registration results). T1 map was produced using data from an inversion-recovery snapshot experiment with TR/TE = 5000/3 ms, TI (inversion time) = (8 TIs logarithmically spaced between 200 – 4000 ms), and two averaged excitations.

DW-MRI was acquired using a pulsed fast spin echo diffusion sequence in three orthogonal diffusion encoding directions with b-values of 0, 300, 500, 700, 900, and 1100 s/mm2, and Δ/δ = 25 ms/2 ms. The apparent diffusion coefficient (ADC) was estimated on a voxel basis using a two parameter fit of the DW-MRI data [32]. To determine N ([x with macron],t), the ADC values from the DW-MRI data are then transformed to estimate cell number [33,34] using Eq. (3):

N(x¯,t)=θ(ADCw-ADC(x¯,t)ADCw-ADCmin),
(3)

where θ represents the tumor cell carrying capacity, ADCw is the ADC of free water at 37° C (2.5 × 10−3 mm2/s) [32], ADC([x with macron],t) is the ADC value at position [x with macron] and time t, and ADCmin is the minimum ADC value which corresponds to the voxel with the largest number of cells. The carrying capacity, θ, was calculated for each imaging voxel assuming spherical tumor cells with a packing density of 0.7405 [35] with an average cell volume of 908 μm3 [36].

Tumor regions-of-interest (ROI) were manually placed at each time point using the T1 maps. ADC measurements within these ROI’s were then transformed to tumor cell number using Eq. (3). T1 maps were used to define tumor, white, and gray matter regions in the MR images.

Similar to the in silico experiments, three sets of model parameters (P1, P2, and P3) were then estimated for each rat. We note that the time and spatial origin of the tumor (as mentioned in Hogea et al [26]) was not determined as tumor growth simulations were initialized with tumor cellularity measurements from DW-MRI. The estimated parameters for each rat (parameters for rat 1; PR1,1, PR1,2, and PR1,3) were used in a FD simulation of Eq. (1) to “grow” simulated tumors from day 4 to day 10 (6 days; 600 iterations) resulting in predicted tumor cell distributions for each parameter set (predicted N ([x with macron],t) for rat 1; NR1, pred,1 ([x with macron],t5–10), NR1, pred,2 ([x with macron],t5–10),and NR1, pred,3 ([x with macron],t5–10), and, respectively).

The three different time point combinations from the in vivo data sets were also fit to a model substituting a spatially invariant proliferation rate (kROI) for the spatially variant proliferation rate (k([x with macron])); i.e., k([x with macron]) [equivalent] kROI in Eq. (1).

2.4. Numerical Methods

A Levenberg-Marquardt weighted least squares nonlinear optimization, implemented with a regularization parameter described in Joachimowicz et al [37,38], was used to estimate model parameters (P1, P2, P3) from tumor cell distribution measurements. All parameters were constrained to non-negative values. Prior to estimating P a 3 × 3 Gaussian filter was applied to each slice of N ([x with macron],t) to reduce the effects of noise within individual voxels. During the optimization scheme, k was estimated voxel-wise in areas within the tumor ROI and assigned 0 elsewhere. Additionally, Dwm and Dgm values were assigned region-wise using a white and gray matter map. The optimized parameters were determined when the objective function, Eq. (4), was minimized:

t=titf((x¯=1x¯i=n(N(x¯,t)))-1·(x¯=1x¯=n(Nest(x¯,t)-N(x¯,t))2)),
(4)

where ti is the initial time point, tf is the final time point, n is the total number of voxels within the tumor, and Nest ([x with macron],t) is the model estimate of N using the current parameter set. For P1 and P2 optimizations, ti and tf were equal to day 4. For P3 optimization ti was equal to day 2 and tf was equal to day 4.

To assess the effect of noise in ADC measurements on estimates of P, in silico parameter optimization was repeated (N = 100) for each set of parameters (P1, P2, P3) with noise added to N ([x with macron],t) from a normal distribution with a zero mean and a standard deviation of 3.3% of the carrying capacity (selected based on the reproducibility of ADC measurements in vivo [32]).

After optimization of P, these values were used in a FD implementation of Eq. (1), initialized with the tumor cell distribution at day 4 (N ([x with macron],t4)), to “grow” a predicted tumor from day 4 thru 10 (600 iterations). Throughout the FD simulation, as the tumor expanded into regions where an estimate of k was unavailable, k was assigned using a local average of available non-zero k’s within a 3 × 3 × 3 kernel.

The accuracy of P estimated from the in silico dataset was evaluated by computing the percent error between the true, Ptrue, and the estimated parameter sets (P1, P2, and P3), the Pearson correlation coefficient (PCC), and the CCC (similar to the PCC but with a penalty for data that do not lie on the line of unity) [39]. A Bland-Altman analysis was also performed between the true, Ptrue, and the estimated (P1, P2, and P3) parameter sets. For the in vivo study, agreement between P1, P2, and P3 estimates of k([x with macron]) was assessed by calculating the PCC and CCC for P1 and P2, P1 and P3, and P2 and P3. For both the in vivo and in silico studies error between the predicted and true tumor growth was assessed at the global (i.e., volume) and local (i.e., voxel) levels. For the in vivo analysis, the tumor ROIs and measured cellularity from DW-MRI are taken as true tumor growth. At the global level, error was assessed by calculating the percent error in tumor volume, the nRMS error, and the Dice similarity coefficient (a measure of spatial overlap between two data sets ranging from 0 (no overlap) to 1 (complete overlap); [40]). The percent error in tumor volume and nRMS error were computed by comparing the true tumor volume and the tumor volumes predicted from FD simulations using P1, P2, and P3 at days 5 thru 10. The Dice similarity coefficient was computed by comparing the spatial overlap between the true tumor ROIs and the tumor ROIs predicted from FD simulations using P1, P2, and P3 at days 5 thru 10. At the local level, error was assessed by computing the CCC between N and Npred at days 5 thru 10. For the in vivo study, paired t-tests were used to evaluate the differences between percent error in tumor volume, Dice, and CCC results observed with the spatially variant k and spatially invariant k models.

3. Results

3.1. In Silico Results

Illustrative results of the in silico experiments are shown Figures 35 and summarized in Table 1. Figure 3 shows the true distribution of k (panel (a)) used to “grow” the in silico tumor, the estimated values of k (panels (b – d)), and plots of the true voxel values of k against the estimated values of k (panels (e – g)). Panels (b) and (e) show the results from parameters estimated from days 0 and 4. Both a low level of agreement (CCC = 0.38) and a weak linear relationship (PCC = 0.54) is observed between the true k and the k estimated using the P1 data sets (kP1). Additionally 88% of voxels in kP1 are overestimated. Parameters estimated using days 2 and 4 (panels (c) and (f)) showed an improved level of agreement (CCC = 0.74), stronger linear relationship (PCC = 0.80), and fewer overestimated voxels (72%) compared to kP1. Using all three time points (panels (d) and (g)) resulted in the best level of agreement (CCC = 0.84), the strongest linear relationship (PCC = 0.87), and the fewest overestimated voxels (58%). Panels (h–j) show the Bland-Altman analysis for P1, P2, and P3 compared to Ptrue. The black lines represent the mean difference while the gray lines represent the 95% confidence interval (CI). A larger 95% CI was observed for kP1kPtrue (−0.35 to 0.91) compared to kP2kPtrue (−0.25 to 0.46) and kP3kPtrue (−0.26 to 0.32).

Figure 3
True and estimated proliferation rate maps from in silico study
Figure 5
Global and local level error analysis for in silico study
Table 1
Parameter estimation error from in silico experiments

Figure 4 shows the true (panel (a)) and the predicted tumor cell distributions (top rows in panels (b) – (d)) and the percent difference between the true and predicted distributions (bottom rows in panels (b) – (d)). Panel (b) shows the predicted N at day 5 using P1, P2, and P3. The highest error between N and Npred at day 5 was observed for parameters P1 (mean ± standard error; 31.15 ± 0.82%) compared to P2 (19.67 ± 0.64%) and P3 (16.84 ± 0.58%). As the tumor continues to grow, increased error is observed between N and Npred. Generally, this error is increased at the periphery (greater than 100% error) relative to the interior (less than 40%). The predicted N at day 8 (panel (c)) resulted in increased mean error relative to day 5 for predictions using P1 (41.38 ± 0.87%), P2 (28.03 ± 0.68%) and P3 (26.03 ± 0.64%). At the final time point (panel (d)), P3 based predictions had a mean error (29.51 ± 0.65%) lower than P1 (45.93 ± 0.90%) and P2 (30.83 ± 0.67%) based predictions. The lowest cumulative error was observed for P3 based predictions (nRMS error; mean = 0.062, standard error = 1.33 × 10−2) compared to P1 based predictions (mean = 0.289, standard error = 2.51 × 10−2) or P2 based predictions (mean = 0.102, standard error = 1.19 × 10−2).

Figure 4
True and predicted tumor cell distributions for in silico study

The results of the ROI and voxel level analysis are shown in Figure 5. Error in tumor volume generally increases the further out in time a prediction is made (panel (a)). Percent error in tumor volume ranged from 11.9 – 36.4% for P1 based predictions, 3.1 – 13.6% for P2 based predictions, and 0.8 – 8.8% for P3 based predictions. All parameter sets, however, resulted in Dice values (panel (b)) greater than 0.83 at days 5 thru 10. An increased level of agreement was observed at the voxel level (panel (c)) for P3 based predictions (CCC: 0.80 – 0.99) relative to P1 based predictions (CCC: 0.53 – 0.93) and P2 based predictions (CCC: 0.74 – 0.98).

Table 1 shows the mean percent error between the true values and estimated values of k, Dwm, and Dgm. Using all three time points resulted in less than 1.0% error in estimates of Dwm, whereas using two time points (days 0 and 4 or days 2 and 4) resulted in greater than 49.5% error. Less than 6.2% error was observed in estimates of Dgm when using days 2 and 4 or days 0, 2, and 4, while 13.7% error was observed when using days 0 and 4 to estimate model parameters. Similarly, the highest mean error in k was observed for parameters estimated from days 0 and 4 (28.3 ± 2.9%), while lower error was observed for the approaches using days 2 and 4 (13.1 ± 2.1%) and days 0, 2, and 4 (5.9 ± 1.8%).

3.2. In Vivo Results

The results of the in vivo experiments are shown in Figures 68 and Tables 23. Figure 6 shows the PCC (a) and CCC (b) analysis for all nine rats, as well as kP1, kP2, and kP3 from rats 3 (c–e) and 6 (f–h). The green bars in panels (a–b) represent the comparison of kP1 to kP2, while the blue and red bars represent the kP1 to kP3 and kP2 to kP3 comparisons, respectively. A high level of correlation existed between kP1 and kP3 (mean PCC = 0.75, standard error = 0.05) and between kP1 and kP2 (mean PCC = 0.72, standard error = 0.09) compared to kP2 to kP3 (mean PCC = 0.46, standard error = 0.08). Similar comments apply to the CCC trends. Panels (c–h) show estimated kP1 (panels (c) and (f)), kP2 (panels (d) and (g)), and kP3 (panels (e) and (h)) for rats 3 and 6. Rat 3 (c–e) is an example of a rat with a high level of correlation (PCC: 0.68 to 0.88) but a low level of agreement (CCC: 0.08 to 0.34). Rat 6 (f–h), however, is an example of a rat with both a high level of correlation (PCC: 0.72 to 0.93), and a high level of agreement (CCC: 0.60 to 0.84).

Figure 6
Estimated proliferation rate maps from in vivo study
Figure 8
Global and local level error analysis for in vivo study
Table 2
Average diffusion parameter values and nRMS error from in vivo experiments
Table 3
Average k values from in vivo experiments

Figure 7 shows the true (panel (a)) and the predicted tumor cell distributions (top rows in panels (b – d)) and the percent difference between the true and predicted distributions (bottom rows in panels (b – d)) for rat 1. The left column in panel (a) shows T2-weighted images with a white box indicating the simulation domain and a black outline around the tumor. Panel (b) shows estimated k and the predicted N at day 5 using P1, P2, and P3. Decreased kP1 (mean ± standard error; 0.08 ± 0.02 day−1), kP2 (0.13 ± 0.03 day−1), and kP3 (mean ± standard error; 0.08± 0.02 day−1) were observed for the low cell density regions relative to the rest of the tumor (kP1; 1.51 ± 0.06 day−1, kP2; 1.87 ± 0.10 day−1, kP3; 1.32 ± 0.05 day−1). The highest error between N and Npred at day 5 was observed for parameters P3 (mean ± standard error; 17.80 ± 0.50%) compared to P1 (17.07 ± 0.50%) and P2 (17.12± 0.51%). Increased error (greater than or equal to 100% error) is observed between N and Npred at both the periphery and in regions where N has low cell numbers (less than 50% of a voxels carrying capacity). This pattern is observed also in panels (c) and (d). Increased error at day 8 (panel (c)) was observed relative to day 5 (panel (b)) with P1 predictions having the highest error (32.58 ± 0.63%) compared to P2 (28.27 ± 0.57%) and P3 (31.43 ± 0.61%). The predicted N at day 10 (panel (d)) shows an overestimation of tumor size for P2 based predictions compared to P1 or P3 based predictions. However, at the voxel level the highest mean error was observed for P1 based predictions (35.80 ± 0.54%) relative to P2 (35.03 ± 0.53%) and P3 (33.45 ± 0.51 %) based predictions.

Figure 7
True and predicted tumor cell distributions for rat 1

Figure 8 presents the global and local level error analysis for the in vivo experiments. Panels (a–c) show the results when a spatially variant k (i.e., k ([x with macron])) is estimated and panels (d–f) show the results for the spatially invariant estimated k (i.e., kROI). For the spatially variant k, percent error in tumor volume (panel (a)) for P1, P2, and P3 based predictions ranged from 14 – 34%, 16 – 50%, and 12 – 29%, respectively. Lower Dice values (panel (b)) were observed compared to the in silico study and ranged from 0.67–0.81 for all approaches. Similarly, the in vivo study had decreased level of agreement (panel (c)) compared to the in silico study with CCC’s less than 0.33 for all approaches. For the spatially invariant k, percent error in tumor volume (panel (d)) for P1, P2, and P3 based predictions ranged from 36 – 58%, 36 – 77%, and 32 – 54%, respectively. The Dice values were also lower than the in silico study and ranged from 0.66–0.79. Lower agreement (CCC < 0.25) at the voxel level was also observed for the spatially invariant k approach compared to the spatially variant k approach. The spatially invariant k’s percent error in tumor volume was significantly greater (P<0.05) than the spatially variant k’s results for all parameter sets. Similarly, both the Dice and CCC values were significantly smaller (P<0.05) for the spatially invariant k predictions compared to the spatially variant k predictions.

The average diffusion estimates and the nRMS error for both the spatially variant and invariant k models are shown in Table 2. For the spatially variant k model fits, the mean Dwm ranged from 1.49 × 104 to 1.57 × 104 μm2/day, and the mean Dgm per animal ranged from 1.58 × 104 to 1.99 × 104 μm2/day. The cumulative error (nRMS error) was lowest for P3 based predictions (0.49) compared to P1 (0.54) or P2 (0.81) based predictions. The spatially invariant k model optimization resulted in higher diffusion values for both Dwm (greater than 3.05 × 104) and Dgm (greater than 3.50 × 104) compared to the spatially variant k results. Increased cumulative error (nRMS error > 1.04) was also observed compared to the spatially variant k results.

Table 3 shows the average proliferation rate for each rat from the spatially variant (k ([x with macron]), where k(x¯)¯ is the average voxel-wise k estimated within the tumor) and the estimated spatially invariant proliferation rate (kROI). k(x¯)¯ ranged from 0.51 to 4.06 day−1, while kROI ranged from 0.94 to 9.94 day−1. kROI was larger than k(x¯)¯ for six rats for parameter estimates using day 0 and day 4, nine rats for days 2 and 4, and three rats when all three time points were used.

4. Discussion

The results of the in silico experiments indicate that the parameters within the reaction-diffusion equation (i.e., Dwm, Dgm, and k) can be accurately estimated and then used to accurately predict future tumor growth at the local and global levels, provided the tumor’s growth is described by the reaction-diffusion equation. While parameters estimated from data with experimentally observed noise does increase the error between the true and estimated values, when three time points are used the error between the true and observed parameters is less than 5.86%. Thus, the addition of a third time point decreases the sensitivity of the parameter optimization algorithm to approximately the same level of noise that is present in the measurement. The increased error observed for Dwm relative to Dgm may be explained due to only 4% of the voxels in the domain being identified as white matter. The limited number of voxels containing white matter most likely makes the model less sensitive to changes in Dwm compared to Dgm.

The results of the in silico study show, the strength of both the estimated parameters and the forward evaluation algorithm, as exhibited in the high level of overlap of tumor volumes (Dice values greater than 0.83 for P1, P2, and P3) and strong agreement at the local level (CCC values greater than 0.80 for P3; greater than 0.53 for P1 and P2) between N and Npred (Figure 5). The largest disagreements occur at the tumor edges. The simple propagation (local average) of proliferation rates outside of the parameter estimation region propagates errors and may need to be improved to incorporate additional information (e.g., local cellularity, distance from vasculature, and nutrient concentration; importantly such data is also available from clinically relevant, non-invasive imaging studies [41]). The significant differences in prediction errors (globally and locally) between the two time point approaches (P1 and P2) suggest the parameter optimization approach is sensitive to the spacing between measurements (P1: 4 days, P2: 2 days). Additionally, the lower error in P2 based predictions compared to P1 based predictions suggest that the tumor growth between days 2 and 4 is more representative of future growth than growth between days 0 and 4.

The in vivo experiments demonstrated greater error at both the global and local level compared to the in silico experiments. The increased error suggests that the reaction-diffusion equation is an incomplete description of C6 biology. The overestimation of tumor volume estimates suggests that overall tumor growth properties are changing between the estimation time points and the prediction time points. At the global level, the expansion of the tumor may be less restricted at earlier time points compared to later time points. At the local level, these changes may be the result of an increase or decrease in proliferation due to changes in the viability of the cells within a particular voxel. The very poor CCCs (less than 0.33, Figure 8) similarly suggest that the reaction-diffusion equation provides a poor description of local properties. The reaction-diffusion equation does, however, provide tumor growth predictions that co-localize (Dice values greater than 0.62) with the true tumor volumes. Different from the in silico study, P1 predictions had lower percent error in tumor volume (less than 34.4%) compared to P2 based predictions (less than 50.1%). This suggests that the tumor growth over days 0 and 4 are more representative of future in vivo growth than the tumor growth over days 2 and 4. The larger distance between measurements allows potentially inconsistent growth rates between days 0 and 2 and days 2 and 4 to be averaged over 4 days, lessening the effects of non-representative volumetric growth on model estimates and predictions. Replacing the voxel-specific proliferation, k ([x with macron]), with a tumor-specific proliferation rate, kROI, resulted in larger global (increased percent error in tumor volume, decreased Dice values) and local (decreased CCCs) level errors. The decreased agreement between predicted and observed tumor ROIs (decreased Dice values) using kROI suggest that a spatially variant k is an important factor in predicting tumor geometry by allowing variations in regional tumor expansion. This factor may also contribute to the increased tumor volume error due to a more uniform expansion of tumor growth.

There are several limitations in this current approach. One limitation is the assumption that all tumor cells within a given voxel (spatially variant k) or within the tumor (spatially invariant k) follow the same proliferation rules. Within tumors there may be groups of actively proliferating cells as well as cells that are quiescent or necrotic [42]. In particular, necrotic tissues (which can be relatively large compared to total tumor volume) can strongly influence tumor growth [42] and patient prognosis [43,44]. Models incorporating different proliferation rules have the potential to more accurately describe in vivo proliferation [45,46]; however, it is challenging to initialize these models using non-invasive measurements. Although the proliferation model in this approach is limited, the spatially variant k lessens the potential error (relative to kROI) in this assumption by discretizing the tumor into individual regions that can have a proliferation rate that more closely captures local behavior. A second limitation is that proliferation rates, k, estimated from early time points are assumed to be constant for the remainder of the tumor growth. The logistic growth term in Eq. (1) allows for temporally variation of the instantaneous growth rate as cell density increases or decreases, but it does not allow for the individual voxel proliferation rates to vary temporally over the course of the experiment; this is a fundamental limitation of this model as formulated. As in vivo tumors expand, the characteristics of a cell’s environment will change and either accelerate or slow future cell proliferation. A more realistic model would incorporate these phenomena [12,13] and adjust the spatio-temporal distribution of k throughout the simulation. These approaches, however, require model parameters that are extraordinarily difficult to measure non-invasively thereby limiting their application to subject-specific model predictions. This approach also does not consider the impact that tumor necrosis and edema may have on tumor cell proliferation [42,47,48] and the estimation of cellularity from ADC measurements (e.g., increased water diffusion may be observed in necrosis or necrosis adjacent regions due to breakdown of barriers to free water movement). However, there currently is not a widely accepted or validated method for segmenting necrosis and edema a priori in brain tumors. A third limitation is that as the tumor expands into voxels not included in the parameter estimation procedure (i.e., voxels outside of the tumor ROI determined at day 4), the proliferation rate in that voxel is then assigned as the average of the nearby known proliferation rates. This average value of the local k does not account for differences in environmental conditions, cell distributions, or cell phenotypes that may alter a voxel’s k. A fourth limitation of Eq. (1) is that D is temporally fixed which results in tumor growth that is unrestricted (i.e., the model assumes that the tumor is growing into an empty space and does not incorporate the effects of the surrounding tissue [49]) and unresponsive to changes in microenvironment properties (e.g., extracellular matrix components, growth factors, tumor necrosis factor, matrix-metalloproteinase [50,51]). This limitation can be amended through including more realistic terms such as mass effect [26,52] to temporally adjust tumor migration behavior which may increase the predictive accuracy of the model. We do note, though, that the models which include the spatial-temporal evolution of D and k must carefully consider how the model parameters will be initialized (i.e., how are the values for the additional model parameters assigned) to provide subject-specific tumor growth predictions. The poor predictive strength of the current model and the temporally constant proliferation rates hinders the reliability of both untreated (and, most likely, treated) tumor growth predictions. Expanding the model to include an additional term to describe the effect of treatment (i.e., death rate as a function of drug dose or radiation dose) or fitting for a post-treatment proliferation rate would provide a platform to compare observed treatment response to predicted treatment response.

In conclusion, parameters can be accurately estimated and used to predict future tumor growth with low error at the global and local levels, provided that the tumor’s growth is described by the reaction-diffusion equation. However, the in vivo experiments suggest that the reaction-diffusion model consisting of just tumor cell diffusion and logistic growth described by Eq. (1) provides an incomplete description of tumor growth and must be amended to provide better descriptions of in vivo C6 glioma growth in rats.

Supplementary Material

pb046006_suppdata.pdf

Acknowledgments

The authors thank Zou Yue for performing the animal surgeries. We thank the National Institutes of Health for funding through NCI R01CA138599, NCI R21CA169387, NCI U01CA174706, NCI R25CA092043, and the Vanderbilt-Ingram Cancer Center Support Grant (NIH P30CA68485). We thank the Kleberg Foundation for the generous support of our institute’s imaging program. ECR holds a Career Award from the BWF. This work was partially supported by a pilot project from the Physical Sciences in Oncology Center at the H. Lee Moffitt Cancer Center and Research Institute (PIs: Robert Gatenby, M.D., Robert Gillies, Ph.D.).

References

1. Yankeelov TE, Atuegwu N, Hormuth DA, Weis JA, Barnes SL, Miga MI, Rericha EC, Quaranta V. Clinically Relevant Modeling of Tumor Growth and Treatment Response. Sci Transl Med. 2013;5(187):187ps9–187ps9. [PMC free article] [PubMed]
2. Gammon K. Mathematical modelling: Forecasting cancer. Nature. 2012;491(7425):S66–7. [PubMed]
3. Savage N. Modelling: Computing cancer. Nature. 2012;491(7425):S62–3. [PubMed]
4. Kim B, Chenevert TL, Ross BD. Growth kinetics and treatment response of the intracerebral rat 9L brain tumor model: a quantitative in vivo study using magnetic resonance imaging. Clin Cancer Res. 1995;1(6):643–50. [PubMed]
5. Atuegwu NC, Arlinghaus LR, Li X, Welch EB, Chakravarthy BA, Gore JC, Yankeelov TE. Integration of diffusion-weighted MRI data and a simple mathematical model to predict breast tumor cellularity during neoadjuvant chemotherapy. Magn Reson Med. 2011;66(6):1689–96. [PMC free article] [PubMed]
6. Paulus W, Peiffer J. Intratumoral Histologic Heterogeneity of Gliomas - A Quantitative Study. Cancer. 1989;64(2):442–7. [PubMed]
7. Wagner M, Nafe R, Jurcoane A, Pilatus U, Franz K, Rieger J, Steinbach JP, Hattingen E. Heterogeneity in malignant gliomas: a magnetic resonance analysis of spatial distribution of metabolite changes and regional blood volume. J Neurooncol. 2011;103(3):663–72. [PubMed]
8. Russell SM, Elliott R, Forshaw D, Golfinos JG, Nelson PK, Kelly PJ. Glioma vascularity correlates with reduced patient survival and increased malignancy. Surg Neurol. 2009;72(3):242–6. [PubMed]
9. Brown JM, Wilson WR. Exploiting tumour hypoxia in cancer treatment. Nat Rev Cancer. 2004;4(6):437–47. [PubMed]
10. Suit H, Skates S, Taghian A, Okunieff P, Efird JT. Clinical implications of heterogeneity of tumor response to radiation therapy. Radiother Oncol. 1992;25(4):251–60. [PubMed]
11. Au JL, Jang SH, Zheng J, Chen C-T, Song S, Hu L, Wientjes MG. Determinants of drug delivery and transport to solid tumors. J Control Release. 2001;74(1–3):31–46. [PubMed]
12. Gerlee P, Anderson ARA. An evolutionary hybrid cellular automaton model of solid tumour growth. J Theor Biol. 2007;246(4):583–603. [PMC free article] [PubMed]
13. Cai Y, Xu S, Wu J, Long Q. J Theor Biol. 1. Vol. 279. Elsevier; 2011. Coupled modelling of tumour angiogenesis, tumour growth and blood perfusion; pp. 90–101. [PubMed]
14. Anderson ARA, Chaplain MAJ. Continuous and discrete mathematical models of tumor-induced angiogenesis. Bull Math Biol. 1998;60(5):857–99. [PubMed]
15. McDougall SR, Anderson ARA, Chaplain MAJ. Mathematical modelling of dynamic adaptive tumour-induced angiogenesis: Clinical implications and therapeutic targeting strategies. J Theor Biol. 2006;241(3):564–89. [PubMed]
16. Koh DM, Collins DJ. Am J Roentgenol. 6. Vol. 188. American Roentgen Ray Society; 2007. Diffusion-Weighted MRI in the Body: Applications and Challenges in Oncology; pp. 1622–35. [PubMed]
17. Yankeelov TE, Gore JC. Dynamic Contrast Enhanced Magnetic Resonance Imaging in Oncology: Theory, Data Acquisition, Analysis, and Examples. Curr Med Imaging Rev. 2009;3(2):91–107. [PMC free article] [PubMed]
18. Østergaard L. Principles of cerebral perfusion imaging by bolus tracking. J Magn Reson Imaging. 2005;22(6):710–7. [PubMed]
19. Padhani A, Krohn K, Lewis J, Alber M. Imaging oxygenation of human tumours. Eur Radiol. 2007;17(4):861–72. [PMC free article] [PubMed]
20. Christen T, Lemasson B, Pannetier N, Farion R, Segebarth C, Rémy C, Barbier EL. Evaluation of a quantitative blood oxygenation level-dependent (qBOLD) approach to map local blood oxygen saturation. NMR Biomed. 2011;24(4):393–403. [PubMed]
21. Kubota K. From tumor biology to clinical Pet: a review of positron emission tomography (PET) in oncology. Ann Nucl Med. 2001;15(6):471–86. [PubMed]
22. Atuegwu NC, Gore JC, Yankeelov TE. The integration of quantitative multi-modality imaging data into mathematical models of tumors. Phys Med Biol. 2010;55(9):2429–49. [PMC free article] [PubMed]
23. Weis JA, Miga MI, Arlinghaus LR, Li X, Chakravarthy AB, Abramson V, Farley J, Yankeelov TE. A mechanically coupled reaction-diffusion model for predicting the response of breast tumors to neoadjuvant chemotherapy. Phys Med Biol. 2013;58(17):5851–66. [PMC free article] [PubMed]
24. Liu Y, Sadowski SM, Weisbrod AB, Kebebew E, Summers RM, Yao J. Patient specific tumor growth prediction using multimodal images. Med Image Anal. 2014;18(3):555–66. [PMC free article] [PubMed]
25. Neal ML, Trister AD, Cloke T, Sodt R, Ahn S, Baldock AL, Bridge CA, Lai A, Cloughesy TF, Mrugala MM, et al. PLoS One. 1. Vol. 8. Public Library of Science; 2013. Discriminating Survival Outcomes in Patients with Glioblastoma Using a Simulation-Based, Patient-Specific Response Metric; p. e51951. [PMC free article] [PubMed]
26. Hogea C, Davatzikos C, Biros G. An image-driven parameter estimation problem for a reaction-diffusion glioma growth model with mass effects. J Math Biol. 2008;56(6):793–825. [PMC free article] [PubMed]
27. Ellingson B, LaViolette P, Rand SD, Malkin MG, Connelly JM, Mueller WM, Prost RW, Schmainda KM. Spatially quantifying microscopic tumor invasion and proliferation using a voxelwise solution to a glioma growth model and serial diffusion MRI. Magn Reson Med. 2011;65(4):1131–43. [PMC free article] [PubMed]
28. Jbabdi S, Mandonnet E, Duffau H, Capelle L, Swanson KR, Pélégrini-Issac M, Guillevin R, Benali H. Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging. Magn Reson Med. 2005;54(3):616–24. [PubMed]
29. Swanson KR, Rostomily RC, Alvord EC. Br J Cancer. 1. Vol. 98. Cancer Research UK; 2008. A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle; pp. 113–9. [PMC free article] [PubMed]
30. Garg I, Miga MI. Preliminary investigation of the inhibitory effects of mechanical stress in tumor growth. Proc SPIE. 2008:69182L–69182L–11.
31. Colvin DC, Loveless ME, Does MD, Yue Z, Yankeelov TE, Gore JC. Earlier detection of tumor treatment response using magnetic resonance diffusion imaging with oscillating gradients. Magn Reson Imaging. 2011;29(3):315–23. [PMC free article] [PubMed]
32. Whisenant JG, Ayers GD, Loveless ME, Barnes SL, Colvin DC, Yankeelov TE. Assessing reproducibility of diffusion-weighted magnetic resonance imaging studies in a murine model of HER2+ breast cancer. Magn Reson Imaging. 2014;32(3):245–9. [PMC free article] [PubMed]
33. Anderson AW, Xie J, Pizzonia J, Bronen RA, Spencer DD, Gore JC. Effects of cell volume fraction changes on apparent diffusion in human cells. Magn Reson Imaging. 2000;18(6):689–95. [PubMed]
34. Atuegwu NC, Colvin DC, Loveless ME, Xu L, Gore JC, Yankeelov TE. Incorporation of diffusion-weighted magnetic resonance imaging data into a simple mathematical model of tumor growth. Phys Med Biol. 2012;57(1):225–40. [PMC free article] [PubMed]
35. Martin I, Dozin B, Quarto R, Cancedda R, Beltrame F. Computer-based technique for cell aggregation analysis and cell aggregation in in vitro chondrogenesis. Cytometry. 1997;28(2):141–6. [PubMed]
36. Rouzaire-Dubois B, Milandri JB, Bostel S, Dubois JM. Control of cell proliferation by cell volume alterations in rat C6 glioma cells. Pflugers Arch. 2000;440(6):881–8. [PubMed]
37. Joachimowicz N, Pichot C, Hugonin J-P. Inverse scattering: an iterative numerical method for electromagnetic imaging. IEEE Trans Antennas Propag. 1991;39(12):1742–53.
38. Ou J, Ong RE, Yankeelov TE, Miga MI. Evaluation of 3D modality-independent elastography for breast imaging: a simulation study. Phys Med Biol. 2008;53(1):147. [PubMed]
39. Lin LI. A concordance correlation coefficient to evaluate reproducibility. Biometrics. 1989;45(1):255–68. [PubMed]
40. Dice LR. Ecology. 3. Vol. 26. Ecological Society of America; 1945. Measures of the Amount of Ecologic Association Between Species; pp. 297–302.
41. Yankeelov TE, Abramson RG, Quarles CC. Quantitative multimodality imaging in cancer research and therapy. Nat Rev Clin Oncol. 2014;11(11):670–80. [PMC free article] [PubMed]
42. Freyer JP. Role of Necrosis in Regulating the Growth Saturation of Multicellular Spheroids. Cancer Res. 1988;48(9):2432–9. [PubMed]
43. Lacroix M, Abi-Said D, Fourney DR, Gokaslan ZL, Shi W, DeMonte F, Lang FF, McCutcheon IE, Hassenbusch SJ, Holland E, et al. A multivariate analysis of 416 patients with glioblastoma multiforme: prognosis, extent of resection, and survival. J Neurosurg. 2001;95(2):190–8. [PubMed]
44. Iliadis G, Kotoula V, Chatzisotiriou A, Televantou D, Eleftheraki A, Lambaki S, Misailidou D, Selviaridis P, Fountzilas G. Volumetric and MGMT parameters in glioblastoma patients: Survival analysis. BMC Cancer. 2012;12(1):3. [PMC free article] [PubMed]
45. Gevertz JL, Gillies GT, Torquato S. Simulating tumor growth in confined heterogeneous environments. Phys Biol. 2008;5(3):036010. [PubMed]
46. Gyllenberg M, Webb G. Quiescence as an explanation of Gompertzian tumor growth. Growth Dev Aging. 1989;53(1–2):25–33. [PubMed]
47. Sherar MD, Noss MB, Foster FS. Ultrasound backscatter microscopy images the internal structure of living tumour spheroids. Nature. 1987;330(6147):493–5. [PubMed]
48. Folkman J, Hochberg M. Self-regulation of growth in three dimensions. J Exp Med. 1973;138(4):745–53. [PMC free article] [PubMed]
49. Helmlinger G, Netti PA, Lichtenbeld HC, Melder RJ, Jain RK. Solid stress inhibits the growth of multicellular tumor spheroids. Nat Biotech. 1997;15(8):778–83. [PubMed]
50. Giese A, Bjerkvig R, Berens ME, Westphal M. Cost of Migration: Invasion of Malignant Gliomas and Implications for Treatment. J Clin Oncol. 2003;21(8):1624–36. [PubMed]
51. Demuth T, Berens M. J Neurooncol. 2. Vol. 70. Kluwer Academic Publishers; 2004. Molecular Mechanisms of Glioma Cell Migration and Invasion; pp. 217–28. [PubMed]
52. Wong KL, Summers R, Kebebew E, Yao J. Tumor Growth Prediction with Hyperelastic Biomechanical Model, Physiological Data Fusion, and Nonlinear Optimization. In: Golland P, Hata N, Barillot C, Hornegger J, Howe R, editors. Medical Image Computing and Computer-Assisted Intervention – MICCAI 2014 SE - 4. 2014. pp. 25–32. [PMC free article] [PubMed]