A Mathematical Model for Tumor Cell Population Dynamics Based on Target Theory and Tumor Lifespan
Abstract
Radiation Therapy (XRT) is one of the most common cancer treatment methods. In this paper, a new mathematical model is proposed for the population dynamics of heterogeneous tumor cells following external beam radiation treatment. According to the Target Theory, the tumor population is divided into different subpopulations based on the diverse effects of ionizing radiation on human cells. A hybrid model consists of a system of differential equations with random variable coefficients representing the transition rates between subpopulations is proposed. This model is utilized to simulate the dynamics of cell subpopulations within a tumor. The model also describes the cell damage heterogeneity and the repair mechanism between two consecutive dose fractions. As such, a new definition of tumor lifespan based on population size is introduced. Finally, the stability of the system is studied by using the Gershgorin theorem. It is proven that the probability of target inactivity post radiation plays the most important role in the stability of the system.
1 Introduction
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 wellknown 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], mixedeffects 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
radiosensitive 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 radiationinduced 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 intratumor 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 discretetime 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 , which is related to the effect of radiation on cells by dividing them into three subpopulations: cells with no effect (), cells with singlestrand break () and cells with doublestrand breaks (). 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.
2 Modeling assumptions and framework
We have considred the following assumptions in our modeling framework:

Cells have the same phenotype but they act independently.

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

Each cell consists of targets, which may be deactivated with probability after each dose fraction. represents the treatment probability matrix in the transition from to inactive targets, i.e., deactivating targets when targets has been disabled before [19]. is written as:
(2.1) 
Each target may be revived with probability . As described previously [19], represents the repair probability matrix in the transition from to , as given by:
(2.2) where and for .

indicates the cell subpopulation with deactivated target(s), where . For , each cell can move from to with the constant timeindependent transition rate of .

A cell will undergo radiationinduced apoptosis if all targets are deactivated. The cell death rate in subpopulation is considered as constant, .

Cells can reproduce if all targets become active. For simplicity, we assume that just before the repair mechanism acts, cells in subpopulation can give birth to new cells proportional to subpopulation with a constant rate of . As such, each cell in subpopulation can divide into exactly two daughter cells with probability or it can remain unchanged with probability between two consecutive dose fractions.
3 Model derivation
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 ’s, is written as follows. For :
(3.1) 
and for :
(3.2) 
These equations produce the following system:
(3.3)  
4 Model Calibration
The probability that a cell will remain in after radiation is . Therefore, the average number of births in one day after applying the dose fraction and just before the dose fraction is equal to:
(4.1)  
As seen in Eq. (4.1), the newborn cells’ population size is proportional to . Therefore, the birth rate can be taken as:
(4.2) 
Lemma 4.1.
Suppose that a cell has deactivated target(s) just before the application of a dose fraction and . After treatment and repair,

represents the probability that a cell with deactivated target(s) just before the application of a dose fraction has deactivated target(s) right before the application of the next dose fraction.

An average number of cells move from into .

For fixed and , the map is increasing.
Proof.

Now assume that a cell has deactivated targets just before applying a dose fraction. After radiation and right before the repair mechanism, this cell may remain in subpopulation with probability , or it may move to subpopulation , , with probability . Following repair, this cell may move from subpopulation to subpopulation with probability . Therefore, the probability of transitioning from subpopulation into subpopulation after treatment and repair (one day) is .

The effect of treatment and repair on one cell is independent of the rest of the cells. Therefore, the average number of cells moving from subpopulation into subpopulation is equal to .

See [20].
∎
The following corollary is a direct consequence of lemma (4.1).
Corollary 4.1.
With the same assumptions described in Lemma (4.1):

The cells’ transition rate from subpopulation into subpopulation is equal to .

The death rate of subpopulation is .
Remark 4.1.
According to lemma 4.1, we can separate the tumor cells into different subpopulations according to their sensitivity to the radiation. Therefore, the death rate in a subpopulation with more deactivated targets is higher than a subpopulation with fewer deactivated targets, which can be interpreted as the treatment heterogeneity in the model.
Now, starting with subpopulation , cells give birth at a constant rate of . In addition, cells transit from subpopulation into subpopulation at a rate of . Conversely, cells move from subpopulation into subpopulation at a rate of or may die at a rate of , where is the transition matrix. Hence, for we have:
(4.5) 
By substituting Eq. (4.5) in Eq. (3.1) we get:
(4.6) 
Same analysis shows that for :
(4.7) 
5 Tumor lifespan
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:
(5.1) 
where
(5.2) 
6 Stability Analysis
Suppose that is an arbitrary integer. System (4.9) can be written as:
(6.1) 
where matrix is described as:
(6.2) 
where and . Therefore
and for
Let be a complex matrix. For let denote the sum of the absolute values of the nondiagonal entries in the th row and be the closed disc centered at with radius , which is known as Gershgorin disc. Eigenvalue of lies within at least one of the Gershgorin discs (Gershgorin Theorem [25]).
Lemma 6.1.
Suppose that . If defines as
(6.4) 
then

(6.5) 
(6.6)
Proof.

