On Study of Immune Response to Tumor Cells in Prey-Predator System

Abstract

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.

Gurpreet Kaur 1 and Naseem Ahmad 2

Academic Editor:Francesco Pappalardo

1, CIRBSc., Jamia Millia Islamia, Jamia Nagar, New Delhi 110025, India
2, Department of Mathematics, Jamia Millia Islamia, Jamia Nagar, New Delhi 110025, India

Received 23 April 2014; Revised 9 June 2014; Accepted 9 June 2014; 22 September 2014

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

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. [1]. Not much is known about its mechanisms of establishment and destruction. It is the second fatal disease after the cardiovascular diseases [2]. 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 [3].

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 [14] 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. [7], Saleem and Agrawal [15], De Pillis et al. [16, 17], Zhivkov and Waniewski [18], Galach [19], Babbs [20], and Ahmad and Kaur [21], by using ordinary differential equations, and tried to investigate the interaction among tumor cells and different type of immune cells. Kuznetsov et al. [7] 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 [19] 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 [22]. El-Gohary [23] extended the work of Sarkar and Banerjee [22] 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 [23] 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 [figure omitted; refer to PDF] be the concentration of tumor cell in the given physiologic space at time [figure omitted; refer to PDF] , [figure omitted; refer to PDF] the concentration of hunting T cells, and [figure omitted; refer to PDF] the concentration of resting T cells. The model governing the interaction between [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] is given by the following system of equations: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is conversion rate of normal cells into tumor cells, [figure omitted; refer to PDF] is growth rate of tumor cells, [figure omitted; refer to PDF] is maximum carrying capacity of tumor cells, [figure omitted; refer to PDF] is rate of killing of tumor cells by hunting cells, [figure omitted; refer to PDF] is conversion rate of resting cells into hunting cells, [figure omitted; refer to PDF] is apoptosis rate of hunting cells, [figure omitted; refer to PDF] is rate of killing of hunting cells by tumor cells, [figure omitted; refer to PDF] is growth rate of resting cells, [figure omitted; refer to PDF] is maximum carrying capacity of resting cells, [figure omitted; refer to PDF] is apoptosis rate of resting cells, [figure omitted; refer to PDF] is proliferation rate of resting cells, and [figure omitted; refer to PDF] is half-saturation for proliferation term.

Define the following dimensionless variables to reduce the number of the system parameters: [figure omitted; refer to PDF] Incorporating the dimension variables in the system of (1) and suppressing the star [figure omitted; refer to PDF] for our convenience, we get [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF]

The solution of (7) in the [figure omitted; refer to PDF] space gives us the equilibrium points and coexistence equilibrium [figure omitted; refer to PDF] . We observe the following.

(i) The equilibrium point lies on the boundary on the positive octant; that is, [figure omitted; refer to PDF]

: This will always exist.

(ii) The equilibrium point [figure omitted; refer to PDF] lies in a [figure omitted; refer to PDF] - [figure omitted; refer to PDF] plane; that is, [figure omitted; refer to PDF]

: where [figure omitted; refer to PDF] exists only when [figure omitted; refer to PDF] for all positive parameters.

(iii): The equilibrium point [figure omitted; refer to PDF] is in the [figure omitted; refer to PDF] - [figure omitted; refer to PDF] plane, that is, [figure omitted; refer to PDF] In this case, we see that [figure omitted; refer to PDF] is negative, so this point may not be taken into account.

(iv) Coexisting equilibrium point is [figure omitted; refer to PDF] .

: Furthermore, eliminating [figure omitted; refer to PDF] and [figure omitted; refer to PDF] between the system of (7), we get a polynomial of degree three; that is, [figure omitted; refer to PDF]

: where [figure omitted; refer to PDF]

The coefficients are not numbers here; they are functions of some model parameters like [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and so forth, so the zero of this polynomial may not be obtained directly. Here, [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] ) is given by [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] is [figure omitted; refer to PDF]

The eigenvalues of the [figure omitted; refer to PDF] are [figure omitted; refer to PDF] (<0), [figure omitted; refer to PDF] (<0), and [figure omitted; refer to PDF] (>0). Observing the sign of eigenvalues, we can state the following theorem.

Theorem 1.

[figure omitted; refer to PDF] is always a saddle point if [figure omitted; refer to PDF] exists with locally unstable manifold along [figure omitted; refer to PDF] direction and with the local stable manifold in [figure omitted; refer to PDF] - [figure omitted; refer to PDF] plane.

