A Mathematical Model for Tumor Cell Population Dynamics Based on Target Theory and Tumor Lifespan

It has been experimentally verified that the ionization process initiated by radiating particles leads to lesions in the cells [1]. The negative effect on DNA structure, makes lesions the most harmful consequence of radiotherapy [2], [3]. Substantial progress has been made both in the classification and evaluation of XRT treatment planning through probabilistic modeling. The most well-known models are Tumor Control Probability (TCP) [4], [5], [6] and Normal Tissue Complication Probability (NTCP) [7], [8].

There have been numerous advances in stochastic modeling of tumor response to radiation treatment, including the linear quadratic model [4], [9], cell population dynamics models [10], [11], [12], mixed-effects behavioral models [13] and cell cycle models [14]. However, many of these models are constructed to evaluate certain important features of XRT, but they do not incorporate biological tumor damage heterogeneity, which is the focus of our study. We refer the reader to Michelson and Leith [15] for further information on different types of heterogeneity.

The key concept in understanding XRT biology is the target theory [16]. A target is a radio-sensitive site within a cell. Each cell contains a certain number of targets, which may be deactivated after being hit by radiation particles. Moreover, between two consecutive dose fractions, each target may become active again following immune system reaction [17]. Despite the development of several complex interpretations of the target theory, the essential principle entails radiation-induced apoptosis of the organism on account of target(s) inactivation within the organism. Although targets are considered functioning biological units [18], the number of targets and their locations in an organism are not always clear. With regard to cell sensitivity, the majority of models usually assume that cell sensitivity is constant during radiation [19], [20], [21]. The same assumption is also made for cell populations, in that the viability of a surviving cell is similar to an irradiated cell, i.e., all cells are assumed to have the same survival probabilities. However, theses assumptions may not be completely accurate, as there is strong evidence that damaged cells are unable to resist radiation [19], [20].

The clinical significance of the intra-tumor heterogeneity of cell phenotypes and cell damage is discussed in [22], and [23]. As such, providing a definition of a suitable treatment duration is rather a clinical challenge, especially when considering therapeutic response variability. In this regard, Keinj et al. developed a discrete-time Markov chain multinomial model for tumor response [19], which employs the target theory. This model inspects the number of surviving cells in the tumor but does not consider the tumor lifespan to be able to measure the tumor’s response to treatment.

In this study, we model tumor population dynamics via a system of ordinary differential equations. Thereafter, we evaluate the transition rates using a Markov chain. The model is then applied to the special case of $m=3$, which is related to the effect of radiation on cells by dividing them into three subpopulations: cells with no effect ($x_0$), cells with single-strand break ($x_1$) and cells with double-strand breaks ($x_2$). In addition, we analyze the system’s stability in this case as well as the system bifurcation with two parameters.

The paper is organized as follows: Section (2) introduces the general theory and preliminary findings. In section (3), the tumor growth model is discussed terms of a system of ordinary differential equations with Markov chain coefficients. The model calibration is presented in section (4). Thereafter, a new definition for tumor lifespan is proposed in section (5). Three targets in each cell and the model parameters are employed to analytically investigate the obtained ODE system stability in section (6). Finally, section (7) concludes the study.

We have considred the following assumptions in our modeling framework:

1. Cells have the same phenotype but they act independently.

2. In the radiotherapy process, the magnitude of each dose fraction ($u_0$) is constant during treatment (i.e. $u_0=2Gy$). The time lag between two consecutive dose fractions is 24 hours.

3. Each cell consists of $m$ targets, which may be deactivated with probability $q$ after each dose fraction. $\text{P}(i,j)$ represents the treatment probability matrix in the transition from $i$ to $j$ inactive targets, i.e., deactivating $j$ targets when $i$ targets has been disabled before [19]. $\text{P}(i,j)$ is written as:

   $$\text{P}(i,j)=\{\frac{(\frac{m-i}{j-i})q^{j-i}(1-q{)}^{m-j}}{0}\frac{i\le j}{j<i}$$

   (2.1)

4. Each target may be revived with probability $r$. As described previously [19], $\text{R}(i,j)$ represents the repair probability matrix in the transition from $i$ to $j$, as given by:

   (2.2)

   where $\text{R}(m,m)=1$ and $\text{R}(m,j)=0$ for $m\ne j$.

5. $x_i$ indicates the cell subpopulation with $i$ deactivated target(s), where $i=0,...,(m-1)$. For $i\ne j$, each cell can move from $x_i$ to $x_j$ with the constant time-independent transition rate of $\alpha (i,j)$.

6. A cell will undergo radiation-induced apoptosis if all targets are deactivated. The cell death rate in subpopulation $x_i$ is considered as constant, $D_i$.

7. Cells can reproduce if all targets become active. For simplicity, we assume that just before the repair mechanism acts, cells in subpopulation $x_0$ can give birth to new cells proportional to subpopulation $x_0$ with a constant rate of $\beta$. As such, each cell in subpopulation $x_0$ can divide into exactly two daughter cells with probability $\mu$ or it can remain unchanged with probability $(1-\mu )$ between two consecutive dose fractions.

As indicated in Fig. (1), tumor dynamics is generally described as the effect of radiotherapy on the different tumor cell subpopulations. The conservation law for subpopulation $x_i$’s, $i=0,\dots ,m-1$ is written as follows. For $i=0$:

and for $i\ne 0$:

The following corollary is a direct consequence of lemma (4.1).

What dose magnitude is required to remove the tumor completely? A small number of cells may still remain after resection, that are not visible and detectable by MRI. Therefore, it is crucial to know how many dose fractions must be applied to eliminate the remaining cancerous cells.
The tumor lifespan is defined as the minimum number of dose fractions required to remove the entire tumor [24]. Therefore, based on tumor population dynamics, the tumor lifespan is defined as:

where

Let $A_{m\times m}=\left(A_{tk}\right)$ be a complex matrix. For $t\in \{1,...,m\}$ let $R_t=\sum _{k\ne t}\left|A_{tk}\right|$ denote the sum of the absolute values of the non-diagonal entries in the $t$-th row and $D(A_{tt},R_t)$ be the closed disc centered at $A_{tt}$ with radius $R_t$, which is known as Gershgorin disc. Eigenvalue of $A$ lies within at least one of the Gershgorin discs $D(A_{tt},R_t)$ (Gershgorin Theorem [25]).

The main result of this section is written as follows:

In this study, the population dynamics of tumor cells in the process of radiotherapy was examined. A system of differential equations with random variable coefficients was introduced to capture the heterogeneity of cell damage and the repair mechanism between two consecutive dose fractions. Subsequently, a new definition for tumor lifespan was introduced based on tumor population size. Based on the tumor lifespan, the effects of the probability that a target will be inactive after a dose fraction (q) and the probability that a target will reactivate after the repair mechanism (r) were investigated numerically. Our results are in good agreement with previously presented results [20].

The first author appreciates Dr. Ivy Chung and Dr. Ung Ngie Min from Faculty of Medicine, university of Malaya and Prof. Fazlul Sarkar from School of Medicine, Wayne State University for their constructing comments with regard to the manuscript. This study was financially supported by FRGS grant number FP015-2015A, University of Malaya.