Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2012, Article ID 475720, 19 pages
http://dx.doi.org/10.1155/2012/475720
Research Article

Qualitative and Computational Analysis of a Mathematical Model for Tumor-Immune Interactions

1Department of Mathematical Sciences, Faculty of Science, United Arab Emirates University, Al-Ain 17551, United Arab Emirates
2Department of Mathematics, Faculty of Science, Mansoura University, Mansoura 35516, Egypt
3Department of Mathematics, Faculty of Science, Helwan University, Cairo 11790, Egypt

Received 7 July 2011; Revised 15 October 2011; Accepted 28 October 2011

Academic Editor: Pedro Serranho

Copyright © 2012 F. A. Rihan et al. 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.

Abstract

We provide a family of ordinary and delay differential equations to model the dynamics of tumor-growth and immunotherapy interactions. We explore the effects of adoptive cellular immunotherapy on the model and describe under what circumstances the tumor can be eliminated. The possibility of clearing the tumor, with a strategy, is based on two parameters in the model: the rate of influx of the effector cells and the rate of influx of IL-2. The critical tumor-growth rate, below which endemic tumor does not exist, has been found. One can use the model to make predictions about tumor dormancy.

1. Introduction

Cancer is one of the most difficult diseases to be treated clinically, and one of the main causes of death. It is the second fatal disease after the cardiovascular diseases. The World Health Organization estimates that the annual cancer-induced mortality number exceeds six million people. Accordingly, the fight against cancer is of major public health interest. For this and other economy-related reasons, a great research effort is being devoted to understand the dynamics of cancer and to predict the impact of any changes on the system reactors. Hence, mathematical models are required to help design therapeutic strategies.

In cancer modeling, we have to care about the scaling problem, where the class of equations, used to describe the model, are to be determined. Indeed, there are three natural scales, which are connected to the different stages of the disease and have to be identified. The first is the subcellular (or molecular) scale, where we focus on studying the alterations in the genetic expressions of the genes contained in the nucleus of a cell, as a result of some special signals, which are received by the receptors on the cell surface and transmitted to the cell nucleus. The second is the cellular scale, which is an intermediate level between the molecular and the macroscopic scale. The third is the macroscopic scale, where we deal with heterogeneous tissues. In the heterogeneous tissues, some of the layers (the external proliferating layer, the intermediate layer, and the inner zone with necrotic cells) constituting the tumor may occur as islands, leading to a tumor comprised of multiple regions of necrosis, engulfed by tumor cells in a quiescent or proliferative state [1]. In case of macroscopic scale, we focus on the interaction between the tumor and normal cells (e.g., immune cells and blood vessels) in each of the three layers. For more details about description of the scaling problem and the passage from each scale to another, we refer to Bellomo et al. [1, 2].

A great research effort is being devoted to understand the interaction between the tumor cells and the immune system. Mathematical models, using ordinary, partial, and delay differential equations [3], play an important role in understanding the dynamics and tracking tumor and immune populations over time. Although the theoretical study of tumor immune dynamics has a long history [4, 5], the multifaceted nature of cancer requires sophisticated, nonlinear mathematical models to capture more realistic growth dynamics.

Many mathematical models have been proposed to model the interactions of cytotoxic T lymphocyte (CTL) response and the growth of an immunogenic tumor (see, e.g., [611]). The model by Kuznetsov et al. [7] takes into account the penetration of the tumor cells by the effector cells, which simultaneously causes the inactivation of effector cells. However, the model of Matzavinos et al. [9] describes the growth of a solid tumor in the presence of an immune system response, with special focus on the attack of tumor cells by the tumor-infiltrating cytotoxic lymphocytes (TICLs) in a small, multicellular tumor, without necrosis and at some stages prior to angiogenesis. The analysis shows that the TICLs can play an important role in the control of cancer dormancy.

The treatment of cancer is then one of the most challenging problems of modern medicine. The treatment should satisfy two basic conditions: first, it should destroy cancer cells in the entire body. Second, it should distinguish between cancerous and healthy cells. Other treatments such as surgery and/or chemoand radiotherapies have played key roles in treatment [12], but in many cases they do not represent a cure. Immunotherapy seems to be the method that best fulfils both of these requirements [7, 13, 14].

Numerous research papers have been made to explore the effects of the immune system in eliminating the tumor cells in the host, by stimulating the host's own immune response to kill cancer cells [15]. When tumor cells appear in a body, the immune system tries to identify and then eliminate them. Immunotherapy refers to the use of cytokines usually together with Adoptive Cellular Immunotherapy (ACI). Cytokines are protein hormones that mediate both natural and specific immunity. They are produced mainly by activated T cells (lymphocytes) during cellular-mediated immunity. Interleukin-2 (IL-2) is the main cytokine responsible for lymphocyte activation, growth, and differentiation. IL-2 has been shown to enhance Cytotoxic T Cells (CTL) activity at different disease stages. However, ACI refers to the injection of cultured immune cells that have antitumor reactivity into tumor bearing host. This interaction is analyzed and studied in various levels of biomathematical researches. They commonly focused on the models on ODEs over time. For example, in 1985, DeBoer et al. [4] suggested a mathematical model which contains eleven ordinary differential equations with five algebraic equations to describe antitumor response with IL-2 taken into account. A simple version of this model is proposed by Kirschner and Panetta [14]. The model is only based on three differential equations. Further analysis by several authors has also been done; see [5, 6, 1618].

Immunotherapy models and their predictions have been extensively studied in [9, 14, 19]. In [14], Kirschner et al. explored the role of cytokine in the disease dynamics and studied the long-term tumor recurrence and short-term tumor oscillations. However, in [19] Kuznetsov and Knott presented a mathematical model for the growth and suppression of the tumor. They showed that the model can describe the regrowth of a dormant tumor by two distinct mechanisms. One explanation for the tumor regrowth is based on a single clone model, while the other is based on a two-clone model. They fitted their ODE models to the data and obtained several curves for the tumor regrowth. They compared their predicted results with clinical and experimental observations, where both results confirm that intensive limited-term immunotherapy does not provide complete tumor elimination. The simulations show that medium-term control of cancer is exhibited when long-life immune memory cells are activated, but long-term control results from reducing the cancer growth rate.