(ii) The variational matrix at [figure omitted; refer to PDF] is [figure omitted; refer to PDF]

The eigenvalues of this matrix are [figure omitted; refer to PDF] (<0), [figure omitted; refer to PDF] = [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] (<0 is the existence condition of [figure omitted; refer to PDF] ). So, we state the following theorem.

Theorem 2.

If [figure omitted; refer to PDF] , then [figure omitted; refer to PDF] is stable asymptotically in [figure omitted; refer to PDF] - [figure omitted; refer to PDF] plane, but if [figure omitted; refer to PDF] does not hold, then it will be saddle point with stable manifold in [figure omitted; refer to PDF] - [figure omitted; refer to PDF] plane and with unstable manifold in [figure omitted; refer to PDF] direction.

(iii) The variational matrix at [figure omitted; refer to PDF] is [figure omitted; refer to PDF]

As [figure omitted; refer to PDF] , therefore, the eigenvalues of [figure omitted; refer to PDF] are obtained as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Consider [figure omitted; refer to PDF] Define [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] + [figure omitted; refer to PDF] - [figure omitted; refer to PDF] . Hence, we get [figure omitted; refer to PDF] + [figure omitted; refer to PDF] + [figure omitted; refer to PDF] = 0. Now, following Routh-Hurwitz criteria, the [figure omitted; refer to PDF] 's are negative if [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . We considered the following cases:

(a) [figure omitted; refer to PDF] [figure omitted; refer to PDF] .

