A Mathematical Model for Tumor Cell Population Dynamics Based on Target Theory and Tumor Lifespan
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.
It has been experimentally verified that the ionization
process initiated by radiating particles leads to lesions in
the cells . The negative effect on DNA
structure, makes lesions the most harmful consequence of
radiotherapy , . 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) , ,  and Normal Tissue Complication Probability (NTCP) , .
There have been numerous advances in stochastic modeling of tumor response to radiation treatment, including the linear quadratic model , , cell population dynamics models , , , mixed-effects behavioral models  and cell cycle models . 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  for further information on different types of heterogeneity.
The key concept in understanding XRT biology is the target theory . 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 . 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 , 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 , , . 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 , .
The clinical significance of the intra-tumor heterogeneity of cell phenotypes and cell damage is discussed in , and . 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 , 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 single-strand break () and cells with double-strand 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 . is written as:
Each target may be revived with probability . As described previously , represents the repair probability matrix in the transition from to , as given by:
where and for .
indicates the cell subpopulation with deactivated target(s), where . For , each cell can move from to with the constant time-independent transition rate of .
A cell will undergo radiation-induced 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 :
and for :
These equations produce the following system:
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:
As seen in Eq. (4.1), the newborn cells’ population size is proportional to . Therefore, the birth rate can be taken as:
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.
Suppose that . Therefore:
Eq. (2.2) shows that for . Therefore:
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 .
The following corollary is a direct consequence of lemma (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 .
According to lemma 4.1, we can separate the tumor cells into different sub-populations 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:
Same analysis shows that for :
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 . Therefore, based on tumor population dynamics, the tumor lifespan is defined as:
6 Stability Analysis
Suppose that is an arbitrary integer. System (4.9) can be written as:
where matrix is described as:
where and . Therefore
Let be a complex matrix. For let denote the sum of the absolute values of the non-diagonal 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 ).
Suppose that . If defines as
According to (6.2), for
First consider that and . Therefore:
The main result of this section is written as follows:
For any , and , the system is stable at equilibrium point , where .
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),
Note that the function is an increasing function for and is an integer. Therefore, for we have:
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.
Suppose that is an integer, and the set denotes the value such that the system (6.1) is stable corresponding to all and . Then:
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 .
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.
-  Curtis S B. Lethal and potentially lethal lesions induced by radiation–a unified repair model. 1986. Radiat. Res. 106, 252–279.
-  Wyman C, Kanaar R. DNA double-strand break repair: all’s well that ends well. Annu. Rev. Genet. 40, 363–383. 2006.
-  Hoeijmakers J H. Genome maintenance mechanisms for preventing cancer. Nature 411, 366–374. 2001.
-  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.
-  Dawson A, Hillen T. Derivation of the tumour control probability(TCP) from a cell cycle model. Comput. Math. Methods Med.7, 121–141. 2006.
-  Gay H A, Niemierko A. A free program for calculating EUD-based NTCP and TCP in external beam radiotherapy. Phys.Med.23, 115–125. 2007.
-  Lyman J T. Complication probability as assessed from dose volume histograms. Radiat. Res. 104, 513–519. 1985.
-  Kallman P, Agren A, Brahme A. Tumour and normal tissue responses to fractionated non-uniform dose delivery. Int.J.Radiat.Biol. 62, 249–262. 1992.
-  Fowler J, The linear-quadratic formula and progress in fractionated radiotherapy. Br.J.Radiol.62,679–694. 1989.
-  Quinn T, Sinkala Z. Dynamics of prostate cancer stem cells with diffusion and organism response. BioSystems, 96(1), 69–79. 2009.
-  Sachs R K, Hlatky L R, Hahnfeldt P. Simple ODE models of tumor growth and anti-angiogenic or radiation treatment. Math. Comput. Modell. 33, 1297–1305. 2001.
-  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), 172-177. 2009.
-  Bastogne T,Samson A,Keinj R,Vallois P,Wantz-MÃ©ziÃ¨res S,Pinel S,Bechet D,Barberi-Heyob M. Phenomenological modeling of tumor diameter growth based on a mixed effects model. J. Theor. Biol. 262, 544–55. 2010.
-  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.
-  Michelson S, Leith J T. Tumor Heterogeneity and Growth Contr. A Survey of Models for Tumor-Immune System Dynamics. Birkuser, pp. 295–333. 1997.
-  Rédei G P. Encyclopedia of genetics, genomics, proteomics, and informatics. Springer Science Business Media, Volume 2. 2008.
-  Turner M E Some classes of hit-theory models. Mathematical Biosciences, 23(3), 219-235. 1975.
-  Nomiya T. Discussions on target theory: past and present. Journal of radiation research, 54(6), 1161-1163. 2013.
-  Keinj R, Bastogne T, Vallois P. Multinomial model-based formulations of TCP and NTCP for radiotherapy treatment planning. J.Theor.Biol. 279, 55–62. 2011.
-  Keinj R, Bastogne T, Vallois P. Tumor growth modeling based on cell and tumor lifespans. J.Theor.Biol. 312, 76–86. 2012.
-  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.
-  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.
-  Durrett, R., Foo, J., Leder, K., Mayberry, J., & Michor, F. Intratumor heterogeneity in evolutionary models of tumor progression. Genetics, 188(2), 461–477. 2011.
-  Oroji, A., Omar, M., & Yarahmadian, S. An Ito stochastic differential equations model for the dynamics of the MCF-7 breast cancer cell line treated by radiotherapy. Journal of theoretical biology, 407, 128–137. 2016.
-  Varga, R. S. Geršgorin and his circles (Vol. 36). Springer Science & Business Media. 2010.