In this paper, we investigate mathematical models for the dynamics between tumor cells, immune-effector cells, and the cytokine interleukin-2 (IL-2). It is worth stressing that we operate at a supermacroscopic scale, namely, by ordinary differential equations. However, the link to lower cellular scale is represented by the delay. The delay differential equations have long been used in modeling cancer phenomena [2026]. It should be noted that the heterogeneity, mutations, and link with the lower molecular scale are neglected. These topics are documented in [17, 27, 28].

The organization of this paper is as follows: in Section 2, we provide different models, using ODEs and DDEs, with interaction functions in the Lotka-Volterra form to describe the response of the effector cells to the growth of tumor cells. In Section 3, we study the local stability of the steady states for tumor-free and endemic persistence. Bifurcation analysis for a three-equations model and finding regions of existence of the equilibria are discussed in Section 4. In Section 5, we discuss the conditions that ensure tumor-clearance possibilities and conclude in Section 6.

2. The Model

The model of Kuznetsov et al. [7] describes the response of the effector cells (ECs) to the growth of tumor cells (TCs). In this model, it has been taken into account the penetration of TCs by ECs, which simultaneously causes the inactivation of ECs. It is assumed that interactions between ECs and TCs are in vitro such that 𝐸, 𝑇, 𝐶, 𝐸, and 𝑇 denote the local concentrations of ECs, TCs, EC-TC conjugates, inactivated effector cells, and “lethally hit” TCs, respectively. The rate of binding of ECs to TCs and the rate of separation of ECs from TCs without damaging them are denoted by 𝑘1 and 𝑘1, respectively. The rate at which EC-TC integrations program for lysis is denoted by 𝑘2, while the rate at which EC-TC interaction inactivate ECs is denoted by 𝑘3. The model takes the form 𝑑𝐸𝑑𝑡=𝑠+𝐹(𝐶,𝑇)𝑑1𝐸𝑘1𝑘𝐸𝑇+1+𝑘2𝐶,𝑑𝑇𝑑𝑡=𝑎𝑇(1𝑏𝑇)𝑘1𝑘𝐸𝑇+1+𝑘3𝐶,𝑑𝐶𝑑𝑡=𝑘1𝑘𝐸𝑇1+𝑘2+𝑘3𝐶,𝑑𝐸𝑑t=𝑎𝑘3𝐶𝑑3𝐸,𝑑𝑇𝑑𝑡=𝑘2𝐶𝑑3𝑇.(2.1) Here, the parameter 𝑠 represents the normal rate (not increased by the presence of the tumor) of the flow of adult ECs into the tumor site and 𝐹(𝐶,𝑇) describes the accumulation of ECs in the tumor site, while 𝑑1, 𝑑2, and 𝑑3 are the coefficients of the processes of destruction and migration of 𝐸, 𝐸, and 𝑇, respectively. The maximal growth of tumor is represented by the coefficient 𝑎, and 𝑏 is the environment capacity. It was suggested in [7] that the function 𝐹 takes the form 𝐹(𝐶,𝑇)=𝐹(𝐸,𝑇)=𝑝𝐸𝑇,𝑟+𝑇(2.2) where 𝑝 and 𝑟 are positive constants. This term is the Michaelis-Menten form to indicate the saturated effects of the immune response.

The idea in this paper is to simplify the above model and reduce it into a two- or three-equation model to describe the interactions of three types of cell populations: the activated immune-system cells, 𝐸(𝑡) (or effector cells such as cytotoxic T-cells, macrophages, and natural killer cells that are cytotoxic to the tumor cells); the tumor cells, 𝑇(𝑡); the concentration of IL-2 in the single tumor-site compartment, 𝐼𝐿(𝑡). The above model can then be governed by the following three equations (see [14]):𝑑𝐸=𝑑𝑡𝑐𝑇𝜇1𝐸+𝜃1𝐸𝐼𝐿+𝑠1,𝑑𝑇=𝑑𝑡𝑟2𝑇1𝑏𝑇𝛼𝐸𝑇,𝑑𝐼𝐿=𝑑𝑡𝜃2𝐸𝑇𝜇2𝐼𝐿+𝑠2,(2.3) with initial conditions 𝐸(0)=𝐸0, 𝑇(0)=𝑇0, 𝐼𝐿(0)=𝐼𝐿0, where 𝑐 is the antigenicity rate of the tumor, 𝑠1 is the external source of the effector cells, with rate of death 𝜇1, whereas the parameter 𝑟2 incorporates both multiplication and death of tumor cells. The maximal carrying capacity of the biological environment for tumor cell is 𝑏1, 𝜃1 is considered as the cooperation rate of effector cells with Interleukin-2 parameter, 𝛼 is the rate of tumor cells, and 𝜃2 is the competition rate between the effector cells and the tumor cells. External input of IL-2 into the system is 𝑠2, and the rate loss parameter of effector cells is 𝜇2.

2.1. Nondimensionalization

