Abstract
In this paper, a new mathematical model is developed to represent the interaction between healthy and cancer cells in the human body, focusing on the role of environmental factors and quality of life in the progression of cancer. We have investigated the dynamic effects of inputs on cancer growth, and provide an explanation of how cancer has variable behavior patterns throughout the lives of different patients. The behavior of the system with input and its trajectory patterns are investigated using trajectory patterns and stability analysis. The analysis suggests that a proper treatment method should change the dynamics of the cancer instead of only reducing the population of cancer cells and treatment burden.
Background
Some of the existing studies on cancer therapy are based on the assumption that cancer growth is a time invariant dynamic system [1]. They have been focused on the following objectives:
(i) To reduce and control the tumor mass such that a specified volume is obtained at the end of the treatment [2].
(ii) To lessen the treatment burden of patients. This method considers some constraints on the treatment policy [3].
(iii) To evaluate the number of injected cells that affect the equilibrium points of the immune system and thus may ultimately be dangerous [4].
Previous studies have investigated the effects of therapeutic inputs, which are considered to have direct effects on the system states [5]. However, the behavior of cancer changes as the disease progresses [6]. External stresses that represent destructive inputs, such as environmental and quality of life factors, can cause disability in the DNA repair genes [7]. They can also interfere with, and alter, the functions of regulatory growth signals (TGFα), growthinhibiting signals (TGFβ), and apoptosis (TP53) [6].
We are interested in analyzing how inputs alter the dynamics of the human body, and whether this is the main factor involved in the occurrence of cancer. This paper deals with a comparison between the responses of human body cells in a set of twin identical brothers who live under two different conditions, or inputs. In the next section, a modified model for the human body cells is presented. In the third section, the equilibrium points of the system and its linear approximation are calculated. The fourth section investigates the stability analysis of the system, using stability theorems and trajectory patterns. The fifth and final section discusses the results and summarizes the conclusions.
Model
We have considered a set of identical twin brothers who were born with the same genetic structures. It is assumed that, after the birth, they reside and grow up in two different locations, and under different environmental conditions. One of them, say A, lives in a polluted environment and under external stresses, but brother B does not.
The presented population model originates from [8]. The modifications in this model are described here. The dynamic behavior of an organ of the body that is affected by the cancer is given by the following timevariant equations.
Where:
The auxiliary equation of the system is
In Eq. 1, x and y are the healthy and cancerous cell concentrations, respectively. Because the state variables are physiologically possible, their values are nonnegative, i.e., x ≥ 0 and y ≥ 0. The coefficients a_{1 }and a_{2 }represent the growth rates of the healthy and cancer cells, respectively. The growth of both healthy tissue and tumor decelerates as the concentrations of both the healthy and tumor tissues approach the carrying capacities K_{1 }and K_{2}, respectively [9]. The effect of the immune system is to kill the mutated and cancer cells in proportional rates d_{1 }and d_{2}. The immune system agents force the cancer cells to suicide through apoptosis [10]. The coefficient c represents the proportion of healthy cells whose genome has been disordered by external stresses. These cells initiate the neoplastic transformation, and are added to the tumor cells [11].
The tumor competes with healthy tissue for resources, such as blood, nutrients, and space, so the organ "feels" the tumor [12]. Moreover, these same cancer cells compete with each other. The competition coefficients between different cells are b_{1}, b_{2}, and g.
The effects of the input on the dynamics of the system are introduced by the registration of input history in the system coefficients. The destructive inputs, in addition to increasing cell mutation, predispose the cancer cells to intensification of the growth rate, and reduce the death rate by conferring the ability to evade the intracellular controlling mechanism and the immune system.
Cancer progression represents a macroevolutionary process where karyotype change or genome replacement plays the key dominant role [13]. In the present work, the transformation rate of healthy cells to cancer cells is assumed to be proportional to the input magnitude u. The biotransformation coefficient saturates at a definite limit K_{3}, which is related to the biological limits of body organs and the accumulation of external effects. The other variant parameters have similar formulations, but with different values and rates. When the destructive inputs terminate, their effects remain in the body and may or may not be compensated by the therapeutic inputs or by the body recovery.
The parameter values are listed in Table 1. These values should not be used for clinical applications, and are merely a means to approximate the cancer dynamics in the numerical analysis [8]. It may be shown by implementation of the Lipschitz theorem that Eq. 1 has a unique solution. The stability of the system is analyzed at the specified time intervals when the parameter variation can be ignored owing to low rates.
Table 1. Parameter estimations
Some important properties of the model are presented here, as described in the following two statements.
Statement 1. Let be the solution of Eq. 1. Then the nonnegative orthant is invariant, i.e., .
Proof:
i) Notice that input u has no direct effect on the variables x and y, it can be concluded that if , then for all t ≥ 0.
ii) If the solution of Eq.1 approaches the horizontal axis from nonnegative orthant then, , and y does not decrease, so its value becomes negative.
iii) If the solution of Eq.1 approaches the vertical axis from nonnegative orthant then,, and x does not decrease, so its value becomes negative.
Hence, there is no solution that exits the first orthant. □
Statement 2. The sum of state variables in Eq. 1 is exponent convergent in the first orthant. The region of convergence is contained in
Proof:
Let W = x + y. Then
Also the gradient of W_{1 }may be shown as;
W_{1 }decreases when x and y converge to infinity, thus it has a maximum value.
If the chosen value r is large enough, the above inequality is valid over the whole domain of analysis. The first element of is negative, thus (x_{ext}, y_{ext}) is the single maximum of W_{1 }and the maximum value is W_{1max }= W_{1}(x_{ext}, y_{ext}). After some simplifications and letting
the following inequality is achieved:
It may be shown [14] that for Eq.7 we have;
Finally, with respect to Eqs. 6 and 8, the exponent convergence region is . □
Equilibrium points
Four equilibrium points of Eq. 1 with u = 0 are calculated as:
Where:
The linear approximation of Eq.1 is given by
Where:
The higher order terms are neglected around the origin (x*, y*) = (0, 0) and the last term C is equal to zero at the equilibrium points.
Stability analysis
We study the stability of the equilibrium points of Eq. 1 in this section. The results of the analysis are stated as follows;
i. If the equilibrium point 3 is located in the first orthant and the equilibrium point 2 is not, then the state variables of Eq. 1 will converge to equilibrium point 3 (the healthy state).
ii. If the equilibrium points 2 and 3 are located in the first orthant, then the state variables of Eq. 1 will converge to one of these two points.
iii. If the equilibrium point 2 is located in the first orthant and the equilibrium point 3 is not, then the state variables of Eq. 1 will converge to the equilibrium point 2 (the cancer state).
Proof:
1) At the equilibrium point 1, x = 0, y = 0, the eigenvalues of A are;
As shown in Table 1, the values of a_{1 }are larger than the sum of d_{1 }and c. Also a_{2 }is larger than the sum of d_{2 }and g. Thus the eigenvalues of equilibrium point 1 are always positive and the origin is an unstable node.
2) At the equilibrium point 2, , the eigenvalues of A are;
In Eq.10 y is positive, then β_{2}/α_{2 }< 0. Also, we notice from Table 1 that if the value of β_{2}/α_{2 }is larger than β_{1}, then the equilibrium point 2, if it exists, is a stable node.
3) At the equilibrium point 3, x = α_{6 }+ α_{7}, y = α(α_{6 }+ α_{7}) + β_{1}.
Noticing that x and y in Eq. 11 are in the first orthant, then α_{6 }+ α_{7 }> 0 and α_{1}(α_{6 }+ α_{7}) + β_{1 }> 0. The principal minors of A at this equilibrium point are as follows:
This means that the coefficient matrix A is negative definite at this equilibrium point and it is a stable node.
4) At the equilibrium point 4, x = α_{6 } α_{7}, y = α_{1}(α_{6 } α_{6}) + β_{1}.
Noticing that x and y in Eq.12 are in the first orthant, then α_{6 } α_{7 }> 0 and α_{1}(α_{6 } α_{7}) + β_{1 }> 0. The principal minors of A at this equilibrium point are as follows:
Therefore, if this equilibrium point exists, then it is an unstable saddle point. □
Numerical analysis
In this section, the trajectory patterns of the dynamic systems for the twin brothers are evaluated. The life spans of both brothers are divided into three stages and the dynamic behaviors of the systems are analyzed before and after the presence of inputs. The first stage is from the birth to 850 weeks, the second stage is from 850 to 1950 weeks, and the third stage is from 1950 weeks to 2900 weeks.
The following figures illustrate the trajectory patterns of the dynamic systems representing twin brothers A and B at the beginning of the three stages of their life.
Figure 1 shows the trajectory patterns at the beginning of the first stage of life (at birth) for both brothers A and B. It shows that the equilibrium points 2 and 4 are not located in the positive orthant. Therefore, they are not feasible in any organs of the newborn bodies. This figure also indicates that the only stable equilibrium point of the system is point 3 and the equilibrium point 1 is unstable. The system converges to the equilibrium point 3, if the initial values are in the positive orthant for two newborn brothers.
Figure 1. Trajectory pattern at the beginning of the first stage of life (at birth) for both brothers A and B.
Figure 2 shows the trajectory pattern at the beginning of the second stage of life for the twin brothers, where brother A has been under the effect of destructive inputs during the first stage of life, but brother B has not.
Figure 2. Trajectory patterns at the beginning of the second stage of life, (a) for brother A (affected by cancer) and (b) for brother B (healthy brother).
Figure 2a shows that for brother A, the equilibrium points 2 and 4 are brought to the positive orthant by receiving the destructive input for a limited time, and the trajectories converge to the healthy equilibrium point 3 from a large portion of the state space. However, the behavior of the system for brother B, as shown in Figure 2b, is similar to Figure 1, with the positions of the equilibrium points 2 and 4 showing only minor changes.
Figure 3 shows the trajectory pattern at the beginning of the third stage of life for the twin brothers, where brother A has been under the effect of destructive inputs during the second stage of life, but brother B has not. Figure 3a indicates that the attraction area of equilibrium point 2 is larger than 3 for brother A. The position of equilibrium point 3 remains approximately constant, but the equilibrium points 2 and 4 change as A becomes older and encounters more destructive inputs. These changes lead to an increase in the probability of the incidence of cancer. The destructive input affects the system dynamics and the effects are accumulated in the body of A. The input is assumed to end, but the behavior of the system is changed.
Figure 3. Trajectory patterns at the beginning of the third stage of life, (a) for brother A (affected by cancer) and (b) for brother B (healthy brother).
The body organ of B is changed a little, but the system does not pose any unstable behavior in the first orthant, as is seen in Figure 3b. The small variation is due to aging and the probable effect of limited destructive inputs.
Figure 4 shows the trajectory pattern at the end of the third stage of life for the twin brothers, where brother A has been under the effect of destructive inputs during the third stage of life, but brother B has not. The only stable equilibrium is equilibrium 2 which is located in the nonnegative orthant at the end of the third stage of life for brother A, as is shown in Figure 4a. This brother will die of cancer if the dynamics of the system is not modified. All trajectories in the positive orthant converge on to the cancer equilibrium point 2. Therefore, the cancer is not treated by minimizing the number of cancer cells. There is a low rate of dynamics change in the body of B. The behavior of the system does not transform such that equilibrium point 2 is located in the nonnegative orthant and thus affects the health condition of B. However the equilibrium points 2 and 4 near the positive orthant of state space. The health margin of B is not as it was at birth.
Figure 4. Trajectory patterns at the end of the third stage of life, (a) for brother A (affected by cancer) and (b) for brother B (healthy brother).
Finally, a situation is analyzed in which B (the healthy brother) encounters destructive inputs from 2900 to 3200 weeks.
Figure 5 shows the trajectory patterns for brother B, where he has been under the effect of destructive inputs from 2900 to 3200 weeks. It shows that the destructive inputs change the dynamics of the body organ such that the equilibrium points 2 and 4 enter the nonnegative orthant. The behavior of the system in this figure is similar to Figure 2a for brother A, but the dynamics change in Figure 2a occurs over a period of 850 weeks, while Figure 5 shows the change within 300 weeks.
Figure 5. Trajectory patterns at the 3200th week of of life, for brother B (healthy brother) after the presence of destructive inputs for some years.
Discussion and conclusion
The stable dynamics for the newborn brothers at birth means that they do not have cancer. If the equilibrium points of the system are in the nonnegative orthant then the system will settle in the equilibrium points 2 or 3. The history of the body organ's dynamics during the life of brother B indicates his healthy condition over his life (Figures 1, 2b, 3b, and 4b). This brother is not afflicted by cancer since he has not faced the destructive inputs. The variation of the system dynamics where the destructive inputs are sufficiently effective, demonstrate cancer formation and its progress in the body of brother A (Figures 1, 2a, 3a, and 4a). The stable behavior for A at birth is gradually converted to less stable conditions. The attraction domain of equilibrium point 2, indicating the cancer cells, becomes larger than equilibrium 3, indicating healthy cells. Thus, the cancer incidence is increased. In such a system, if the number of cancer cells reaches a specified limit, then the patient will die. In our model, brother A reaches the stage in which no treatment is able to cure the patient completely and he dies.
This research shows that focusing on the "indelible changes of the system dynamics" is the best way to describe the cancer formation process. The other major conclusions are as follows:
a) The changes in the dynamics of the system occur gradually in the body of the patient. Thus, the duration of effective destructive inputs is directly related to the probability of the cancer occurrence. This can be seen in Figures 2a and 3a.
b) As shown in Figures 2b, 3b, and 4b, the behavior of the system for the healthy brother B, who is not subjected to the destructive inputs, also changes as he becomes older.
c) The effects of environmental and life quality factors are the main causes for the onset of cancer, in bodies where there are no dominant hereditary genetic disorders.
d) The sensitivity of the old healthy brother B to the destructive inputs is greater than that of the young affected brother A, and the system dynamics changes more rapidly for him.
Our next objective is to find a proper therapeutic input that can move equilibrium point 2 out of the nonnegative space. This treatment method would guarantee the impossibility of tumor recurrence. The undesirable changes of the system should be modified by using the corrective inputs, and the treatment not restricted to the system outputs. Furthermore, in future work, the parameters, especially the input, should relate to the qualitative conditions of a patient's life and to the clinical data.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
AG and MK participated in numerous discussions leading to the development of the idea, jointly wrote the manuscript and edited the final version of the paper. All authors read and approved the final manuscript.
Acknowledgements
The authors wish to acknowledge Dr. T. Sadeghi in the Medical Science University of Tehran for her useful comments; and Eng. E.A. Shavazi in the Chemical Engineering Faculty of Newcastle University in Australia for his assistance.
References

