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.

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:

logistic growth function is assumed for the growth of tumor cells in the absence of hunting CTL cells;

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;

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;

resting cells also follow logistic growth in absence of tumor cells;

once a hunting T cell has been converted, it will never return to the resting stage;

resting cells also were stimulated due to the presence of tumor cells, and this is considered by the Michaelis-Menten function.

Let T(t) be the concentration of tumor cell in the given physiologic space at time t, H(t) the concentration of hunting T cells, and R(t) the concentration of resting T cells. The model governing the interaction between T(t), H(t), and R(t) is given by the following system of equations: where q is conversion rate of normal cells into tumor cells, r 1 is growth rate of tumor cells, k 1 is maximum carrying capacity of tumor cells, α 1 is rate of killing of tumor cells by hunting cells, β is conversion rate of resting cells into hunting cells, d 1 is apoptosis rate of hunting cells, α 2 is rate of killing of hunting cells by tumor cells, r 2 is growth rate of resting cells, k 2 is maximum carrying capacity of resting cells, d 2 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 THR space gives us the equilibrium points and coexistence equilibrium E 4(T 4, H 4, R 4). We observe the following.The coefficients are not numbers here; they are functions of some model parameters like a 0, a 1, a 2, and so forth, so the zero of this polynomial may not be obtained directly. Here, C 0 is negative clearly while other coefficients are not defined in the sign, but the model parameters determine the sign of these coefficients.

The equilibrium point lies on the boundary on the positive octant; that is,

This will always exist.

The equilibrium point E 2 lies in a T-R plane; that is,

where E 2 exists only when a 5 + (a 8 T 2/(T 2 + K)) > a 7 for all positive parameters.

The equilibrium point E 3 is in the T-H plane, that is, In this case, we see that T 3 is negative, so this point may not be taken into account.

Coexisting equilibrium point is E 4(T 4, H 4, R 4).

Furthermore, eliminating H and R between the system of (7), we get a polynomial of degree three; that is,

where

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 E 0 (T, H, R) 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

The variational matrix at E 1 is

The eigenvalues of the V 1 are (<0), (<0), and (>0). Observing the sign of eigenvalues, we can state the following theorem.

Theorem 1 .

E 1 is always a saddle point if E 2 exists with locally unstable manifold along R direction and with the local stable manifold in T-H plane.

(ii) The variational matrix at E 2 is

The eigenvalues of this matrix are (<0), λ 2 = a 2 R 2 − a 3 − a 4 T 2, and λ 3 = −(a 5 − a 7 + (a 8 T 2/(T 2 + K))) (<0 is the existence condition of E 2). So, we state the following theorem.

Theorem 2 .

If λ 2 < 0, then E 2 is stable asymptotically in T-R plane, but if λ 2 < 0 does not hold, then it will be saddle point with stable manifold in T-R plane and with unstable manifold in H direction.

(iii) The variational matrix at E 4 is

As R 4 = (a 3 + a 4 T 4)/a 2, therefore, the eigenvalues of V 4 are obtained as follows: where A = a 1 − 2a 1 T 4 − H 4 and B = a 5 − 2a 5 R 4 − a 6 H 4 − a 7 + (a 8 T 4/(T 4 + K)). Consider Define A 1 = −(A + B), A 2 = AB + a 2 a 6 R 4 H 4 − a 4 H 4 T 4, and A 3 = a 4 BH 4 T 4  +  (a 2 a 8 KH 4 T 4 R 4/(T 4+K)2)  −  a 2 a 6 AH 4 R 4. Hence, we get λ 3  +  A 1 λ 2  +  A 2 λ + A 3 = 0. Now, following Routh-Hurwitz criteria, the λ's are negative if A 1 > 0, A 3 > 0, and A 1 A 2 − A 3 > 0. We considered the following cases:

A 1 > 0⇒−(A + B) > 0.

Which is only possible when A < 0 and B < 0. When

A < 0⇒a 1 − 2a 1 T 4 − H 4 < 0  that  is;  a 1 < 2a 1 T 4 + H 4;

B < 0⇒a 5 − 2a 5 R 4 − a 6 H 4 − a 7 + (a 8 T 4/(T 4 + K)) < 0  that  is;  a 5 + (a 8 T 4/(T 4 + K)) < 2a 5 R 4 + a 6 H 4 + a 7;