System (2.3) is an example of stiff  (One definition of the stiffness is that the global accuracy of the numerical solution is determined by stability rather than local error and implicit methods are more appropriate for it.) model, in the sense that it has properties that make it slow and expensive to solve using explicit numerical methods. Stiffness often appear due to the differences in speed between the fastest and slowest components of the solutions, and stability constraints. The efficient use of reliable numerical methods, that is based in general on implicit formulae, for dealing with stiff problems involves a degree of sophistication not necessarily available to nonspecialists [29]. In addition, the state variables of these types of models are very sensitive to small perturbations (or changes) in the parameters occuring in the model. Consequently, the parameter estimates are also sensitive to the noisy data and observations. To ease the analysis and stability of the steady states with meaningful parameters and less sensitive (or rubus) model, we nondimensionalize the bilinear model (2.3), by taking the following rescaling: 𝐸𝑥=𝐸0𝑇,𝑦=𝑇0𝐼,𝑧=𝐿𝐼𝐿0,𝜃1=𝜃1𝐼𝐿0𝑡𝑠,𝜃2=𝜃2𝐸0𝑇0𝑡𝑠𝐼𝐿0,𝜇1=𝜇1𝑡𝑠,𝜇2=𝜇2𝑡𝑠,𝑏=𝑏𝑇0,𝑐=𝑐𝑇0𝑡𝑠𝐸0,𝛼=𝛼𝐸0𝑡𝑠,𝜏=𝑡𝑠𝑡,𝑟2=𝑟2𝑡𝑠,𝑠1=𝑠1𝑡𝑠𝐸0,𝑠2=𝑠2𝑡𝑠𝐼𝐿0.(2.4) Therefore, after the above substitution into (2.3) and replacing 𝜏 by 𝑡, the model becomes𝑑𝑥𝑑𝑡=𝑐𝑦𝜇1𝑥+𝜃1𝑥𝑧+𝑠1,𝑑𝑦𝑑𝑡=𝑟2𝑦(1𝑏𝑦)𝛼𝑥𝑦,𝑑𝑧𝑑𝑡=𝜃1𝑥𝑦𝜇2𝑧+𝑠2,(2.5) with initial conditions 𝑥(0)=𝑥0, 𝑦(0)=𝑦0, and 𝑧(0)=𝑧0. Here 𝑥(𝑡), 𝑦(𝑡), and 𝑧(𝑡) denote the dimensionless density of ECs, TCs, and LI-2, respectively. In model (2.5), there are four possible cases of treatments, according the values of 𝑠1 and 𝑠2: (i) notreatment case (𝑠1=𝑠2=0), (ii) adoptive cellular immunotherapy case (𝑠1>0, 𝑠2=0), (iii) interleukin-2 case (𝑠1=0, 𝑠2>0), (iv) and immunotherapy with both adoptive cellular immunotherapy (ACI) and IL-2 (𝑠1>0, 𝑠2>0).

Yafia [10] considered system (2.5) in the absence of immunotherapy with IL-2,𝑑𝑥𝑑𝑡=𝜔𝑥𝑦𝜇𝑥+𝑠,𝑑𝑦𝑑𝑡=𝑟𝑦(1𝑏𝑦)𝑥𝑦,(2.6) where 𝜔 is immune response to the appearance of the TCs, 𝑠 has the same meaning of 𝑠1, 𝑟 has the meaning of 𝑟2, and 𝜇 has the meaning of 𝜇1 in the above model. If we consider a time delay 𝜏>0 in (2.6) due to the time-lag in the interaction between ECs and TCs, the model takes the form𝑑𝑥𝑑𝑡=𝜔𝑥(𝑡)𝑦(𝑡𝜏)𝜇𝑥(𝑡)+𝑠,𝑑𝑦𝑑𝑡=𝑟𝑦(𝑡𝜏)(1𝑏𝑦(𝑡))𝑥(𝑡)𝑦(𝑡).(2.7) In this model, we only consider the time delay in the dependent variable 𝑦 (representing tumor) of the nonlinear term. Of course other models assume time delays in both variables 𝑥 and 𝑦 [6]. Further models that consider time delays when modeling tumor growth are discussed in [3032]. To solve model (2.7), we should provide an initial function with initial function 𝑦(𝑡)=𝜓(𝑡), 𝑡[𝜏,0] instead of the initial value 𝑦(0) at 𝑡=0 (see [3]). It has been shown that model (2.7) has visible and bounded solution (see [6, 11]). When the time delay is included in the simplified model (2.6), the state of returning tumor cells can be observed, as DDE models have richer dynamics than do ODE models; see the graphs displayed in Figures 1, 2, and 3. We next study the stability of the steady states of the above models, according the values of the parameters given in Table 1.

tab1
Table 1: The nondimensionalization parameters of bilinear model (2.5).
fig1
Figure 1: Solution of the DDEs (2.7) when 𝜔=0.01184, 𝜇=0.3747, 𝑠=0.1181, 𝑟=1.636, 𝑏=0.002, and 𝜏=0.8. This shows an unstable endemic equilibrium.
fig2
Figure 2: Solution of the DDEs (2.7) when 𝜔=0.01184, 𝜇=0.3747, 𝑠=0.2181, 𝑟=1.636, 𝑏=0.002 and 𝜏=0.8. This shows a stable endemic equilibrium.
fig3
Figure 3: Solution of the DDEs (2.7) when 𝜔=0.04184, 𝜇=0.03747 (a), and 𝜇=0.3747 (b), 𝑠=0.2181, 𝑟=1.636, 𝑏=0.002, and 𝜏=0.8. The tumor-free equilibrium is asymptotically stable in the left banner and unstable in the right banner.

3. Steady States and Stability

The solutions of practical interest should have nonnegative population 𝑥, 𝑦, and 𝑧. However, it is hard to find a closed analytical solution for the above nonlinear models, instated we can study their qualitative behavior by studying the stability of the steady states. We then assume that the parameters occuring in the models are also nonnegative.

3.1. Tumor-Free Equilibrium and Its Stability

