Mechanically Coupled Reaction-Diffusion Model to Predict Glioma Growth: Methodological Details

Abstract

Biophysical models designed to predict the growth and response of tumors to treatment have the potential to become a valuable tool for clinicians in care of cancer patients. Specifically, individualized tumor forecasts could be used to predict response or resistance early in the course of treatment, thereby providing an opportunity for treatment selection or adaption. This chapter discusses an experimental and modeling framework in which noninvasive imaging data is used to initialize and parameterize a subject-specific model of tumor growth. This modeling approach is applied to an analysis of murine models of glioma growth.

Keywords: Biophysical stress; Cancer; Diffusion; Finite difference method; Invasion; MRI.

Fig.1

Experimental timeline and estimation of in vivo cell number from DW-MRI data. (a) On day 0, rats are injected intracranially with 10⁵ C6 glioma cells. (b) Jugular catheters are then inserted on day 8. (c) On days 10 through 20, rats are imaged with MRI with 3D gradient echo, DW-MRI, and CE-MRI. (d) CE-MRI is used to identify tumor tissue by subtracting pre-contrast image from the post-contrast image. (e) ADG(x, y, z, t) is then estimated from DW-MRI data. Finally, N(x, y, z, t) is estimated (f) within the tumor tissue using Eq. 9 and ADC(x, y, z, t)

Fig. 2

Tumor growth modeling and prediction flow chart. DW-MRI and CE-MRI data is first acquired in rats at days t_i to t_f. A subset of the total data (t_i to t_n is used to first estimate model parameters using an iterative optimization algorithm. The optimized model parameters are then used in a forward evaluation of the model system to predict tumor growth at the remaining data points (t_{n + 1} to t_f). The error is then assessed between the model and measured values of N(x, y, z, t)

Fig. 3

Algorithm and example forward evaluation of mechanical and tumor cell model. The mechanical model is first solved to calculate the tissue displacement vector {U} due to N(x, y, z, t), Eq. 21. {U} is then used to calculate strain, stress, and σ(x, y, z, t). The new value of D(x, y, z, t)is calculated using Eq. 2 and σ_vm(x, y, z, t). Finally, D(x, y, z, t) is used in Eq. 6 to calculate the value of N(x, y, z, t + 1)

Fig. 4

Iterative parameter optimization approach. A schematic is shown above for the iterative parameter optimization algorithm using the Levenberg-Marquardt method [32, 33]. The model is first evaluated with an initial guess of model parameters, line 1. The objective function is then evaluated with the current set of model parameters, line 2. The optimal model parameters are then determined in an iterative “while-loop” which ceases when stopping criteria are met. At the beginning of each iteration, the Jacobian is built, line 3, and is used to solve for the new guess of model parameters, line 4. The model is then re-evaluated with the new model parameters, line 5, and the objective function is calculated, line 6. Finally, the error is compared to the previously observed lowest error to determine if the new parameter values are acceptable. The optimization ceases when the stopping criteria are met