A 3 > 0⇒a 4 BH 4 T 4 + (a 2 a 8 KH 4 T 4 R 4/(T 4+K)2) − a 2 a 6 AH 4 R 4 > 0.

As A < 0 and B < 0. Therefore, A 3 > 0 for a 4 BH 4 T 4 + (a 2 a 8 KH 4 T 4 R 4/(T 4+K)2) > a 2 a 6 AH 4 R 4;

A 1 A 2 − A 3 > 0⇒−(A + B)(AB + a 2 a 6 R 4 H 4 − a 4 H 4 T 4)  −  (a 4 BH 4 T 4 + (a 2 a 8 KH 4 T 4 R 4/(T 4+K)2) − a 2 a 6 AH 4 R 4)  >  0⇒a 4 AH 4 T 4 > A 2 B + AB 2 + a 2 a 6 BR 4 H 4 + (a 2 a 8 KH 4 R 4 T 4/(T 4+K)2). Hence we arrive at a conclusion.

Theorem 3 .

If A 1 > 0, A 3 > 0, and A 1 A 2 − A 3 > 0, then E 4 is locally asymptotically stable in T-H-R plane.

2.3. Global Stability of the System

The system considered in this paper is determined by three nonnegative equilibrium points E 1, E 2, and E 4. Out of these points, E 1 is unstable; therefore, our aim is to verify the global stability of the system through E 2 and E 4. Hence, we check the global stability of these points as follows.

Theorem 4 .

If the equilibrium point E 2 is locally asymptotically stable in the interior of a positive quadrant of T-R plane, then it will be globally asymptotically stable there.

Proof

Define Dulac function H 1(T, R) = (1/TR); T, R > 0. Using the system of (3) and (5), let

Now, It is observed here that Div(H 1 f, H 1 g) does not change its sign in the positive quadrant of T-R plane. By Bendixson-Dulac criterion, there is no limit cycle or homoclinic connection observed in the positive quadrant of the T-R plane. Hence, if E 2 is locally asymptotically stable, then it will be globally asymptotically stable in the interior of a positive quadrant of T-R plane.

Theorem 5 .

If the equilibrium point E 4 is locally asymptotically stable in the interior of the positive octant of THR space then it will be globally asymptotically stable there.

Consider the following Lyapunov function around E 4(T 4, H 4, R 4): Differentiating V with respect to t, we get Using (3)–(5) in , we have Thus, is a quadratic form which can be expressed as , where X = (T − T 4, H − H 4, R − R 4) and A is a symmetric matrix given by with a 11 = (a 1 + (1/TT 4)),  a 12 = (1 + a 4)/2,  a 13 = −(a 8 K/2(T + K)(T 4 + K)), a 22 = 0, a 23 = (a 6 − a 2)/2, and a 33 = a 5.

is negative definite if a 11 > 0, a 33 > 0, and a 13 2 − a 11 a 33 ≤ 0. Hence, the V is a Lyapunov function with respect to E 4.

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 a 5 by taking all others a's fixed to study the behavior of the model.

Figure 1

Tumor cells, hunting cells, and resting cells with time for a 1 = 1.82, a 2 = 0.239, a 3 = 0.2, a 4 = 0.04, a 6 = 0.5, a 7 = 0.01, a 8 = 2, k = 1, (a) a 5 = 0.0191, (b) a 5 = 0.0691, (c), a 5 = 0.191, and (d) a 5 = 0.6291.

Figure 2

Phase portraits corresponding to the system (3)–(5) for the following parameter values: a 1 = 1.82, a 2 = 0.239, a 3 = 0.2, a 4 = 0.04, a 6 = 0.5, a 7 = 0.01, a 8 = 2, k = 1, (a) a 5 = 0.0191, (b) a 5 = 0.0691, (c) a 5 = 0.191, and (d) a 5 = 0.6291.

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 a 5. The level of interaction is maximum according to Figure 1(a) when a 5 = 0.0191 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 a 5. 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 a 5 = 0.046039. 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 a 5, 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 E 2 is established using Bendixson-Dulac criteria and coexisting steady state E 4 is done by Lyapunov function. Using the Michealis-Menten function approach, it has been observed that as a 5, 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.