To ease the analysis, we start with the 2-population model (2.6). The steady states of the reduced model (2.6) are the intersection of the null-clines 𝑑𝑥/𝑑𝑡=0 and 𝑑𝑦/𝑑𝑡=0. If 𝑦=0, the free tumor equilibrium is at (𝑥,𝑦)=(𝑠/𝜇,0). This steady state always exists, since 𝑠/𝜇>0. It is clear that the tumor-free equilibrium 𝐸0=(𝑥,𝑦)=(𝑠/𝜇,0) of the model (2.6) is asymptotically stable if 𝑟𝜇<𝑠 and unstable if 𝑟𝜇>𝑠. Whoever, when we consider the DDEs model, the characteristic equation of the linearized model of (2.7) at 𝐸0=(𝑠/𝜇,0) takes the form(𝜆+𝑟)𝜆𝑒𝜆𝜏+𝑠𝜇=0.(3.1) When 𝜏=0, it is clear that 𝐸0 is asymptotically stable when 𝑟𝜇<𝑠 and unstable otherwise. However, if 𝜏>0, (3.1) has a negative real root 𝜆=𝑟 and roots of𝜆𝑒𝜆𝜏+𝑠𝜇=0.(3.2) Puting 𝜆=𝜉𝑖 in (3.2) and separating real and imaginary parts yields 𝜉2=𝑟2𝑠𝜇2.(3.3) Therefore, when |𝑟𝜇|<|𝑠| there are no positive real root 𝜉. This shows that all the roots of (3.1) have negative real parts and 𝐸0 is asymptotically stable.

In case of the three-equation model (2.5), and at the equilibrium points, we have 0=𝑐𝑦𝜇1𝑥+𝜃1𝑥𝑧+𝑠1,0=𝑟2𝑦(1𝑏𝑦)𝛼𝑥𝑦,0=𝜃2𝑥𝑦𝜇2𝑧+𝑠2.(3.4) Putting 𝑦=0 yields the tumor-free equilibrium, namely,𝐸0=𝑠1𝜇2𝜇1𝜇2𝜃1𝑠2𝑠,0,2𝜇2.(3.5) It is clear that the infection-free equilibrium 𝐸0 exists if and only if 𝑠2<𝜇1𝜇2/𝜃1. Therefore, we restrict our analysis to the case where 𝑠2<𝜇1𝜇2/𝜃1. To study its stability, we consider the corresponding Jacobian matrix 𝐽𝐸0=𝜇1+𝜃1𝑠2/𝜇2𝑐𝜃1𝑠1𝜇2/𝜇1𝜇2𝜃1𝑠20𝑟2𝛼𝑠1𝜇2/𝜇1𝜇2𝜃1𝑠200𝜃2𝑠1𝜇2/𝜇1𝜇2𝜃1𝑠2𝜇2.(3.6) It has the eigenvalues 𝜇1+𝜃1𝑠2/𝜇2, 𝑟2𝛼𝑠1𝜇2/(𝜇1𝜇2𝜃1𝑠2), and 𝜇2. Therefore, 𝐸0 is locally asymptotically stable if and only if 𝑟2<𝛼𝑠1𝜇2/(𝜇1𝜇2𝜃1𝑠2), while otherwise it is an unstable saddlepoint.

3.2. Endemic Equilibrium and Its Stability

Consider again the two-equation model (2.6). If 𝑦0, the steady states are obtained by solving 𝜔𝑟𝑏𝑦2𝑟(𝜔+𝜇𝑏)𝑦+𝜇𝑟𝑠=0. In this case, we have two endemic equilibria 𝑃1=(𝑥1,𝑦1) and 𝑃2=(𝑥2,𝑦2), where 𝑥1=𝑟(𝑏𝜇𝜔)Δ2𝜔,𝑦1=𝑟(𝑏𝜇+𝜔)+Δ,𝑥2𝑟𝑏𝜔2=𝑟(𝑏𝜇𝜔)+Δ2𝜔,𝑦2=𝑟(𝑏𝜇+𝜔)Δ,2𝑟𝑏𝜔(3.7) with Δ=𝑟2(𝑏𝜇𝜔)2+4𝜔𝑟𝑏𝜇>0. The Jacobian matrix of the system (2.6) at the endemic equilibrium 𝑃1 is 𝐽endemic=𝜔𝑦1𝜇𝜔𝑥1𝑦1𝑟2𝑏𝑟𝑦1𝑥1.(3.8)

Proposition 3.1. If the endemic equilibrium 𝑃1 exists and has nonnegative coordinates, then 𝑡𝑟(𝐽endemic)>0 and 𝑃1 is unstable.

Proof. Since 𝐽𝑡𝑟endemic=𝜔2𝜔(𝑟𝑏+𝑏𝜇)𝑟𝑏2𝜇+2𝑏𝜔𝜔𝑟𝑏2𝑟𝑏𝜔𝑟2(𝑏𝜇+𝜔)24𝑟𝑏𝜔(𝑟𝜇𝑠),(3.9) then inequality 𝑡𝑟(𝐽endemic)>0 is true if 𝑟𝜔2𝜔(𝑟𝑏+𝑏𝜇)𝑟𝑏2𝜇>(𝑟𝑏𝜔)𝑟2(𝑏𝜇+𝜔)24𝑟𝑏𝜔(𝑟𝜇𝑠).(3.10) Therefore, when 𝑟𝜇<𝑠 and 𝜔<𝑏𝜇, we have 𝜔2𝜔𝑏(𝑟+𝜇)𝑟𝜇𝑏2>0 and hence both sides of the inequality are positive. Therefore, if the point 𝑃1 exists and has nonnegative coordinates, then 𝑡𝑟(𝐽endemic)>0 and the point 𝑃1 is unstable whenever 𝜔<𝑏𝜇 and 𝑟𝜇<𝑠.

Similarly, it is easy to prove the following proportion.

Proposition 3.2. If the point 𝑃2 exists and has nonnegative coordinates, then it is asymptotically stable.

