Model-Based Evaluation of Spontaneous Tumor Regression in Pilocytic Astrocytoma

Pilocytic astrocytoma (PA) is the most common brain tumor in children. This tumor is usually benign and has a good prognosis. Total resection is the treatment of choice and will cure the majority of patients. However, often only partial resection is possible due to the location of the tumor. In that case, spontaneous regression, regrowth, or progression to a more aggressive form have been observed. The dependency between the residual tumor size and spontaneous regression is not understood yet. Therefore, the prognosis is largely unpredictable and there is controversy regarding the management of patients for whom complete resection cannot be achieved. Strategies span from pure observation (wait and see) to combinations of surgery, adjuvant chemotherapy, and radiotherapy. Here, we introduce a mathematical model to investigate the growth and progression behavior of PA. In particular, we propose a Markov chain model incorporating cell proliferation and death as well as mutations. Our model analysis shows that the tumor behavior after partial resection is essentially determined by a risk coefficient γ, which can be deduced from epidemiological data about PA. Our results quantitatively predict the regression probability of a partially resected benign PA given the residual tumor size and lead to the hypothesis that this dependency is linear, implying that removing any amount of tumor mass will improve prognosis. This finding stands in contrast to diffuse malignant glioma where an extent of resection threshold has been experimentally shown, below which no benefit for survival is expected. These results have important implications for future therapeutic studies in PA that should include residual tumor volume as a prognostic factor.


Supplementary Material
Transition rates Formally, the TGP process (X t ) t≥0 is modeled as a Markov process on the state space S = {0, 1, 2, ..., N, E}. The dynamics is determined by the rate matrix Q = (q(k, l)) k,l∈S with q(k, l) = where u and v represent the mutation probabilities from wild-type to type-I cells and from type-I cells to type-II cells, respectively. These rates induce that the states N and E are absorbing states of the process.

Absorption probabilities
In order to calculate the absorption probabilities of the Markov process determined by the rates in equation (S1), a sub-process is investigated. This sub-process is characterized by the state space S = {1, 2, 3, ...., N, E} and rate matrix Q = ( q(k, l)) k,l∈S which is obtained from the original rates given in (S1) by eliminating state 0 and setting u = 0 such that By dening q(k) := − q(k, k), we get We further dene transition probabilities which look as follows.
The absorption probabilities of the processes determined by the rate matrix q and the transition matrix P are equal. This holds due to the fact that P implies an equivalent process where we only eliminated transitions into the same state which only inuences the time-scale of the process. The absorption probabilities for the process which is determined by the transition matrix P = (p i,j ) i,j∈ S is obtained as follows. Denote by the vector α N = α N (i, v) i∈ S the absorption probabilities, where α N (i, v) equals the absorption probability in state N starting from state i. First step analysis yields It holds that α N (E, v) = 0, α N (N, v) = 1 and therefore By multiplying each line with its denominator, one gets an equivalent system P α N = b for a (N − 1) × (N − 1) matrix P and α N := (α N (i, v)) i=1,...,N −1 . This system reads in tableau form as follows.
We are interested in the absorption probability α N (1, v) = α N (1, v), i.e. the probability of getting absorbed in state N when the process is started with a single type-I cell. We use Cramer's rule which reads where P 1 is the matrix formed by replacing the rst column of P by the column vector b. We calculate det P rst. By induction over N the general structure can be inferred. It holds that Furthermore, det P = 120(v 5 + 25v 4 + 100v 3 + 100v 2 + 25v + 1) = 5!(v 5 + 5 2 v 4 + 10 2 v 3 + 10 2 v 2 + 5 2 v + 1) for N = 6.
Therefore, we conclude that the general form of det P is given by Here, P N (x) denotes the Legendre polynomials [2] which are the particular solutions to the Legendre dierential equation The second determinant is calculated as follows. The matrix P 1 has the structure Therefore, the determinant can be calculated by applying Laplace expansion along the rst column and evaluating the determinant of the remaining triangular matrix.

Asymptotic absorption probabilities
Here, we derive the absorption probability in dependency of the risk coecient γ as the system size N tends to innity. The second assumption on the parameter regime of the TGP model reads One important property of the Legendre polynomials is the integral representation (see [2]) This integral representation is used in order to calculate the limit for the denominator in (S8) as N goes to innity.
The exchange of the limit and the integral is justied by using Lebesgue's dominated convergence theorem since there is an integrable majorant which can be derived as follows.
for N suciently large and since the sequence 1 + 1 N N is monotonously increasing.
The limit in (S10) can be calculated as follows.
lim N →∞ . Hence, we obtain for (S10) where I n denotes the modied Bessel function of the rst kind. This function can be expressed by the series where Γ(x) denotes the Gamma function, see [2] for details. Now, equation (S9) implies A comparison between simulation results produced by sampling trajectories of the process (X t ) t≥0 , absorption probabilities obtained by using the exact formula (S7) and the asymptotic equation (S12) is given in Table S1.

Derivation of the regression function
The regression function β N γ (ρ) can be estimated by a diusion approximation of the TGP process. In order to achieve this, a derivation in [1] is used and extended. There, it is shown that it suces to investigate a modied process (Y t ) t≥0 . This process is determined by the original rates given by (S1) with the following modication. The rate for a type-I mutation equals u = 0 in q(k, l) if k > 0. Hence, type-I mutations are not allowed if type-I cells are already present in the system. The modication can be justied by the assumption u 1 N which allows to treat each mutant lineage independently. Note that this decomposition was already used in the calculation of the absorption probabilities. There, the state space has been reduced by state 0 since the occurrence of the successful mutant is assumed at the beginning. For the process (Y t ) t≥0 , this reduction of the state space is not sensible since we want to investigate the regression probability, i.e. the probability of reaching state 0. Hence, both modied processes only dier in the possibility of reaching state 0. Under appropriate time and space scaling, (Y t ) t≥0 can be asymptotically approximated for N → ∞ by the Wright-Fisher diusion process Z t on [0, 1]. The details of this construction can be found in [1]. The important connection between the processes (Y t ) t≥0 and (Z t ) t≥0 is that Z t = 0 implies X t = 0. Therefore, we approximate β N γ (ρ) ≈ β γ (ρ) where β γ (ρ) is the probability that Z t reaches 0 when starting in ρ, β γ (ρ) := P(Z t = 0 for some t > 0|Z 0 = ρ), 0 ≤ ρ ≤ 1.

It holds that
where the constant c is determined by the condition β γ (0) = 1. See [1, Lemma 6.9] for details and a rigorous derivation of this approximation.
Here, we express the series representation (S13) in terms of Bessel functions and derive the constant c as follows. In the rst step, the indices of the sum are adjusted. In the second step, we used that Γ(n) = (n − 1)! for n ∈ N.