: Which is only possible when [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . When

(i) [figure omitted; refer to PDF] ; [figure omitted; refer to PDF] ;

(ii) [figure omitted; refer to PDF] ; [figure omitted; refer to PDF] ;

(b) [figure omitted; refer to PDF] .

: As [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Therefore, [figure omitted; refer to PDF] for [figure omitted; refer to PDF] ;

(c) [figure omitted; refer to PDF] - [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] . Hence we arrive at a conclusion.

Theorem 3.

If [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , then [figure omitted; refer to PDF] is locally asymptotically stable in [figure omitted; refer to PDF] - [figure omitted; refer to PDF] - [figure omitted; refer to PDF] plane.

2.3. Global Stability of the System

The system considered in this paper is determined by three nonnegative equilibrium points [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Out of these points, [figure omitted; refer to PDF] is unstable; therefore, our aim is to verify the global stability of the system through [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Hence, we check the global stability of these points as follows.

Theorem 4.

If the equilibrium point [figure omitted; refer to PDF] is locally asymptotically stable in the interior of a positive quadrant of [figure omitted; refer to PDF] - [figure omitted; refer to PDF] plane, then it will be globally asymptotically stable there.

Proof.

Define Dulac function [figure omitted; refer to PDF] . Using the system of (3) and (5), let [figure omitted; refer to PDF]

Now, [figure omitted; refer to PDF] It is observed here that [figure omitted; refer to PDF] does not change its sign in the positive quadrant of [figure omitted; refer to PDF] - [figure omitted; refer to PDF] plane. By Bendixson-Dulac criterion, there is no limit cycle or homoclinic connection observed in the positive quadrant of the [figure omitted; refer to PDF] - [figure omitted; refer to PDF] plane. Hence, if [figure omitted; refer to PDF] is locally asymptotically stable, then it will be globally asymptotically stable in the interior of a positive quadrant of [figure omitted; refer to PDF] - [figure omitted; refer to PDF] plane.

Theorem 5.

If the equilibrium point [figure omitted; refer to PDF] is locally asymptotically stable in the interior of the positive octant of [figure omitted; refer to PDF] space then it will be globally asymptotically stable there.

Proof.

Consider the following Lyapunov function around [figure omitted; refer to PDF] : [figure omitted; refer to PDF] Differentiating [figure omitted; refer to PDF] with respect to [figure omitted; refer to PDF] , we get [figure omitted; refer to PDF] Using (3)-(5) in [figure omitted; refer to PDF] , we have [figure omitted; refer to PDF] Thus, [figure omitted; refer to PDF] is a quadratic form which can be expressed as [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] = [figure omitted; refer to PDF] and [figure omitted; refer to PDF] is a symmetric matrix given by [figure omitted; refer to PDF] with [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] .

[figure omitted; refer to PDF] is negative definite if [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Hence, the [figure omitted; refer to PDF] is a Lyapunov function with respect to [figure omitted; refer to PDF] .

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 [figure omitted; refer to PDF] by taking all others a's fixed to study the behavior of the model.

Tumor cells, hunting cells, and resting cells with time for [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , (a) [figure omitted; refer to PDF] , (b) [figure omitted; refer to PDF] , (c), [figure omitted; refer to PDF] , and (d) [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

Phase portraits corresponding to the system (3)-(5) for the following parameter values: [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , (a) [figure omitted; refer to PDF] , (b) [figure omitted; refer to PDF] , (c) [figure omitted; refer to PDF] , and (d) [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] . The level of interaction is maximum according to Figure 1(a) when [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] . 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 [figure omitted; refer to PDF] , 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.

4. Conclusion

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 [figure omitted; refer to PDF] is established using Bendixson-Dulac criteria and coexisting steady state [figure omitted; refer to PDF] is done by Lyapunov function. Using the Michealis-Menten function approach, it has been observed that as [figure omitted; refer to PDF] , 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.

References

[1] Cancer Facts & Figures , American Cancer Society, 2014.

[2] F. A. Rihan, M. Safan, M. A. Abdeen, D. Abdel Rahman, "Qualitative and computational analysis of a mathematical model for tumor-immune interactions," Journal of Applied Mathematics , vol. 2012, 2012.

[3] "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/

[4] A. K. Abbas, A. H. Lichtman, J. S. Pober Celluar and Molecular Immunology , WB Saunders, Philadelphia, Pa, USA, 1997., 3rd.

[5] I. M. Roitt, J. Brostoff, D. K. Male Immunology , Mosby, St.Louis, Miss, USA, 1998., 5th.

[6] F. Nani, H. I. Freedman, "A mathematical model of cancer treatment by immunotherapy," Mathematical Biosciences , vol. 163, no. 2, pp. 159-199, 2000.

[7] V. A. Kuznetsov, I. A. Makalkin, M. A. Taylor, 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.

[8] M. Robertson-Tessi, A. El-Kareh, A. Goriely, "A mathematical model of tumor-immune interactions," Journal of Theoretical Biology , vol. 294, pp. 56-73, 2012.

[9] R. Eftimie, J. L. Bramson, 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.

[10] C. L. Jorcyk, M. Kolev, K. Tawara, B. Zubik-Kowal, "Experimental versus numerical data for breast cancer progression," Nonlinear Analysis. Real World Applications , vol. 13, no. 1, pp. 78-84, 2012.

[11] 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.

[12] C. Bianca, F. Chiacchio, F. Pappalardo, M. Pennisi, "Mathematical modeling of the immune system recognition to mammary carcinoma antigen," BMC Bioinformatics , vol. 13, no. 17, article S21, 2012.

[13] L. G. De Pillis, A. Radunskaya, K. J. Bathe, 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, of Computational Fluid and Solid Mechanics, 2003.

[14] A. L. Woelke, M. S. Murgueitio, R. Preissner, "Theoretical modeling techniques and their impact on tumor immunology," Clinical and Developmental Immunology , vol. 2010, 2010.

[15] M. Saleem, T. Agrawal, "Chaos in a tumor growth model with delayed responses of the immune system," Journal of Applied Mathematics , vol. 2012, 2012.

[16] L. G. de Pillis, A. E. Radunskaya, 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.

[17] L. G. De Pillis, W. Gu, 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.

[18] P. Zhivkov, 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.

[19] 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.

[20] 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.

[21] Gurpreet Kaur, 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.

[22] R. R. Sarkar, S. Banerjee, "Cancer self remission and tumor stability-a stochastic approach," Mathematical Biosciences , vol. 196, no. 1, pp. 65-81, 2005.

[23] 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.

[24] A. Albert, M. Freedman, 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.

[25] B. Joshi, X. Wang, S. Banerjee, H. Tian, A. Matzavinos, 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.

[26] A. Usman, T. Jackson, 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, pp. 23, 2005.

[27] W. M. Yokoyama, S. Kim, A. R. French, "The dynamic life of natural killer cells," Annual Review of Immunology , vol. 22, pp. 405-429, 2004.

[28] D. Kirschner, J. C. Panetta, "Modeling immunotherapy of the tumo-immune interaction," Journal of Mathematical Biology , vol. 37, no. 3, pp. 235-252, 1998.

[29] R. R. Sarkar, 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.