We extend the above analysis to the case of the three-equation model (2.5). The tumor-persistent solutions are obtained by putting 𝑦0 and omitting 𝑥 and 𝑧 in (3.4) to get a scalar equation in the variable 𝑦 which reads 𝐹𝑟2,𝑦=𝐶3𝑦3+𝐶2𝑦2+𝐶1𝑦+𝐶0=0,(3.11) where 𝐶3=𝜃1𝜃2𝑟22𝑏2𝐶>0,2=2𝜃1𝜃2𝑏𝑟22𝐶<0,1=𝑟22𝜃1𝜃2+𝑟2𝑏𝜇1𝜇2𝑠2𝜃1𝛼+𝑐𝜇2𝛼2,𝐶0𝑟=𝛼2𝜃1𝑠2𝜇1𝜇2+𝜇2𝑠1𝛼.(3.12) Since 𝑠2<𝜇1𝜇2/𝜃1, then the coefficient 𝐶1 is always positive, while the coefficient 𝐶0 can take positive and negative values depending on the values of the model parameters. Also, (3.11) is welldefined for all 𝑦[0,1/𝑏]. Its left-hand side is a polynomial of degree three, and its zeros are not easy to be obtained in a closed-form. However, some conditions in the parameters occur in the model to ensure the existence of its solutions could be deduced. Equation (3.11) can also be seen as a bifurcation equation in 𝑟2 and 𝑦, where we keep all other parameters fixed. Once a solution 𝑦>0 of this equation has been obtained, we could find positive 𝑥 and 𝑧 from the other equations in (3.4). Therefore, there is a one-to-one correspondence between the solutions of (3.11) and the endemic stationary solutions. For more insights, we next study the bifurcation analysis.

4. Bifurcation Analysis of Model (2.5)

The bifurcation analysis gives a deeper analysis about the model. It answers the query that “how does the behavior of the solutions change as parameters change.” We restrict ourselves to only study the bifurcation analysis of ODEs models rather than DDEs models.

4.1. Bifurcation Points for the Parameter 𝑟2

In this subsection, we then analyze the bifurcation so that the tumor growth rate 𝑟2 acts as a bifurcation parameter. Therefore, to find the bifurcation point(s), we put 𝑦=0 in (3.11) to get 𝑟2=𝛼𝜇2𝑠1/(𝜇1𝜇2𝜃1𝑠2)=𝑟2. Hence, there is only one transcritical bifurcation point at𝑟2,𝑦=𝛼𝜇2𝑠1𝜇1𝜇2𝜃1𝑠2,0.(4.1) Now, we compute the direction of bifurcation at (𝑟2,0) so that 𝑑𝑦𝑑𝑟2||||(𝑟2,0)𝐹=𝑟2𝐹𝑦||||(𝑟2,0),(4.2) where 𝐹𝑟2||(𝑟2,0)𝜇=𝛼1𝜇2𝜃1𝑠2<0,𝐹𝑦||(𝑟2,0)=𝑐+𝑏𝑠1+𝜃1𝜃2𝜇2𝑠21𝜇1𝜇2𝜃1𝑠22𝜇2𝛼2>0.(4.3) Hence, the bifurcation at the point (𝑟2,0) is forward, irrespective of the values of the model-parameters. We notice that the model we consider here has only one bifurcation point (𝑟2,0) at which the bifurcation is forward.

4.2. Bifurcation Diagrams for the Parameter

The parameter 𝛼 is very important in the model that plays an effective role to define cancer behavior. We investigate numerically, in this subsection, the bifurcation of the model for the parameter 𝛼. We consider the four cases: no treatment case (𝑠1=𝑠2=0), adoptive cellular immunotherapy case (𝑠1>0, 𝑠2=0), interleukin-2 case (𝑠1=0, 𝑠2>0), and immunotherapy with both ACI and IL-2 case (𝑠1>0, 𝑠2>0).

If we solve 𝐹(𝑦,𝛼)=𝐹(𝛼,𝑦)=0 in 𝛼, we have 𝛼+=𝑟2𝑏(𝑦1/𝑏)2𝜇2𝑐𝑦+𝑠1𝜃1𝑠2𝜇1𝜇2+𝜇1𝜇2𝑠2𝜃124𝜃1𝜃2𝜇2𝑦𝑐𝑦+𝑠1,𝛼=𝑟2𝑏(𝑦1/𝑏)2𝜇2𝑐𝑦+𝑠1𝜃1𝑠2𝜇1𝜇2𝜇1𝜇2𝑠2𝜃124𝜃1𝜃2𝜇2𝑦𝑐𝑦+𝑠1.(4.4) To plot (𝛼,𝑦) in the interval 0<𝑦<1/𝑏, under the conditions that 𝜇1𝜇2𝑠2𝜃124𝜃1𝜃2𝜇2𝑦𝑐𝑦+𝑠1>0,4𝜃1𝜃2𝑐𝜇2𝑦24𝜃1𝜃2𝜇2𝑠1𝜇𝑦+1𝜇2𝑠2𝜃12>0,(4.5) we have 𝑦2+𝑠1𝑐𝜇𝑦1𝜇2𝑠2𝜃124𝜃1𝜃2𝑐𝜇2<0.(4.6) Then 12𝑠1𝑐𝑠1𝑐2+𝜇1𝜇2𝑠2𝜃12𝜃1𝜃2𝑐𝜇21<𝑦<2𝑠1𝑐+𝑠1𝑐2+𝜇1𝜇2𝑠2𝜃12𝜃1𝜃2𝑐𝜇2.(4.7) For 0<𝑦<𝑦+, where 𝑦+=12𝑠1𝑐+𝑠1𝑐2+𝜇1𝜇2𝑠2𝜃12𝜃1𝜃2𝑐𝜇2,(4.8) we have two cases 𝑦+<1/𝑏 or 𝑦+>1/𝑏. The graphs in Figure 4, which are obtained numerically, display the bifurcation diagrams for different cases, where (i) 𝑠1=𝑠2=0, (ii) 𝑠1=10, 𝑠2=0, (iii) 𝑠1=0, 𝑠2=40, and (vi) 𝑠1=10, 𝑠2=40.

fig4
Figure 4: shows the bifurcation diagrams for the bilinear model (2.5) for the parameter 𝛼: [—] represents the stable equilibrium, [- - -] represents the unstable equilibrium, [] is the stable limit cycles, while [∘] is the saddle node bifurcation. [•] is the transcritical bifurcation and [□] the supercritical Hopf bifurcation. The values of parameters are given in Table 1.

