* Manuscript Numerical solution of a model for brain cancer progression after therapy

Abstract We present a numerical scheme used to investigate a mathematical model of tumor growth which incorporates multiple disparate timescales. We simulate the model with different initial data. The initial conditions explored herein correspond to a small remnant of tumor tissue left after surgical resection. Our results indicate that tumor regrowth begins at the pre‐surgery tumor‐healthy tissue interface and penetrates back into the original tumor area. This growth is rate‐limited by the reformation of the tumor vascular network.


Introduction
We propose a mathematical model for tumor growth based upon the assumption that growth and motility are mutually exclusive phenomena [1,4,10]. Cells exist in one of two states, and can transit between those states at density dependent per capita rates. We use the Holling type III functional response from the mathematical population theory as the form of these transitions. Growth of tumor cells is limited by a single nutrient, here considered to be oxygen. The model presented here is a simplified variation of the one detailed in [14]. For simplicity, we recreate the heuristic derivation here. More detailed descriptions and justifications of the model and its function forms may be found in [14]. A much more complicated three-dimensional model which is essentially and expanded version of this model is implemented on real brain geometry and validated by a set of clinical data in [5], indicating that the model framework proposed in [14] and used in this paper is robust and plausible.
The two states of this model represent cells locked in a proliferative state and cells locked in a migratory state. Net cellular proliferation is a function of available nutrient, as experimentally quantified in [7,15]. Migratory cells do not grow, but die at a constant per capita rate. Migration is assumed to be a random walk through the tumor region, with solid boundaries which permit no cells to invade the surrounding healthy tissue.
Based on these assumptions, we give a word model which describes the system and then formalize the model into a set of partial differential equations. The initial conditions considered in this work correspond to a novel setting for tumor growth. Rather than investigating early tumor development, we assume that such a tumor has already grown, been detected, and removed. Since surgical resection never guarantees removal of all transformed cells [3,8], we consider initial conditions which dictate different distributions of possible remaining cells and study the nature of tumor regrowth, if in fact it does even regrow. In this setting, we assume spherical symmetry in the region of the excised tumor, and we assume that all vasculature has been removed along with the tumor tissue. Thus, the only source of nutrient to the system is from the exterior boundary.
The organization of this paper is as follows. In Section 2 we provide a mathematical description of the model with new initial conditions corresponding to the distribution of the remaining cancer cells after the therapy. The spatial discretization of the model by pseudospectral method is then described in Section 3. In Section 4 we investigate the equation of the model for the tumor cell density in the proliferation state. In Section 5 we investigate the influence of initial conditions on the stability of computations. The numerical experiments with the model are presented in Section 6. Finally, in Section 7 some concluding remarks are given.

Mathematical Model
Let ρ M (r, t) be the tumor cell density in migratory state, ρ G (r, t) be the tumor cell density in proliferation state, and O(r, t) be the oxygen concentration. Here, r denotes radius and t denotes time. We consider the following model equations, which are introduced in [14] for the cell densities and the oxygen concentration: with the total cell density ρ(r, t) = ρ M (r, t)+ρ G (r, t) and the oxygen dependent growth and consumption functions are Here, A and B are the maximal oxygen dependent growth and death rates, respectively. Constants C 1 and C 2 are the oxygen thresholds and p and q are the cooperativity coefficients for their respective processes. The constant σ is considered to be a survival scaling factor, although its biological interpretation is vague. The constant γ is a conversion factor relating oxygen concentration to the cell density. The Laplacian ∆ r ρ M is given by with ∆ r ρ G defined in a similar way. The system (2.1) has to be supplemented with the initial and boundary conditions. As in [14] we will always assume that the boundary conditions take the form where C 0 > 0 is a given constant. We consider a modification of the initial conditions considered in [14] which correspond to a specific tumor development. These conditions are given by As already mentioned in Section 1 in this paper we also consider a novel setting and assume that the tumor has already grown, has been detected and removed.
This corresponds to different initial conditions which dictate different distributions of the remaining cells, and we assume that where O(r, 0) is defined as in [14] by (2.5). Notice that the initial functions (2.5), (2.6) and (2.7) are even with respect to r and can be symmetrically extended over the interval [−R, R]. The overall new model consisting of (2.1), boundary conditions (2.2), and the initial conditions (2.5), (2.6) and (2.7) will be discretized with respect to the space variable r by pseudospectral method and then integrated in time by the code ode15s from Matlab ODE suite [12].