Bellomo N, Li NK, Maini PK: On the foundations of cancer modeling: selected topics, speculations, and perspectives.
Mathematical Models and Methods in Applied Sciences 2008, 18(4):593646. Publisher Full Text

Ribba B, Colin T, Schnell S: A multiscale mathematical model of cancer, and its use in analyzing irradiation therapies.

Ghaffari A, Naserifar N: Optimal therapeutic protocols in cancer immunotherapy.
Computers in Biology and Medicine 2010, 40(3):261270. PubMed Abstract  Publisher Full Text

Fister KR, Panetta JC: Optimal control applied to cellcyclespecific cancer chemotherapy.
Siamese Journal of Applied Mathematics 2000, 60(3):10591072. Publisher Full Text

De Pillis LG, Fister KR, Gu W, Collins C, Daub M, Gross D, Moore J, Preskill B: Seeking bangbang solutions of mixed immunechemotherapy of tumors.
Electronic Journal of Differential Equations 2007, 171:124.

Kumar V, Cotran RS, Robbins SL: Robbins Basic Pathology. 7th edition. Philadelphia: W.B. Saunders Co; 2003.

Klein CA, Hölzel D: Systemic cancer progression and tumor dormancy.
Cell Cycle 2006, 5(16):17881798. PubMed Abstract  Publisher Full Text