Given the threshold point (𝛼,𝑦)=(𝛼,0), the tumor clearance condition is 𝛼>𝛼, where 𝛼=𝜇1𝜇2𝜃1𝑠2𝑟2𝜇2𝑠1.(4.9) Therefore, when 𝑠1>0, then 𝛼>0 and 𝛼1/𝑠1. Thus, we can arrive to tumor clearance quickly when the value of 𝑠1 increases. We notice from Figure 4 that the locations of saddle node bifurcation points 𝐀 and 𝐁 bridge the one-positive equilibrium to the three-positive equilibria. The supercritical Hopf bifurcation point 𝐂 joints between existence of stable limit cycles and nonexistence of limit cycles. The transcritical bifurcation point 𝐃 at (𝛼,𝑦)=(𝛼,0) bridges the one-positive equilibrium and no positive equilibria.

4.3. Regions of Existence of the Equilibria

In addition to the tumor-free equilibrium, (3.11) may have one to three persistent-tumor equilibria, depending on the values of the model-parameters. However, before we proceed we provide the following proposition, which is helpful in the analysis.

Proposition 4.1. Equation (3.11) does not have two persistent-tumor equilibria if 𝑟2<𝑟2, where 𝑟2 is given in (4.1).

Proof. Since the bifurcation direction at the point (𝑟2,0) is always forward, then (3.11) has two positive roots, for 𝑟2<𝑟2, if and only if 𝐹(𝑟2,𝑦)=𝑦(𝑐3𝑦2+𝐶2𝑦+𝐶1)=0 has two positive zeros, where 𝐶2=2𝑏𝐶3=2𝑏𝜃1𝜃2𝜇2𝑠1𝛼𝜇1𝜇2𝜃1𝑠22,𝐶1=𝜃1𝜃2𝜇2𝑠1𝛼𝜇1𝜇2𝜃1𝑠22+𝜇2𝛼2𝑏𝑠1+𝑐𝜇2.(4.10) However, 𝐶224𝐶1𝐶3=4𝐶3𝜇2𝛼2(𝑏𝑠1+𝑐𝜇2)<0. Therefore, the proof is complete.

Now, to find the conditions, on the model parameters, being required for the existence of the persistent equilibria, we make the use of both Descant's rule of signs and the Sturm sequence [33]. In (3.11), it is clear that the coefficients 𝐶3, 𝐶2, and 𝐶1 have fixed signs, while the coefficient 𝐶0 can take positive or negative values. Hence, the number of feasible tumor-persistent equilibria (on the interval 𝑦[0,1/𝑏]) depends on the difference between the number of sign changes at 𝑦=0 and at 𝑦=1/𝑏 in the Sturm sequence {𝑃0(𝑦),𝑃1(𝑦),𝑃2(𝑦),𝑃3(𝑦)}, such that 𝑃0(𝑦)=𝐹(𝑦)=𝐶3𝑦3+𝐶2𝑦2+𝐶1𝑦+𝐶0,𝑃1(𝑦)=𝐹(𝑦)=3𝐶3𝑦2+2𝐶2𝑦+𝐶1,𝑃22(𝑦)=91𝑅𝑦+9𝑆,𝑃3𝐶(𝑦)=𝑇,where𝑅=223𝐶1𝐶3/𝐶3,𝐶𝑆=2𝐶19𝐶0𝐶3/𝐶3𝑅,𝑇=33𝐶3𝑆𝐶2𝑅23𝐶32.(4.11) Hence the number of sign changes depends on the sign of the coefficient 𝐶0 and the remainders 𝑅, 𝑆, and 𝑇. Table 2 shows the conditions required for the existence of persistent-tumor equilibria as well as their numbers, where we take into account that 𝐶0>0 implies 𝑆<0.

tab2
Table 2: The number of positive steady states (SS) is determined by the signs of the coefficients (3.11) and the signs of the quantities 𝑅, 𝑆, and 𝑇 from the Sturm sequence. Blank entries correspond to coefficients which may take positive, negative, or zero values.

We may note that the tumor-persistent equilibria do not exist for 𝐶0>0, while one or three equilibria exist depending on the other relevant quantities. The interest is to find the area where the tumor-free equilibrium is a global attractor. Based on Proposition 4.1 and Table 2, this area is determined by 𝐶0>0 that is equivalent to𝑟2<𝜇2𝑠1𝛼𝜇1𝜇2𝜃1𝑠2=𝑟2,(4.12) where 𝑟2 is the critical growth rate of the tumor cell population, separating between nonexistence and existence of positive endemic equilibria.

If we consider the general case of immunotherapy with both ACI and IL-2 treatments, then according to the conditions given in Table 2 and data displayed in Figure 4, then Figure 5 displays six stability regions in terms of the two parameters 𝛼, 𝑟2 (according to the number of positive equilibria and the limit cycles). However, Figure 6 shows phase spaces for different equilibria where there are stable steady states, unstable steady states, stable manifold, unstable manifold and initial conditions. For example, Figure 6(f) shows the trajectories where there are three tumor-persistent equilibria for which two of them are locally stable and one, lying in between, for tumor-free equilibrium, is unstable (that correspond to region 3.b in Figure 5). Figure 6(e) shows the trajectories where there are three tumor-persistent equilibria for which one of them is locally stable, one is stable limit cycles, and one, lying in between, for tumor-free equilibrium is unstable (that correspond to region 3.a in Figure 5).

475720.fig.005
Figure 5: The two-dimensional transition structure as a function of 𝛼 (the tumor cell is predated by effector cells rate) and 𝑟2 (the maximal growth rate of the tumor cells population): sd is the saddle node bifurcation, ts is the transcritical bifurcation, and hf is the supercritical Hopf bifurcation. (Immunotherapy with both ACI and IL-2 Case; see Table 1 and 𝑠1=40, 𝑠2=10.).
fig6
Figure 6: The phase spaces at different equilibria where stable steady states [•], unstable steady states [○], stable manifold [- - -], unstable manifold [], and initial conditions [+], exist. (with immunotherapy with both ACI and IL-2, see Table 1 and 𝑠1=40, 𝑠2=10.).

