A Multilayer Grow-or-Go Model for GBM: Effects of Invasive Cells and Anti-Angiogenesis on Growth

Abstract

The recent use of anti-angiogenesis (AA) drugs for the treatment of glioblastoma multiforme (GBM) has uncovered unusual tumor responses. Here, we derive a new mathematical model that takes into account the ability of proliferative cells to become invasive under hypoxic conditions; model simulations generate the multilayer structure of GBM, namely proliferation, brain invasion, and necrosis. The model is able to replicate and justify the clinical observation of rebound growth when AA therapy is discontinued in some patients. The model is interrogated to derive fundamental insights int cancer biology and on the clinical and biological effects of AA drugs. Invasive cells promote tumor growth, which in the long run exceeds the effects of angiogenesis alone. Furthermore, AA drugs increase the fraction of invasive cells in the tumor, which explain progression by fluid-attenuated inversion recovery (FLAIR) signal and the rebound tumor growth when AA is discontinued.

Appendix: Model Parameters and Implementation

The parameters of the simulations of Table 1 are shown in Table 2. For implementation, the same point of view as in Colin et al. (2013) is chosen. The various numerical schemes and the main difficulties are presented.

1.1 Level-Set and Penalty Method

Brain geometry is very complex. Meshing it is difficult and involves a huge number of points. A level-set formulation is used in order to describe its boundary and to impose our boundary conditions with penalty methods (Angot et al. 1999). This allows us to use a Cartesian mesh. The distance to the boundary is chosen as a level-set function:

