On Study of Immune Response to Tumor Cells in Prey-Predator System
This paper aims to develop the mathematical model that explores the immune response to a tumor system as a prey-predator system. A deterministic model defining the dynamics of tumor growth progression and regression has been analyzed. Our analysis indicates the tumor recurring and dormancy on the cellular level in combination with resting and hunting cells. The model considered in the present study is a generalization of El-Gohary (2008) by introducing the Michaelis-Menten function. This function describes the stimulation process of the resting cells by the tumor cells in the presence of tumor specific antigens. Local and global stability analysis have been performed along with the numerical simulation to support our findings.
Cancer is one of the leading causes of death worldwide. Each year, the American Cancer Society projects the number of new cancer cases and death to estimate the contemporary cancer burden. In 2014, there will be an estimated 16 lakh new cancer cases diagnosed and around 5 lakh cancer deaths in the U.S. . Not much is known about its mechanisms of establishment and destruction. It is the second fatal disease after the cardiovascular diseases . According to World Economic Forum (WEF), cancer is among one of the three greatest risks to the global economy due to escalating cost of care, the threat to productivity from death and disability, and the effects of costs on household impoverishment .
The body is made up of many types of cells. Normally, cells grow and divide to produce new cells in a controlled and orderly manner. Sometimes, however, new cells continue to be produced when they are not needed. As a result, a mass of extra tissue called a tumor may develop. Tumor can be cancerous (malignant) or noncancerous (benign). A benign tumor is a known as tumor cell mass that does not fragment and spread beyond its original area of growth. Generally, benign impact on the body is not harmful and easy to be treated. Benign tumor can be harmful by growing large enough to interfere with normal body functions. Malignant tumors are nonencapsulated growths of tumor cells that are harmful; they have no wall or clear-cut border and may spread or invade other parts of the body normal tissue. In course of time, the cancer cells interfere with the normal functioning of organs via lymph or blood vessels [4, 5]. This stage is known as secondary tumor. The development of the cancerous cells is very complex and involves interaction of the various cells. In human anatomy, the immune system is triggered to quest and eradicate the tumor cells when they are detectable as nonself [6–9].
Researches are going in two directions to fight with cancer: one is experimental and the other is theoretical, that is, mathematical. Both experimentalists and theoreticians join hands to get rid of deadly disease, that is, cancer [10–13]. The studies done in  provide a comprehensive overview of the different approaches used in the modelling of the tumor-immune system interaction dynamics. A number of researches had been done so far which we referred to as Kuznetsov et al. , Saleem and Agrawal , De Pillis et al. [16, 17], Zhivkov and Waniewski , Galach , Babbs , and Ahmad and Kaur , by using ordinary differential equations, and tried to investigate the interaction among tumor cells and different type of immune cells. Kuznetsov et al.  study nonlinear dynamics of immunogenic tumors and manifest a number of phenomena as immunostimulation of tumor growth and sneaking through and formation of a dormant tumor. Galach  concluded in his simplified model that tumor was in dormant state, but in time delay model there was recurring of tumor in the presence of immune cells.
For understanding the interaction between tumor and immune cells, several researchers used the concept of prey-predator [2, 7, 15–17, 19, 22, 23]. There is a difference between the classical prey-predator model and tumor-immune system prey-predator model; it is that in the latter model survival of the immune population does not depend on the number of the preys (tumor cells). The immune cells play the role of the predator while tumor cells of prey. The cellular immune response identifies and eliminates the tumor cells from the host because tumor cells produce some antigens on its cell surface. The strength of the immune response depends on the tumor antigenicity [8, 24–26]. The cellular response is carried by T lymphocytes. These T helper cells cannot kill tumor cells, but they send urgent biochemical signals to a special type of effector cells called cytotoxic T lymphocytes. These effector cells eliminate cancerous cells by mounting a cytotoxic reaction that lyses their target [7, 19, 27].
Cancer self-remission model using stochastic approach is studied by Sarkar and Banerjee . El-Gohary  extended the work of Sarkar and Banerjee  by providing an optimal control strategy that turned unstable steady states into asymptotically stable. The goal of the present paper is to modify the existing model of El-Gohary  by introducing a function called Michaelis-Menten. When the resting cells bind to tumor cells, they release cytokines, which act to increase its affinity. We introduce Michaelis-Menten function in this paper and we observe that the introduction of this particular function helps us to get stability as the growth rate of resting cells increases.
2. Mathematical Model
In this section, we considered tumor progression and regression as a prey-predator like system. The predator is the immune system which slaughters the tumor cells (prey). In most of the mathematical models of the tumor-immune system, the response of the immune system is considered as a single population of cells, namely, effector cell [2, 7, 16, 19], which perform the task of destroying cancer cells. This simplifying assumption allows decreasing the complexity of the dynamics of the immune system.
The predator, that is, immune system, is eradicating tumor cells in two stages: one is hunting cells and another is resting cells. Here, we are considering that hunting cells can slaughter tumor cells, but resting cells cannot. The cellular immune response identifies and eliminates the tumor cells from the host because tumor cells produce some antigens on its outer surface. The strength of the immune response depends on the tumor antigenicity [8, 24–26]. The cellular response of the immune system is carried by T lymphocytes. During maturation, T cells surface contains specialized antibody like receptors that see fragments of antigens on the surface of tumor cells. In most of the cases, T cells can recognize only antigen that is bound to a cell membrane protein called major histocompatibility complex (MHC) molecule. MHC molecule is a protein recognized by resting T cells, which distinguish between self and nonself. Resting T cells engulf the tumor cells and then produce various growth factors known collectively as cytokines, but they cannot kill tumor cells. Cytokines are chemical messenger switches which turn on the cytotoxic T lymphocytes (hunting cells). In contrast to the resting T cell, the cytotoxic T lymphocytes generally not only secrete many cytokines but also eliminate tumor cells by mounting a cytotoxic reaction that lyses their target [7, 19, 27].
Considering the above biological mechanism, we have produced a mathematical model of tumor development in immune response. The model involves certain assumptions as follows:(i)logistic growth function is assumed for the growth of tumor cells in the absence of hunting CTL cells;(ii)the tumor cells and hunting cells are being eradicated at a rate proportional to the densities of tumor cells and hunting predator cells according to the law of mass action;(iii)the resting predator cells are converted to the hunting cells, either by direct contact with them or by contact with a fast diffusing substance (cytokines) produced by the hunting cells;(iv)resting cells also follow logistic growth in absence of tumor cells;(v)once a hunting T cell has been converted, it will never return to the resting stage;(vi)resting cells also were stimulated due to the presence of tumor cells, and this is considered by the Michaelis-Menten function.
Let be the concentration of tumor cell in the given physiologic space at time , the concentration of hunting T cells, and the concentration of resting T cells. The model governing the interaction between , , and is given by the following system of equations: where is conversion rate of normal cells into tumor cells, is growth rate of tumor cells, is maximum carrying capacity of tumor cells, is rate of killing of tumor cells by hunting cells, is conversion rate of resting cells into hunting cells, is apoptosis rate of hunting cells, is rate of killing of hunting cells by tumor cells, is growth rate of resting cells, is maximum carrying capacity of resting cells, is apoptosis rate of resting cells, is proliferation rate of resting cells, and is half-saturation for proliferation term.
Define the following dimensionless variables to reduce the number of the system parameters: Incorporating the dimension variables in the system of (1) and suppressing the star for our convenience, we get where
2.1. Determination of Equilibrium
To find the equilibrium point, we put the time rate of change as zero. Therefore, the system of (3)–(5) reduces to
The solution of (7) in the space gives us the equilibrium points and coexistence equilibrium . We observe the following.(i) The equilibrium point lies on the boundary on the positive octant; that is, This will always exist. (ii)The equilibrium point lies in a - plane; that is, where exists only when for all positive parameters.(iii)The equilibrium point is in the - plane, that is, In this case, we see that is negative, so this point may not be taken into account.(iv)Coexisting equilibrium point is . Furthermore, eliminating and between the system of (7), we get a polynomial of degree three; that is, where The coefficients are not numbers here; they are functions of some model parameters like , , , and so forth, so the zero of this polynomial may not be obtained directly. Here, is negative clearly while other coefficients are not defined in the sign, but the model parameters determine the sign of these coefficients.
2.2. Local Stability of the Prospective Equilibrium Points
The local asymptotic stability of each nonnegative equilibrium point has been studied by computing the variational matrix and finding the eigenvalues. For the stability of equilibrium points, the real parts of eigenvalues of variational matrix must be negative.
The variational matrix due to linearization of the system equations (3)–(5) at the arbitrary equilibrium point (, , ) is given by
The local dynamical behavior of equilibrium points are investigated and obtained results by computing the variational matrices corresponding to each equilibrium point. The local asymptotic stability for each equilibrium point has been analyzed as follows(i)The variational matrix at is
The eigenvalues of the are (<0), (<0), and (>0). Observing the sign of eigenvalues, we can state the following theorem.
Theorem 1. is always a saddle point if exists with locally unstable manifold along direction and with the local stable manifold in - plane.
(ii) The variational matrix at is
The eigenvalues of this matrix are (<0), = , and (<0 is the existence condition of ). So, we state the following theorem.
Theorem 2. If , then is stable asymptotically in - plane, but if does not hold, then it will be saddle point with stable manifold in - plane and with unstable manifold in direction.
(iii) The variational matrix at is
As , therefore, the eigenvalues of are obtained as follows: where and . Consider Define , , and + − . Hence, we get + + = 0. Now, following Routh-Hurwitz criteria, the ’s are negative if , , and . We considered the following cases:(a) . Which is only possible when and . When(i); ;(ii); ;(b). As and . Therefore, for ;(c) − . Hence we arrive at a conclusion.
Theorem 3. If , , and , then is locally asymptotically stable in -- plane.
2.3. Global Stability of the System
The system considered in this paper is determined by three nonnegative equilibrium points , , and . Out of these points, is unstable; therefore, our aim is to verify the global stability of the system through and . Hence, we check the global stability of these points as follows.
Theorem 4. If the equilibrium point is locally asymptotically stable in the interior of a positive quadrant of - plane, then it will be globally asymptotically stable there.
Proof. Define Dulac function . Using the system of (3) and (5), let
Now, It is observed here that does not change its sign in the positive quadrant of - plane. By Bendixson-Dulac criterion, there is no limit cycle or homoclinic connection observed in the positive quadrant of the - plane. Hence, if is locally asymptotically stable, then it will be globally asymptotically stable in the interior of a positive quadrant of - plane.
Theorem 5. If the equilibrium point is locally asymptotically stable in the interior of the positive octant of space then it will be globally asymptotically stable there.
Proof. Consider the following Lyapunov function around :
Differentiating with respect to , we get
Using (3)–(5) in , we have
Thus, is a quadratic form which can be expressed as , where = and is a symmetric matrix given by
with , , , , , and .
is negative definite if , , and . Hence, the is a Lyapunov function with respect to .
3. Simulation and Discussion
This section is devoted to the study of a mathematical model presented by the system of (3)–(5). Using Matlab solver and assuming the values of some parameters randomly, we try to reach the conclusions. Figures 1 and 2 have been obtained for different values of by taking all others a’s fixed to study the behavior of the model.
From a numerical simulation, we discovered that by introducing the simulating function in resting cell population due to the presence of tumor; while increase the growth rate of resting cells can control the progression of tumor. In Figure 1, we compare the temporal growth of the number of tumor cells, hunting cells, and resting cells for different values of the parameter . The level of interaction is maximum according to Figure 1(a) when which shows a sustained oscillation for different populations. Observing Figure 1(a) and its phase portrait Figure 2(a), we exhibit a stage of “recurring” tumor [19, 28, 29]. The cyclic fluctuation reduces in different populations while increasing the value of . The sustains oscillation give a way to a stable spiral with quick damping, which leads to the persistent tumor, that could be described as dormant (A small stable tumor does not change in size is referred as dormant tumor). Leading towards the persistent tumor, that could be described as dormant (a small stable tumor which does not change in size is referred to as dormant tumor [19, 29]). This state is desirable clinical condition since the tumor growth is blocked. The phase portrait Figure 2(a) exhibits a trend that the recovery phase is weaker than the tumor phase. It has also been noticed that the sustained oscillation phase tends to damped oscillation phase for . According to the damped phase in Figure 1(b), we see that hunting cells attack tumor cells together with resting cells for a short span of time and later hunting cells do not interfere either with resting cells or with tumor cells. From the rest of Figures 1(c), 1(d) and 2(b) to 2(d), it has been observed that, on increasing , the growth rate of resting cells, the tumor phase tends to recovery phase. Also, it has been seen that the system has not been reached to complete recovery phase; that is, tumor is not eliminated.
An interaction between tumor cells as prey and immune system which consists of resting and hunting cells as predator has been studied. We analyze the model with regard to local and global stability of equilibrium points. Global stability of steady state is established using Bendixson-Dulac criteria and coexisting steady state is done by Lyapunov function. Using the Michealis-Menten function approach, it has been observed that as , the growth of resting cells, increases, the behavior of tumor growth persists from recurring state to a dormant state. A better understanding of the influence of the immune system on dormancy could lead to the development of immunological therapies in order to prevent the progression of the tumour and then the switch from the state of dormancy to tumour growth. Our model may lead to novel detection techniques in the clinic and therapeutic options to prevent deadly tumor recurrence and metastasis.
Conflict of Interests
The authors have no conflict of interests related to the conduct and reporting of this research. They have no financial relationships with any organization.
Cancer Facts & Figures, American Cancer Society, 2014.
F. A. Rihan, M. Safan, M. A. Abdeen, and D. Abdel Rahman, “Qualitative and computational analysis of a mathematical model for tumor-immune interactions,” Journal of Applied Mathematics, vol. 2012, Article ID 475720, 19 pages, 2012.View at: Publisher Site | Google Scholar
“Cancer deaths in India likely to touch 7 lakh by 2015,” http://archive.indianexpress.com/news/-cancer-deaths-in-india-likely-to-touch-7-lakh-by-2015-/1068972/.View at: Google Scholar
A. K. Abbas, A. H. Lichtman, and J. S. Pober, Celluar and Molecular Immunology, WB Saunders, Philadelphia, Pa, USA, 3rd edition, 1997.
I. M. Roitt, J. Brostoff, and D. K. Male, Immunology, Mosby, St.Louis, Miss, USA, 5th edition, 1998.
F. Nani and H. I. Freedman, “A mathematical model of cancer treatment by immunotherapy,” Mathematical Biosciences, vol. 163, no. 2, pp. 159–199, 2000.View at: Publisher Site | Google Scholar | MathSciNet
V. A. Kuznetsov, I. A. Makalkin, M. A. Taylor, and A. S. Perelson, “Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis,” Bulletin of Mathematical Biology, vol. 56, no. 2, pp. 295–321, 1994.View at: Publisher Site | Google Scholar
M. Robertson-Tessi, A. El-Kareh, and A. Goriely, “A mathematical model of tumor-immune interactions,” Journal of Theoretical Biology, vol. 294, pp. 56–73, 2012.View at: Publisher Site | Google Scholar
R. Eftimie, J. L. Bramson, and D. J. D. Earn, “Interactions between the immune system and cancer: a brief review of non-spatial mathematical models,” Bulletin of Mathematical Biology, vol. 73, no. 1, pp. 2–32, 2011.View at: Publisher Site | Google Scholar
C. L. Jorcyk, M. Kolev, K. Tawara, and B. Zubik-Kowal, “Experimental versus numerical data for breast cancer progression,” Nonlinear Analysis. Real World Applications, vol. 13, no. 1, pp. 78–84, 2012.View at: Publisher Site | Google Scholar | MathSciNet
T. L. Jackson, “A mathematical investigation of the multiple pathways to recurrent prostate cancer: comparison with experimental data,” Neoplasia, vol. 6, no. 6, pp. 697–704, 2004.View at: Publisher Site | Google Scholar
C. Bianca, F. Chiacchio, F. Pappalardo, and M. Pennisi, “Mathematical modeling of the immune system recognition to mammary carcinoma antigen,” BMC Bioinformatics, vol. 13, no. 17, article S21, 2012.View at: Publisher Site | Google Scholar
L. G. De Pillis, A. Radunskaya, and K. J. Bathe, “A mathematical model of immune response to tumor invasion,” in Proceedings of the 2nd MIT Conference on Computational Fluid and Solid Mechanics, K. J. Bathe, Ed., Computational Fluid and Solid Mechanics, 2003.View at: Google Scholar
A. L. Woelke, M. S. Murgueitio, and R. Preissner, “Theoretical modeling techniques and their impact on tumor immunology,” Clinical and Developmental Immunology, vol. 2010, Article ID 271794, 11 pages, 2010.View at: Publisher Site | Google Scholar
M. Saleem and T. Agrawal, “Chaos in a tumor growth model with delayed responses of the immune system,” Journal of Applied Mathematics, vol. 2012, Article ID 891095, 16 pages, 2012.View at: Publisher Site | Google Scholar | MathSciNet
L. G. de Pillis, A. E. Radunskaya, and C. L. Wiseman, “A validated mathematical model of cell-mediated immune response to tumor growth,” Cancer Research, vol. 65, no. 17, pp. 7950–7958, 2005.View at: Google Scholar
L. G. De Pillis, W. Gu, and A. E. Radunskaya, “Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations,” Journal of Theoretical Biology, vol. 238, no. 4, pp. 841–862, 2006.View at: Publisher Site | Google Scholar
P. Zhivkov and J. Waniewski, “Modelling tumour-immunity interactions with different stimulation functions,” International Journal of Applied Mathematics and Computer Science, vol. 13, no. 3, pp. 307–315, 2003.View at: Google Scholar | MathSciNet
M. Galach, “Dynamics of the tumor-immune system competition—the effect of time delay,” International Journal of Applied Mathematics and Computer Science, vol. 13, no. 3, pp. 395–406, 2003.View at: Google Scholar | MathSciNet
C. F. Babbs, “Predicting success or failure of immunotherapy for cancer: insights from a clinically applicable mathematical model,” American Journal of Cancer Research, vol. 2, pp. 204–213, 2012.View at: Google Scholar
Gurpreet Kaur and Naseem Ahmad, “A study of population dynamics of normal and immune cells in presence of tumor cells,” International Journal of Scientific & Engineering Research, vol. 4, no. 4, pp. 770–775, 2013.View at: Google Scholar
R. R. Sarkar and S. Banerjee, “Cancer self remission and tumor stability—a stochastic approach,” Mathematical Biosciences, vol. 196, no. 1, pp. 65–81, 2005.View at: Publisher Site | Google Scholar
A. El-Gohary, “Chaos and optimal control of cancer self-remission and tumor system steady states,” Chaos, Solitons and Fractals, vol. 37, no. 5, pp. 1305–1316, 2008.View at: Publisher Site | Google Scholar | Zentralblatt MATH
A. Albert, M. Freedman, and A. S. Perelson, “Tumors and the immune system: the effects of a tumor growth modulator,” Mathematical Biosciences, vol. 50, no. 1-2, pp. 25–58, 1980.View at: Publisher Site | Google Scholar
B. Joshi, X. Wang, S. Banerjee, H. Tian, A. Matzavinos, and M. A. J. Chaplain, “On immunotherapies and cancer vaccination protocols: a mathematical modelling approach,” Journal of Theoretical Biology, vol. 259, no. 4, pp. 820–827, 2009.View at: Publisher Site | Google Scholar
A. Usman, T. Jackson, and C. Cunningham, “Application of the mathematical model of tumor-immune interactions for IL-2 Adoptive Immunotherapy to studies on patients with Metastatic Melanoma or Renal Cell Cancer,” The Rose-Hulman Undergraduate Math Journal, vol. 6, no. 2, p. 23, 2005.View at: Google Scholar
W. M. Yokoyama, S. Kim, and A. R. French, “The dynamic life of natural killer cells,” Annual Review of Immunology, vol. 22, pp. 405–429, 2004.View at: Publisher Site | Google Scholar
D. Kirschner and J. C. Panetta, “Modeling immunotherapy of the tumo—immune interaction,” Journal of Mathematical Biology, vol. 37, no. 3, pp. 235–252, 1998.View at: Publisher Site | Google Scholar
R. R. Sarkar and S. Banerjee, “A time delay model for control of malignant tumor growth,” in Proceedings of the National Conference on Nonlinear Systems and Dynamics, pp. 1–4, 2006.View at: Google Scholar