5. Tumor-Clearance Possibilities

Let us introduce the following definition to facilitate the analysis.

Definition 5.1. The threshold parameter 0 (the minimum tumor-clearance parameter) is the parameter that has the property that if 0<1, then the endemic tumor does not exist, while if 0>1 the tumor persists (see [34]).

The parameter 0 can then be expressed in terms of the ratio between the tumor-growth rate and the critical tumor-growth rate separating between nonexistence and existence of endemic tumor. Now, for the three-equations model, the tumor-free equilibrium is the unique equilibrium if and only if 𝑟2<𝑟2, where it is also locally asymptotically stable. Therefore, the minimum tumor-clearance parameter is 0=𝑟2/𝑟2 and clearing the tumor requires the achievement of the inequality 0<1. It is equivalent to the following set of inequalities:𝑠1>𝜇1𝜇2𝜃1𝑠2𝜇2𝛼𝑟2=𝑠1,𝑠2<𝜇1𝜇2𝜃1=𝑠2.(5.1) Hence clearing the tumor depends mainly on the concentration of treatments: the external source of effector cells 𝑠1 and the treatment 𝑠2, which represents the external input of IL-2. If 𝑠2=0, then the tumor can be cleared by treatment with adoptive cellular immunotherapy alone, 𝑠1>(𝜇1/𝛼)𝑟2. However, for 𝑠1=0, then the inequality 𝑟2<𝑟2 cannot be held and therefore, it is impossible to treat cancer by IL-2 alone. However, a strategy based on using both adoptive immunotherapy and IL-2 with concentrations 𝑠2<𝑠2 and 𝑠1>𝑠1 could be used to clear the tumor; see Figure 7. We arrive to the following corollary.

475720.fig.007
Figure 7: The critical 𝑠1 as a function of 𝑟2 for several levels of 𝑠2. For values of 𝑠1 above the threshold 𝑠1, the tumor cells do not exist (with 𝛼=0.5556), see Table 1.

Corollary 5.2. In the tumor-clearance problem, we have the following three cases:(i)if 𝑠1=0, the tumor could never be cleared,(ii)if 𝑠2=0, the tumor could be cleared by adding an external source of effector cells with concentration slightly above 𝑠1=(𝜇1/𝛼)𝑟2,(iii)if 𝑠10 and 𝑠20, then the tumor could be cleared with concentrations 𝑠2<𝜇1𝜇2/𝜃1 and 𝑠1>((𝜇1𝜇2𝜃1𝑠2)/𝜇2𝛼)𝑟2.

6. Summary and Conclusion

In this paper, we introduced a family of differential models (ODEs and DDEs) to describe the dynamics of cancer. The ODEs models model cancer at supermacroscopic, in the sense that they describe the interaction between the tumor cells and the normal (immune) cells [1]. However, the DDEs models link it with the lower cellular scale. The qualitative and evolution of the models have been displayed with different values of the parameters 𝛼 (the rate of tumor cells predated by the effector cells) and 𝑟2 (the maximal growth rate of the tumor cells population). Although the underlying models are simple, they display very rich dynamics and give a good picture for the phenomena of real interaction of tumor growth and immunotherapy. The minimum tumor-clearance parameter 0 has been expressed in terms of the ratio between the tumor-growth rate and the critical tumor-growth rate. The cases at which the tumor can be cleared are summarized in Corollary 5.2. The obtained results can help to gain a better understanding of interaction mechanisms and make predictions, determine and evaluate control strategies, and convey more general insight to biologists.

The numerical simulations (have been obtained by semi-implicit RK methods [29]) demonstrate that the system with time delay exhibits richer complex dynamics, such as quasiperiodic and chaotic patterns, compared with models without memory or after-effect. The steady states of DDEs models are similar to the steady states of ODEs models. We shall extend this work to investigate the qualitative behavior and bifurcation analysis of more sophisticated models of DDEs in modeling tumor-immune interactions with immunotherapy and control functionals to maximize the effector cells and interleukin-2 concentration and to minimize the tumor cells.

Acknowledgments

The authors would like to thank the referees and Professor Pedro Serranho for their valuable comments on the paper. The first author should thank Emirates Foundation for partially funding this research.

