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.
Some of the existing studies on cancer therapy are based on the assumption that cancer growth is a time invariant dynamic system . 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 .
(ii) To lessen the treatment burden of patients. This method considers some constraints on the treatment policy .
(iii) To evaluate the number of injected cells that affect the equilibrium points of the immune system and thus may ultimately be dangerous .
Previous studies have investigated the effects of therapeutic inputs, which are considered to have direct effects on the system states . However, the behavior of cancer changes as the disease progresses . External stresses that represent destructive inputs, such as environmental and quality of life factors, can cause disability in the DNA repair genes . They can also interfere with, and alter, the functions of regulatory growth signals (TGF-α), growth-inhibiting signals (TGF-β), and apoptosis (TP53) .
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.
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 . 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 time-variant equations.
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 a1 and a2 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 K1 and K2, respectively . The effect of the immune system is to kill the mutated and cancer cells in proportional rates d1 and d2. The immune system agents force the cancer cells to suicide through apoptosis . 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 .
The tumor competes with healthy tissue for resources, such as blood, nutrients, and space, so the organ "feels" the tumor . Moreover, these same cancer cells compete with each other. The competition coefficients between different cells are b1, b2, 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 macro-evolutionary process where karyotype change or genome replacement plays the key dominant role . 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 K3, 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 . 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.
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
Let W = x + y. Then
Also the gradient of W1 may be shown as;
W1 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 (xext, yext) is the single maximum of W1 and the maximum value is W1max = W1(xext, yext). After some simplifications and letting
the following inequality is achieved:
It may be shown  that for Eq.7 we have;
Four equilibrium points of Eq. 1 with u = 0 are calculated as:
The linear approximation of Eq.1 is given by
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.
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).
1) At the equilibrium point 1, x = 0, y = 0, the eigenvalues of A are;
As shown in Table 1, the values of a1 are larger than the sum of d1 and c. Also a2 is larger than the sum of d2 and g. Thus the eigenvalues of equilibrium point 1 are always positive and the origin is an unstable node.
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. □
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 2-a 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 2-b, 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 3-a 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 3-b. 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 4-a. 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 2-a for brother A, but the dynamics change in Figure 2-a 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, 2-b, 3-b, and 4-b). 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, 2-a, 3-a, and 4-a). 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 2-a and 3-a.
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.
The authors declare that they have no competing interests.
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.
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.
Mathematical Models and Methods in Applied Sciences 2008, 18(4):593-646. Publisher Full Text
Siamese Journal of Applied Mathematics 2000, 60(3):1059-1072. Publisher Full Text
Nonlinear Analysis: Real World Applications 2009, 10:2709-2715. Publisher Full Text
Mathematical and Computer Modeling 2001, 33:1297-1305. Publisher Full Text
Nonlinear Analysis: Real World Applications 2010.