Pseudo-spectral Radial Discretization
To compute numerical approximations to the unknown functions ρ M , ρ G , O we discretize the system (2.1) with respect to r ∈ [−R, R] by pseudo -spectral method based on the Chebyshev-Gauss-Lobatto points r i = −R cos iπ N , for i = 0, 1, . . . , N . To omit the point of singularity r = 0 we apply only odd numbers N . We have chosen this method because it is fast, efficient, and very accurate. Pseudo-spectral methods are exponentially convergent, i.e., the error behaves like O(α N ) for some 0 < α < 1, compare [6,9], while the finite difference methods converge only with polynomial rate, i.e., the error behaves like O(1/N p ) for some integer p, the order of the method. Hence, we can employ a significantly smaller number of grid points N to get a comparable accuracy to finite difference methods with a much larger number of grid points.
The first order derivatives ∂ρ M (r i , t)/∂r can be replaced by the spectrally accurate approximations is the differentiation matrix of the first order, compare [2] and [11]. The approximations (3.1) and the boundary conditions (2.2) lead to which results in the approximations for ρ M (r 0 , t) and ρ M (r N , t) of the form Then, the approximations (3.1) and (3.2) lead to The approximations (3.

3) and (3.4) result in
where '·' stands for component-wise multiplication of vectors. It can be checked that the denominators in (3.2) satisfy the inequalities Therefore, computing the approximations ρ 0 (t) and ρ N (t) from (3.2) does not cause significant rounding errors. However, since the boundary conditions in (2.2) are different for the functions ρ M and O, for the Laplacian ∆ r O, we apply a different approximation than (3.5). For this approximation, we take the advantage from the fact that the values of O are known at the boundaries and are given exactly by C 0 (the boundary condition for O in (2.2)). To derive the approximation for ∆ r O we apply (3.1) with ρ M replaced by O and additionally we apply the spectrally accurate approximation is the differentiation matrix of the second order, compare with [2,11].
and the vector (O) rr (t) is approximated by Here, From (3.6) and (3.7) we obtain the approximation Application of (3.5) and (3.8) to (2.1) results in the semi-discrete system The resulting semi-discrete system (3.9) of differential equations is stiff and in the next sections, we apply the code ode15s from the Matlab ODE suite [12] to compute approximations to ρ M , ρ G , and O. This code is designed for stiff differential systems. It is much more efficient than the code ode45 based on embedded pair of explicit Runge-Kutta methods of order p = 4 and p = 5 which is designed for non-stiff equations.
In this section, we investigate the second equation of the model (2.1). It is used to compute the tumor cell density ρ G (r, t) in the proliferation state and it can be written in the form: with the variable coefficient