Freedman HI, Pinho STR: Stability criteria for the cure state in a cancer model with radiation.
Nonlinear Analysis: Real World Applications 2009, 10:27092715. Publisher Full Text

Sachs RK, Hlatky LR, Hahnfeldt P: Simple ODE models of tumor growth and antiangiogenic or radiation treatment.
Mathematical and Computer Modeling 2001, 33:12971305. Publisher Full Text

Wei X, Guo C: Global existence for a mathematical model of the immune response to cancer.
Nonlinear Analysis: Real World Applications 2010.
doi:10.1016/j.nonrwa.2010.02.017

Hirata Y, Bruchovsky N, Aihara K: Development of a mathematical model that predicts the outcome of hormone therapy for prostate cancer.
Journal of Theoretical Biology 2010, 264:517527. PubMed Abstract  Publisher Full Text

Kuang Y, Nagy JD, Elser JJ: Biological stoichiometry of tumor dynamics: mathematical models and analysis.
Dynamics of Continuous, Discrete and Impulsive Systems 2004, (B4):221240.

Heng HHQ, Stevens JB, Bremer SW, Ye KJ, Liu G, Ye CJ: The Evolutionary Mechanism of Cancer.
Journal of Cellular Biochemistry 2010, 109(6):10721084. PubMed Abstract  Publisher Full Text

Slotine JE, Li WL: Applied Nonlinear Control. Englewood Cliffs, New Jersey: PrenticeHall Inc; 1991.