|Home | About | Journals | Submit | Contact Us | Français|
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.
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 [1–3]. 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 [8–11]. Unfortunately, modeling these characteristics frequently requires knowledge of parameters that can only be measured by highly invasive methods or within idealized (in vitro) settings [12–15]. 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 .
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 , blood volume [17,18], blood flow [17,18], hypoxia , oxygen saturation , and metabolism . 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,22–29] have incorporated imaging measurements from MRI, PET, and x-ray computed tomography into mathematical models of tumor growth. Preliminary efforts in both breast  and pancreatic  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 , 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  or simulate a virtual control to assess patient response to radiotherapy . Jbadbi et al  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 . 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  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.
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)):
where N (,t) is the number of tumor cells at three-dimensional position and time t, D() is the tumor cell diffusion coefficient at position , k () is the net tumor cell proliferation at position , and θ is the tumor cell carrying capacity. Note that the proliferation term varies temporally as a function of cell density, N (,t), although it is assumed that the proliferation rate, k (), is temporally constant. As described below, quantitative DW-MRI provides estimates of N (,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 () for each voxel within the tumor, and two D() 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.
Figures 1 and and22 show the approach for the in silico experiments. An initial distribution of tumor cells, N (,t0), was seeded within a rat brain domain. A spatial map of k was determined from DW-MRI estimates of N (,t) from a C6 glioma bearing rat using Eq. (2):
where N (,t0) and N (, t1) represent the distribution of cells at time t0 and t1, respectively, while k () 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 () 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 (,t0) and parameters Ptrue was used to grow an in silico tumor for 10 days (1000 time steps/iterations). Tumor cell distributions, N (,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 (,t0) and N (,t4), respectively). P2 was estimated using tumor cell measurements from days 2 and 4 (N (,t2) and N (,t4), respectively), while P3 was estimated using tumor cell measurements from days 0, 2, and 4 (N (,t0), N (,t2), and N (,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 (,t0) (for P1 and P3) or N (,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 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 (,t5–10), N pred,2 (,t5–10), and N pred,3 (,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).
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 ; 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 . To determine N (,t), the ADC values from the DW-MRI data are then transformed to estimate cell number [33,34] using Eq. (3):
where θ represents the tumor cell carrying capacity, ADCw is the ADC of free water at 37° C (2.5 × 10−3 mm2/s) , ADC(,t) is the ADC value at position 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  with an average cell volume of 908 μm3 .
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 ) 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 (,t) for rat 1; NR1, pred,1 (,t5–10), NR1, pred,2 (,t5–10),and NR1, pred,3 (,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()); i.e., k() kROI in Eq. (1).
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 (,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:
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 (,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 (,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 ).
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 (,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) . 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() 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); ). 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.
Illustrative results of the in silico experiments are shown Figures 3 – 5 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 kP1 – kPtrue (−0.35 to 0.91) compared to kP2 – kPtrue (−0.25 to 0.46) and kP3 – kPtrue (−0.26 to 0.32).
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).
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%).
The results of the in vivo experiments are shown in Figures 6 – 8 and Tables 2–3. 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 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 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 ()) 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 (), where is the average voxel-wise k estimated within the tumor) and the estimated spatially invariant proliferation rate (kROI). ranged from 0.51 to 4.06 day−1, while kROI ranged from 0.94 to 9.94 day−1. kROI was larger than 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.
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 ). 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 (), 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 . In particular, necrotic tissues (which can be relatively large compared to total tumor volume) can strongly influence tumor growth  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 ) 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.
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.).