(4.2)
Since the functions O(r, t) and ρ(r, t) = ρ M (r, t) + ρ G (r, t) are symmetric about the r = 0 axis, λ(r, t) has the same property. Our goal is to investigate the sign of the function λ(r, t) for t ≥ 0 and r ∈ [0, R]. Note that λ(r, t) evolves differently for different initial functions ρ(r, 0) and O(r, 0). Fig. 1 shows different coefficients λ(r, t) computed by the definition (4.2) with the functions ρ(r, t) and O(r, t) obtained from model (2.1)-(2.2) supplemented by three different sets of initial conditions: Case (a) corresponds to early tumor development without medical treatment, compare with [14], where the modification of this condition is considered; while cases (b) and (c) correspond to tumor development after surgical resection of already grown tumors which may remain transformed cells, [3], [8].
For Fig. 1, the numerical approximations to ρ M (r, t), ρ G (r, t), and O(r, t) are computed by solving the semi-discrete system (3.9) for r ∈ [−R, R]. We apply the code ode15s (see [12]) with AbsT ol = 10 −8 and RelT ol = 10 −6 to integrate (3.9) in time and then we use the definition (4.2) to compute the approximations to λ(r, t). To compute the numerical approximations to ρ M (r, t), ρ G (r, t), and O(r, t), for r ∈ [−R, R], we extend the corresponding initial conditions over the interval [−R, R] by using the fact that the all three initial functions in the cases (a), (b), and (c) are even as functions of r.
Note that, since A > 0, by the definition (4.2), the initial functions in cases (a),(b) and (c) imply Since λ(r, t) is a continuous function, this inequality is maintained also for t > 0 which are close to zero. This observation is confirmed by the numerical experiments presented in Fig. 1. However, Fig. 1 shows that the values of the function λ(r, t) start to be negative for t > T , with some T sufficiently large. Moreover, λ(r, t) < 0, for t > T and even all r ∈ [−R, R], in the cases (b) and (c). This shows that the numerical computations for the model (2.1)-(2.2) are stable provided that stable numerical methods are applied to the semi-discrete system (3.9) to integrate it in time t.
The goal of this paper is to investigate the cases (b) and (c), which correspond to the cancer development from the cells remained after surgical resections of previously grown tumors. The case (a) was investigated in [14]. In Section 5, we present results of numerical experiments for the cases (b) and (c). The experiments show stable computations and confirm the above analysis of the coefficient λ(r, t).

Analysis of Initial Conditions for the Cancer Model in Spherical Coordinates
The goal of this section is to investigate the numerical solutions to ρ G (r, t) and ρ M (r, t) started at t = 0 from different initial functions, which correspond to different surgical resections. We consider three different sets of initial conditions which correspond to a small remnants of tumor cells after surgery. These conditions are given by: We would like to mention again that the case (a) in Section 4 corresponds to early tumor development. Some modification of this condition was investigated in [14].  Fig. 3 presents numerical approximations to ρ G (r, t) and ρ M (r, t) at t = 5, 10, 15 in the cases (b) and (c). We observe that both densities ρ G (r, t) and ρ M (r, t) are somewhat higher in case (b) than in case (c). This can be explained

Numerical Computations of Tumor Cell Densities and Oxygen Concentrations
In this section, we investigate the model (2.1)-(2.2) supplemented by the initial conditions (2.5)-(2.7). Fig. 4 presents tumor cell densities in migratory and proliferation states after t = 8, 9,10,11,12,14,35, and 40 days. The pictures of Fig. 4 show differences between the cell densities developed in different time ranges. We observe that the cells grow rapidly between t = 5 and t = 14 days, their growth slows down between t = 14 and t = 35 and after t = 35 days the cell densities do not change. Therefore, we conclude that, for (2.5)-(2.7), the most significant changes in the tumor development occur between the 8th and the 14th day.
The development of the oxygen concentration is presented in Fig. 5. After about t = 12 days the oxygen concentration stays constant like in Fig. 4

Conclusions
As shown in Fig. 3 and 4, we observe a quick (t < 14 days) regrowth of tumor density back into the vacant cavity. This regrowth is driven primarily by leftover migratory cells converting back into a proliferative state in a suddenly nutrient-rich environment. This regrowth occurs at the boundary of the surgical resection site; cells more towards the core have little opportunity to establish a large colony (compare initial condition densities 1 and 2 in Section 5 with resulting steady state densities in Figures 4 and 5). The inward extent of growth is limited by consumption of diffusible nutrients (Fig. 5). Taken together, the data generated by this model suggest that the process of tumor recurrence after surgery is rate-limited by the reconstruction of the angiogenic network needed to support large tumor structures. Given the timescale disparities between clinically observed bulk tumor growths (order weeks to months, see [13] for example) versus small-scale tumor growth (order days), one could approximate the rate of vascularization by the rate of bulk tumor expansion.