$$\begin{aligned} l :x \rightarrow {\left\{ \begin{array}{ll} + d \bigl (\mathbf {x}, \partial \Omega _\text {brain}\bigr ) &{} \text {if }\, \mathbf {x} \in \Omega _\text {brain}, \\ - d\bigl ( \mathbf {x}, \partial \Omega _\text {brain}\bigr ) &{} \text {if }\, \mathbf {x} \notin \Omega _\text {brain}, \end{array}\right. } \end{aligned}$$

where $d$ denotes the classical distance in $\mathbb {R}^2$ or $\mathbb {R}^3$.

All the boundary conditions are no-flux conditions on $\partial \Omega _C$. It is presented here on the equation satisfied by the pressure $\pi$ obtained b taking the divergence of Eq. (11) and using Eq. (14). This equation together with the boundary condition $(K \nabla \pi )\cdot \mathbf {n} = 0$ (where K is a tensor) is replaced by:

$$\begin{aligned} - \nabla \cdot \bigl ( K^\varepsilon \nabla \pi ^\varepsilon \bigr ) = \left( m(C)P- \frac{B}{||B ||_{1, \Omega _\text {brain}}} \int _{\Omega _\text {brain}} m(C) P \text {d}\Omega \right) \text { on }\, \Omega , \end{aligned}$$

with $\pi ^\varepsilon = 0$ on $\partial \Omega ,\, \Omega$ being a box in which $\Omega _\text {brain}$ is embedded. Here, $K^\varepsilon$ is defined by:

$$\begin{aligned} K^\varepsilon = \chi _{\mathbf {x} \in \Omega _\text {brain}}K + \varepsilon \chi _{\mathbf {x} \notin \Omega _\text {brain}}. \end{aligned}$$

and $\chi _{\mathbf {x} \in \Omega _\text {brain}}$ denotes the characteristic function of $\Omega _\text {brain}$. As $\varepsilon \rightarrow 0,\, \pi ^\varepsilon$ converge to $\pi$ solution of (10):

$$\begin{aligned} - \nabla \cdot \bigl ( K \nabla \pi \bigr ) = \left( m(C)P- \frac{B}{||B ||_{1, \Omega _\text {brain}}} \int _{\Omega _\text {brain}} m(C) P \text {d}\Omega \right) \end{aligned}$$

on $\omega _\text {brain}$ with $(K \nabla \pi )\cdot \mathbf {n} = 0$ on $\partial \Omega _\text {brain}$.

1.2 Numerical Scheme

Concerning time discretization, a second-order splitting scheme is used. The main advantage of this strategy is that the resolution of system Eqs. (2), (8), (4), (7), and (15) on one time step is decomposed into a sequence of well-known problems: diffusion, convection, and ODEs. In order to benefit from Cartesian meshes, well-known high-order numerical schemes are used to discretize the equations in space (Jiang and Peng 2000; Eymard et al. 2000).

The fluxes are computed by means of centered discretizations. This provides second-order schemes. We refer to Drblikova and Mikula (2007), Drblikova et al. (2009) for more details on anisotropic diffusion. For diffusion equation, we use a Crank–Nicolson discretization in time in order to ensure second-order accuracy. Due to the diffusion term, the invasive cell density is quite smooth. The chemotaxis part on invasive cells is solved with a WENO scheme (Jiang and Peng 2000).

The coupling between advection and tumor growth on each density is more difficult to solve. Indeed, one has to be careful on the advection and mitosis to avoid a loss of mass: As advection at velocity $\mathbf {v}$ is supposed to balance mass variation, they have to be discretized exactly at he same time. Our ”splitting philosophy” does not give an accurate outcome. The mass conservation could be improved using a RK2 scheme on time to solve Eqs. (2)–(11) as described in Colin et al. (2013).

Here, we choose another strategy. Let us rewrite the advection as a transport term and a ”divergence” term: $\nabla \cdot \bigl ( u \mathbf {v} \bigr ) = \mathbf {v} \cdot \bigl (\nabla u \bigr ) + \bigl (\nabla \cdot \mathbf {v}\bigr ) u$. By observing that both the transport part—”$\partial _{t}u + \mathbf {v} \cdot \bigl (\nabla u \bigl ) = \ldots$”—and the other part (source terms and divergence part - $div(\mathbf {v}) u$) preserve the total density ”$P+I+N+B= 1$” (Eq. (10)), we propose the following discretization:

• 1. compute $\widetilde{\mathcal {M}}_n := m(C_n) P_n$,

• 2. compute the loss $\widetilde{\mathcal {B}}_n := - \frac{B_n}{||B_n ||_{1, \Omega _\text {brain}}} \int _{\Omega _\text {brain}} \mathcal {M}_n \text {d}\Omega$

• 3. solve $\nabla \cdot \bigl ( K \nabla \widetilde{\pi }_{n+1} \bigr ) = \bigl ( \widetilde{\mathcal {M}}_n- \widetilde{\mathcal {B}}_n \bigr )$ with a finite-volume scheme,

• 4. compute $\widetilde{\mathbf {v}}_{n+1} = K \nabla \widetilde{\pi }_{n+1}$,

• 5. compute

  $$\begin{aligned} \widetilde{P}_{n+1} =\,&P_n + \Delta t \bigl [ F_\text {Weno}(P_n,\widetilde{\mathbf {v}}_{n+1}) - \bigl ( \widetilde{\mathcal {M}}_n -\widetilde{\mathcal {B}}_n \bigr )P_n +\widetilde{\mathcal {M}}_n \bigr ] ,\\ \widetilde{u}_{n+1} =\,&u_n + \Delta t \bigl [ F_\text {Weno}(u_n,\widetilde{\mathbf {v}}_{n+1}) -\bigl ( \widetilde{\mathcal {M}}_n -\widetilde{\mathcal {B}}_n \bigr ) u_n \bigr ]&u=I,N, \\ \widetilde{B}_{n+1} =\,&B_n + \Delta t \bigl [ F_\text {Weno}(B_n,\widetilde{\mathbf {v}}_{n+1}) - \bigl ( \widetilde{\mathcal {M}}_n -\widetilde{\mathcal {B}}_n \bigr )B_n - \widetilde{\mathcal {B}}_n \bigr ], \end{aligned}$$

• 6. compute $\mathcal {M}_{n+1} := m(C_n) \widetilde{P}_{n+1}$,

• 7. compute $\mathcal {B}_{n+1} := - \frac{\widetilde{B}_{n+1}}{||\widetilde{B}_{n+1} ||_{1, \Omega _\text {brain}}} \int _{\Omega _\text {brain}} \mathcal {M}_{n+1} \text {d}\Omega$,

• 8. solve $\nabla \cdot \bigl ( K \nabla \pi _{n+1} \bigr ) = \bigl ( \mathcal {M}_{n+1}- \mathcal {B}_{n+1} \bigr )$

• 9. compute $\mathbf {v}_{n+1} = K \nabla \pi _{n+1}$

• 10. compute

  $$\begin{aligned} P_{n+1} =&\frac{P_n+\widetilde{P}_{n+1}}{2} + \frac{\Delta t}{2} \biggl [ F_\text {Weno}(\widetilde{P}_{n+1},\mathbf {v}_{n+1}) - \bigl ( \mathcal {M}_{n+1} - \mathcal {B}_{n+1} \bigr )\widetilde{P}_{n+1}\\&+\mathcal {M}_{n+1} \biggr ] ,\\ I_{n+1} =&\frac{I_n+\widetilde{I}_{n+1}}{2} + \frac{\Delta t}{2} \biggl [ F_\text {Weno}(\widetilde{I}_{n+1},\mathbf {v}_{n+1}) - \bigl ( \mathcal {M}_{n+1} - \mathcal {B}_{n+1} \bigr )\widetilde{I}_{n+1}\biggr ] ,\\ N_{n+1} =&\frac{N_n+\widetilde{N}_{n+1}}{2} + \frac{\Delta t}{2} \biggl [ F_\text {Weno}(\widetilde{N}_{n+1},\mathbf {v}_{n+1}) - \bigl ( \mathcal {M}_{n+1} - \mathcal {B}_{n+1} \bigr )\widetilde{N}_{n+1}\biggr ] ,\\ B_{n+1} =&\frac{B_n+\widetilde{B}_{n+1}}{2} \!+\! \frac{\Delta t}{2} \biggl [ F_\text {Weno}(\widetilde{B}_{n+1},\mathbf {v}_{n+1}) \!-\! \bigl ( \mathcal {M}_{n+1} \!-\! \mathcal {B}_{n+1} \bigr )\widetilde{B}_{n+1} \!-\! \mathcal {B}_{n+1} \biggr ], \end{aligned}$$

where $F_\text {Weno}(u, \mathbf {v})$ denotes the numerical computation of $\mathbf {v}\nabla u$ with an WENO5 Scheme.

This discretization provides a second-order accuracy in time (Heun’s scheme as in Colin et al. 2013). There does not remain any splitting between advection and mass variation due to tumor growth and loss of $B$. The other terms—such as transition between proliferative and invasive states, necrosis, oxygen concentration—are still dealt with by a splitting scheme.