According to (6.2), for
(6.7) 
First consider that and . Therefore:
(6.8) Moreover, for
(6.9)
∎
The main result of this section is written as follows:
Theorem 6.1.
For any , and , the system is stable at equilibrium point , where .
Proof.
It is enough to show that all eigenvalues of matrix have negative real parts . To provide this we will show that for any , and any point of Gershgorin circles have negative real part, where . For this purpose we apply Gershgorin Theorem on matrix . Based on Lemma (6.1),
(6.10) 
Note that the function is an increasing function for and is an integer. Therefore, for we have:
(6.11) 
where .
Consequently,
and where (Figure (2)). This shows that the Gershgorin Circles belong to the left half
of real line. In addition, according to Gershgorin
Theorem, each eigenvalue of matrix
belongs in one of Gershgorin discs.
Therefore, each eigenvalue of
matrix has negative real part.
Hence, every eigenvalue of matrix has negative real part. This completes the proof.
∎
Theorem 6.2.
Suppose that is an integer, and the set denotes the value such that the system (6.1) is stable corresponding to all and . Then:
(6.12) 
7 Conclusion
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].
Acknowledgement
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 FP0152015A, University of Malaya.
References
 [1] Curtis S B. Lethal and potentially lethal lesions induced by radiation–a unified repair model. 1986. Radiat. Res. 106, 252–279.
 [2] Wyman C, Kanaar R. DNA doublestrand break repair: all’s well that ends well. Annu. Rev. Genet. 40, 363–383. 2006.
 [3] Hoeijmakers J H. Genome maintenance mechanisms for preventing cancer. Nature 411, 366–374. 2001.
 [4] Zaider M, Minerbo G N. Tumour control probability:a formulation applicable to any temporal protocol of dose delivery. Phys. Med. Biol. 45, 279–293. 2000.
 [5] Dawson A, Hillen T. Derivation of the tumour control probability(TCP) from a cell cycle model. Comput. Math. Methods Med.7, 121–141. 2006.
 [6] Gay H A, Niemierko A. A free program for calculating EUDbased NTCP and TCP in external beam radiotherapy. Phys.Med.23, 115–125. 2007.
 [7] Lyman J T. Complication probability as assessed from dose volume histograms. Radiat. Res. 104, 513–519. 1985.
 [8] Kallman P, Agren A, Brahme A. Tumour and normal tissue responses to fractionated nonuniform dose delivery. Int.J.Radiat.Biol. 62, 249–262. 1992.
 [9] Fowler J, The linearquadratic formula and progress in fractionated radiotherapy. Br.J.Radiol.62,679–694. 1989.
 [10] Quinn T, Sinkala Z. Dynamics of prostate cancer stem cells with diffusion and organism response. BioSystems, 96(1), 69–79. 2009.
 [11] Sachs R K, Hlatky L R, Hahnfeldt P. Simple ODE models of tumor growth and antiangiogenic or radiation treatment. Math. Comput. Modell. 33, 1297–1305. 2001.
 [12] Gámez, M., López, I., Garay, J., & Varga, Z. Observation and control in a model of a cell population affected by radiation. Biosystems, 96(2), 172177. 2009.
 [13] Bastogne T,Samson A,Keinj R,Vallois P,WantzMÃ©ziÃ¨res S,Pinel S,Bechet D,BarberiHeyob M. Phenomenological modeling of tumor diameter growth based on a mixed effects model. J. Theor. Biol. 262, 544–55. 2010.
 [14] Kirkby N F, Burnet N G, Faraday D B F. Mathematical modelling of the response of tumour cells to radiotherapy. Nucl. Instrum. MethodsPhys. Res. Sect. B188, 210–215. 2002.
 [15] Michelson S, Leith J T. Tumor Heterogeneity and Growth Contr. A Survey of Models for TumorImmune System Dynamics. Birkuser, pp. 295–333. 1997.
 [16] Rédei G P. Encyclopedia of genetics, genomics, proteomics, and informatics. Springer Science Business Media, Volume 2. 2008.
 [17] Turner M E Some classes of hittheory models. Mathematical Biosciences, 23(3), 219235. 1975.
 [18] Nomiya T. Discussions on target theory: past and present. Journal of radiation research, 54(6), 11611163. 2013.
 [19] Keinj R, Bastogne T, Vallois P. Multinomial modelbased formulations of TCP and NTCP for radiotherapy treatment planning. J.Theor.Biol. 279, 55–62. 2011.
 [20] Keinj R, Bastogne T, Vallois P. Tumor growth modeling based on cell and tumor lifespans. J.Theor.Biol. 312, 76–86. 2012.
 [21] O’Rourke S F C, McAneney H, Hillen T. Linear quadratic and tumour control probability modeling in external beam radiotherapy J. Math. Biol. 58, 799–817. 2009.
 [22] Gupta, P. B., Fillmore, C. M., Jiang, G., Shapira, S. D., Tao, K., Kuperwasser, C., & Lander, E. S. Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells. Cell, 146(4), 633–644. 2011.
 [23] Durrett, R., Foo, J., Leder, K., Mayberry, J., & Michor, F. Intratumor heterogeneity in evolutionary models of tumor progression. Genetics, 188(2), 461–477. 2011.
 [24] Oroji, A., Omar, M., & Yarahmadian, S. An Ito stochastic differential equations model for the dynamics of the MCF7 breast cancer cell line treated by radiotherapy. Journal of theoretical biology, 407, 128–137. 2016.
 [25] Varga, R. S. Geršgorin and his circles (Vol. 36). Springer Science & Business Media. 2010.