References

  1. N. Bellomo, N. K. Li, and P. K. Maini, “On the foundations of cancer modelling: selected topics, speculations, and perspectives,” Mathematical Models & Methods in Applied Sciences, vol. 18, no. 4, pp. 593–646, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  2. N. Bellomo, A. Bellouquid, J. Nieto, and J. Soler, “Multiscale biological tissue models and flux-limited chemotaxis for multicellular growing systems,” Mathematical Models & Methods in Applied Sciences, vol. 20, no. 7, pp. 1179–1207, 2010. View at Publisher · View at Google Scholar
  3. F. Rihan, “Delay differential models in dynamic diseases,” in Proceedings of the International Conference on Bioinformatics and Computational Biology, vol. 2, pp. 73–79, Honolulu, Hawaii, USA, 2010.
  4. R. J. DeBoer, P. Hogeweg, F. Dullens, R. D. Weger, and W. DenOtter, “Macrophage T lymphocyte interactions in the antitumor immune response: a mathematical model,” Journal of Immunology, vol. 134, no. 4, pp. 2748–2758, 1985. View at Google Scholar
  5. C. DeLisi and A. Rescigno, “Immune surveillance and neoplasia. I. A minimal mathematical model,” Bulletin of Mathematical Biology, vol. 39, no. 2, pp. 201–221, 1977. View at Google Scholar · View at Zentralblatt MATH
  6. M. Gałach, “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 · View at Zentralblatt MATH
  7. 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 Google Scholar · View at Zentralblatt MATH · View at Scopus
  8. V. A. Kuznetsov, “Mathematical modeling of the development of dormant tumors and immune stimulation of their growth,” Cybernetics and Systems Analysis, vol. 23, no. 4, pp. 556–564, 1987. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  9. A. Matzavinos, M. A. J. Chaplain, and V. A. Kuznetsov, “Mathematical modelling of the spatio-temporal response of cytotoxic T-lymphocytes to a solid tumour,” Mathematical Medicine and Biology, vol. 21, no. 1, pp. 1–34, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  10. R. Yafia, “Hopf bifurcation analysis and numerical simulations in an ODE model of the immune system with positive immune response,” Nonlinear Analysis, vol. 8, no. 5, pp. 1359–1369, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  11. R. Yafia, “Hopf bifurcation in differential equations with delay for tumor-immune system competition model,” SIAM Journal on Applied Mathematics, vol. 67, no. 6, pp. 1693–1703, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  12. 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 · View at Google Scholar
  13. 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 · View at Google Scholar · View at Scopus
  14. D. Kirschner and J. C. Panetta, “Modeling immunotherapy of the tumor—immune interaction,” Journal of Mathematical Biology, vol. 37, no. 3, pp. 235–252, 1998. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus
  15. D. Liu, S. Ruan, and D. Zhu, “Bifurcation analysis in models of tumor and immune system interactions,” Discrete and Continuous Dynamical Systems, vol. 12, no. 1, pp. 151–168, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  16. J. C. Arciero, T. L. Jackson, and D. E. Kirschner, “A mathematical model of tumor-immune evasion and siRNA treatment,” Discrete and Continuous Dynamical Systems, vol. 4, no. 1, pp. 39–58, 2004. View at Google Scholar · View at Zentralblatt MATH
  17. N. Bellomo and M. Delitala, “From the mathematical kinetic, and stochastic game theory to modelling mutations, onset, progression and immune competition of cancer cells,” Physics of Life Reviews, vol. 5, no. 4, pp. 183–206, 2008. View at Publisher · View at Google Scholar · View at Scopus
  18. O. Isaeva and V. Osipov, “Modelling of anti-tumor immune response: immunocorrective effect of weak centimeter electromagnetic waves,” Journal of Mathematical Methods in Medicine, vol. 10, no. 3, pp. 185–201, 2009. View at Google Scholar
  19. V. A. Kuznetsov and G. D. Knott, “Modeling tumor regrowth and immunotherapy,” Mathematical and Computer Modelling, vol. 33, no. 12-13, pp. 1275–1287, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  20. M. Bodnar, U. Foryś, and J. Poleszczuk, “Analysis of biochemical reactions models with delays,” Journal of Mathematical Analysis and Applications, vol. 376, no. 1, pp. 74–83, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  21. U. Foryś, “Multi-dimensional Lotka-Volterra systems for carcinogenesis mutations,” Mathematical Methods in the Applied Sciences, vol. 32, no. 17, pp. 2287–2308, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  22. U. Foryś, M. Bodnar, and J. Poleszczuk, “Negativity of delayed induced oscillations in a simple linear DDE,” Applied Mathematics Letters, vol. 24, no. 6, pp. 982–986, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  23. J. Miekisz, J. Poleszczuk, M. Bodnar, and U. Forys, “Stochastic models of gene expression with delayed degradation,” Bulletin of Mathematical Biology, vol. 73, no. 9, pp. 2231–2247, 2011. View at Publisher · View at Google Scholar
  24. M. J. Piotrowska and U. Foryś, “Analysis of the Hopf bifurcation for the family of angiogenesis models,” Journal of Mathematical Analysis and Applications, vol. 382, no. 1, pp. 180–203, 2011. View at Publisher · View at Google Scholar
  25. M. Piotrowska and U. Forys, “The nature of hopf bifurcation for the gompertz model with delays,” Mathematical and Computer Modelling, vol. 54, pp. 9–10, 2011. View at Google Scholar
  26. J. Poleszczuk, M. Bodnar, and U. Foryś, “New approach to modeling of antiangiogenic treatment on the basis of Hahnfeldt et al. model,” Mathematical Biosciences and Engineering, vol. 8, no. 2, pp. 591–603, 2011. View at Publisher · View at Google Scholar
  27. E. Gabetta and E. Regazzini, “About the gene families size distribution in a recent model of genome evolution,” Mathematical Models & Methods in Applied Sciences, vol. 20, no. 6, pp. 1005–1020, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  28. J. Paulsson, “Models of stochastic gene expression,” Physics of Life Reviews, vol. 2, no. 2, pp. 157–175, 2005. View at Publisher · View at Google Scholar · View at Scopus
  29. F. A. Rihan, E. H. Doha, M. I. Hassan, and N. M. Kamel, “Numerical treatments for Volterra delay integro-differential equations,” Computational Methods in Applied Mathematics, vol. 9, no. 3, pp. 292–308, 2009. View at Google Scholar · View at Zentralblatt MATH
  30. M. Bodnar and U. Foryś, “Behaviour of solutions to Marchuk's model depending on a time delay,” International Journal of Applied Mathematics and Computer Science, vol. 10, no. 1, pp. 97–112, 2000. View at Google Scholar · View at Zentralblatt MATH
  31. M. Bodnar and U. Foryś, “Periodic dynamics in a model of immune system,” Applicationes Mathematicae, vol. 27, no. 1, pp. 113–126, 2000. View at Google Scholar · View at Zentralblatt MATH
  32. H. M. Byrne, “The effect of time delays on the dynamics of avascular tumor growth,” Mathematical Biosciences, vol. 144, no. 2, pp. 83–117, 1997. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  33. R. Beaumont and R. Pierce, The Algebraic Foundations of Mathematics, Addison-Wesely, Reading, Mass, USA, 1963.
  34. M. Safan, H. Heesterbeek, and K. Dietz, “The minimum effort required to eradicate infections in models with backward bifurcation,” Journal of Mathematical Biology, vol. 53, no. 4, pp. 703–718, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH