About this Journal Submit a Manuscript Table of Contents
Journal of Applied Mathematics
Volume 2012 (2012), Article ID 891095, 16 pages
http://dx.doi.org/10.1155/2012/891095
Research Article

Chaos in a Tumor Growth Model with Delayed Responses of the Immune System

Department of Applied Mathematics, Z. H. College of Engineering & Technology, A.M.U, Aligarh 202002, India

Received 3 November 2011; Revised 22 December 2011; Accepted 29 December 2011

Academic Editor: Mohamad Alwash

Copyright © 2012 M. Saleem and Tanuja Agrawal. 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

A simple prey-predator-type model for the growth of tumor with discrete time delay in the immune system is considered. It is assumed that the resting and hunting cells make the immune system. The present model modifies the model of El-Gohary (2008) in that it allows delay effects in the growth process of the hunting cells. Qualitative and numerical analyses for the stability of equilibriums of the model are presented. Length of the time delay that preserves stability is given. It is found that small delays guarantee stability at the equilibrium level (stable focus) but the delays greater than a critical value may produce periodic solutions through Hopf bifurcation and larger delays may even lead to chaotic attractors. Implications of these results are discussed.

1. Introduction

It is well known that cancer is one of the greatest killers in the world, and the control of tumor growth requires great attention. The development of a cancerous tumor is complex and involves interaction of many cell types. Main components of these cells are tumor cells (or abnormal cells also known as bad cells) and immune and healthy tissue cells (or normal cells also known as good cells).

A tumor is a dynamic nonlinear system, in which bad cells grow, spread, and eventually overwhelm good cells in the body. The form of the dynamic nonlinear system modeling the cancer and the class of the equations that describe such a system are related to the scaling problem. Indeed, there are three natural scales that are connected to different stages of the disease and have to be identified. The first is the subcellular (or molecular) scale, where one focuses on studying the alterations in the genetic expressions of the genes contained in the nucleus of a cell. As a result of this, some special signals, which are received by the receptors on the cell surface, are transmitted to the cell nucleus. The second is the cellular scale, which is an intermediate level between the molecular scale and the macroscopic scale to be described in the following. The third is the macroscopic scale, where one deals with heterogeneous tissues. In the heterogeneous tissues, some of the layers (e.g., the external proliferating layer, the intermediate layer, and the inner zone with necrotic cells) constituting the tumor may occur as islands. This leads to a tumor comprising of multiple regions of necrosis engulfed by tumor cells in a quiescent or proliferative state [1]. In case of macroscopic scale, the main focus is 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 one scale to another, one may 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 the tumor and immune system populations over time.

Many authors have used mathematical models to describe the interaction among the various components of tumor microenvironment, (see de Boer et al. [4], Goldstein et al. [5], De Pillis et al. [6], and Kronic et al. [7]). These papers mainly deal with immune response to tumor growth. In the last few years a great deal of human and economical resources is devoted to cancer research with a view to develop different control strategies and drug therapies with main emphasis on experimental aspects and immunology (see Aroesty et al. [8], Eisen [9], Knolle [10], Murray [11], Adam [12], Adam and Panetta [13], Owen and Sherrat [14], De Pillis and Radunskaya [15], Dingli et al. [16], and Menchón et al. [17]). There are many existing reviews of mathematical models of tumor growth and tumor immune system interactions such as Bellomo and Preziosi [18], Araujo and McElwain [19], Nagy [20], Byrne et al. [21], Castiglione and Piccoli [22], Martins et al. [23], Roose et al. [24], Chaplain [25], and Bellomo et al. [26]. Some of these reviews follow a historical approach (Araujo and McElwain [19]), while others focus on multiscale modeling or on particular aspects of tumor evolution (Bellomo and Preziosi [18], Martins et al. [23], and Bellomo et al. [26]). Recently, Bellomo et al. [27] study the competition between tumor and immune cells modeled by a nonlinear dynamical system, which identifies the evolution of the number of cells belonging to different interacting populations such as tumor and immune cells at different scales, namely molecular, cellular and macroscopic. Bellomo and Delitala [28] have applied the methods of the classical mathematical kinetic theory for active particles to study the immune competition with special attention to cancer phenomena. They mainly focus on modeling aspects of the early stage of cancer onset and competition with the immune system.

Several authors have used the concept of prey-predator-type interactions in tumor studies where in general the immune cells play the role of predator and the tumor cells that of prey (see Kuznetsov et al. [29], Kirschner and Panetta [30], Sarkar and Banerjee [31], and El-Gohary [32]). These are mainly ordinary differential equation models, which certainly provide a simpler framework to explore the interactions among tumor cells and the different types of immune and healthy tissue cells. Kuznetsov et al. [29] study nonlinear dynamics of immunogenic tumors with emphasis on parameter estimation and global bifurcation analysis. Immunotherapy of tumor-immune interaction has been studied by Kirschner and Panetta [30]. They indicated that the dynamics between tumor cells, immune cells, and IL-2 can explain both short-term oscillations in tumor size as well as long-term tumor relapse. Sarkar and Banerjee [31] discuss self-remission and tumor stability by taking stochastic approach.

The delay differential equations have long been used in modeling cancer phenomena [3339]. Byrne [40] considers the effect of time delay on the dynamics of avascular tumor growth by incorporating a time-delayed factor into the net proliferation rate of the cells. Burić et al. [41] consider the effects of time delay on the two-dimensional system, which represents the basic model of the immune response. They study variations of the stability of the fixed points due to time delay and the possibility for the occurrence of the chaotic solutions. Recently, Foryś and Kolev [42] propose and study the role of time delay in solid avascular tumor growth. They study a delay model in terms of a reaction-diffusion equation and mass conservation law. Two main processes are taken into account that is, proliferation and apoptosis. Gałach [43] studies a simplified version of the Kuznetsov-Taylor model, where immune reactions are described by a bilinear term with time delay. Yafia [44] analyzes an interaction between the proliferating and quiescent cells tumor with a single delay. He shows the occurrence of Hopf bifurcation as the delay crosses some critical value.

Recently, El-Gohary [32] studied a cancer self-remission and tumor system and provided optimal control strategies that made its unstable steady states asymptotically stable. In the present paper, we modify the model of El-Gohary [32] by introducing a constant time delay 𝑇 in the growth rate of the hunting cells of the immune system. This modification, while on one hand incorporates certain thresholds that may be helpful to control the tumor cell growth, on the other hand hints at the complex dynamics that a tumor may have. It may be mentioned here that by representing tumor growth with ordinary differential equations we indeed operate in the present study at the supermacroscopic scale, while the link with the lower cellular scale is represented by the delay. Of course, we do not consider heterogeneity, mutations, and link with the lower molecular scale in the present paper (for details one may see [27, 28, 45, 46]).

2. The Model and Equilibrium Solutions

El-Gohary [32] considered the following model for cancer self-remission and tumor growth:𝑑𝑀𝑀𝑑𝑡=𝑞+𝑟𝑀1𝑘1𝛼𝑀𝑁,𝑑𝑁𝑑𝑡=𝛽𝑁𝑍𝑑1𝑁,𝑑𝑍𝑍𝑑𝑡=𝑠𝑍1𝑘2𝛽𝑁𝑍𝑑2𝑍.(2.1) In (2.1), different variables and parameters have the following interpretations.

𝑀(𝑡),  𝑁(𝑡),𝑍(𝑡): densities of tumor cells, hunting predator cells, and resting cells at time 𝑡𝑞:conversion of normal cells to malignant cells,𝑟:growth rate of tumor cells,𝑘1:maximum carrying capacity of tumor cells,𝛼:rate of killing of tumor cells by hunting cells,𝛽:conversion rate of the resting cells to tumor cells,𝑑1:natural death rate of hunting cells,𝑠:growth rate of resting cells,𝑑2:natural death rate of resting cells.

Using nondimensional variables and parameters as𝜏=𝑞𝑘11𝑡,𝑥1=𝑘11𝑀,𝑥2=𝛼𝑘1𝑞1𝑁,𝑥3=𝑘21𝑎𝑍,1=𝑟𝑘1𝑞1,𝑎2=𝛽𝑘1𝑘2𝑞1,𝑎3=𝑑1𝑘1𝑞1,𝑎4=𝑠𝑘1𝑞1,𝑎5=𝛽𝛼1,𝑎6=𝑑2𝑘1𝑞1,(2.2) El-Gohary [32] obtained the following nondimensional form of model (2.1): 𝑑𝑥1𝑑𝜏=1+𝑎1𝑥11𝑥1𝑥1𝑥2,𝑑𝑥2𝑑𝜏=𝑎2𝑥2𝑥3𝑎3𝑥2,𝑑𝑥3𝑑𝜏=𝑎4𝑥31𝑥3𝑎5𝑥2𝑥3𝑎6𝑥3.(2.3)

We modify model (2.3) by assuming that there is a constant time delay 𝑇 since the time resting cells give a signal to hunting cells for activation and the mature hunting cells are ready to kill the tumor cells. More specifically, we incorporate this assumption by replacing the growth term 𝑎2𝑥2(𝜏)𝑥3(𝜏) in (2.3) by 𝑎2𝑥2(𝜏𝑇)𝑥3(𝜏𝑇). Thus our model takes the form𝑑𝑥1𝑑𝜏=1+𝑎1𝑥11𝑥1𝑥1𝑥2,𝑑𝑥2𝑑𝜏=𝑎2𝑥2(𝜏𝑇)𝑥3(𝑡𝑇)𝑎3𝑥2,𝑑𝑥3𝑑𝜏=𝑎4𝑥31𝑥3𝑎5𝑥2𝑥3𝑎6𝑥3.(2.4) It is obvious that model (2.4) would have same equilibrium solutions as model (2.3), whose equilibriums have been reported as𝐸1=121+41+𝑎1,𝐸,0,02=121+41+𝑎1𝑎,0,16𝑎4,𝐸3=𝑥1,𝑥2,𝑥3=12𝑎1𝑎1𝑥2+𝑎1𝑥22+4𝑎1,𝑎4𝑎5𝑎13𝑎2𝑎6𝑎5,𝑎3𝑎2,(2.5) in El-Gohary [32], under the biologically feasible conditions as𝑎4>𝑎6,𝑎3𝑎2+𝑎6𝑎4<1.(2.6) The first condition states that the ratio of the natural death rate of resting cell to its growth rate is less than one. The second condition implies that the ratio of the natural death rate of the hunting cell to its growth rate is also less than one.

3. Linear Stability Analysis

3.1. Stability without Delay (i.e., 𝑇=0)

The local stability results for equilibriums of model (2.4) for 𝑇=0 have been reported in El-Gohary [32] as shown in (Table 1).

tab1
Table 1
3.2. Stability with Delay (i.e., 𝑇0)

Assuming small deviations 𝑢𝑖 around the equilibrium 𝐸3=(𝑥1,𝑥2,𝑥3) such that 𝑢𝑖=𝑥𝑖𝑥𝑖, the linearized system of the model (2.4) becomes𝑑𝑢1=𝑎𝑑𝜏112𝑥1𝑥2𝑢1𝑥1𝑢2,𝑑𝑢2𝑑𝜏=𝑎3𝑢2+𝑎2𝑥3𝑢2(𝜏𝑇)+𝑎2𝑥2𝑢3(𝜏𝑇),𝑑𝑢3=𝑎𝑑𝜏5𝑎3𝑢2𝑎2𝑎4𝑎3𝑢3𝑎2.(3.1) In the case of a positive delay (𝑇>0), the characteristic equation for this system can be written as𝑃(𝜆)+𝑄(𝜆)𝑒𝜆𝑇=0,(3.2) where𝑃(𝜆)=𝜆3+𝐴1𝜆2+𝐴2𝜆+𝐴3,𝑄(𝜆)=𝐵1𝜆2+𝐵2𝜆+𝐵3(3.3) such that𝐴1=𝑎3+𝑎4𝑎3𝑎2𝐴+𝑘,2=𝑎4𝑎32𝑎2+𝑘𝑎3+𝑘𝑎4𝑎3𝑎2𝐴3=𝑘𝑎4𝑎23𝑎2,𝐵1=𝑎3,𝐵2𝑎=4𝑎23𝑎2+𝑎5𝑎3𝑥2𝑘𝑎3,𝐵3=𝑘𝑎4𝑎23𝑎2+𝑘𝑎5𝑎3𝑥2,(3.4) with𝑘=𝑎1𝑥22+4𝑎1.(3.5) Now substituting 𝜆=𝑖𝜔 (where 𝜔 is positive) in (3.2) and separating the real and imaginary parts, we obtain the following system of transcendental equations:𝐴1𝜔2𝐴3=𝐵3𝐵1𝜔2cos(𝜔𝑇)+𝐵2𝜔sin(𝜔𝑇),(3.6)𝜔3𝐴2𝜔=𝐵2𝐵𝜔cos(𝜔𝑇)3𝐵1𝜔2sin(𝜔𝑇).(3.7) Squaring and adding (3.6) and (3.7), we get𝐵3𝐵1𝜔22+𝐵22𝜔2=𝐴1𝜔2𝐴32+𝜔3𝐴2𝜔2,(3.8) which can be simplified to𝜔6+𝐶1𝜔4+𝐶2𝜔2+𝐶3=0,(3.9) where𝐶1=𝐴122𝐴2𝐵12,𝐶2=𝐴22𝐵222𝐴1𝐴3+2𝐵1𝐵3,𝐶3=𝐴32𝐵32=𝐴3+𝐵3𝐴3𝐵3.(3.10) Equation (3.9) can be written as a cubic𝜌3+𝐶1𝜌2+𝐶2𝜌+𝐶3=0(3.11) with 𝜌=𝜔2.

For parameter values such that 𝐶1 is positive, the simplest assumption that (3.11) will have a unique positive root is 𝐶3=𝐴32𝐵32<0. Since 𝐴3+𝐵3 is positive, it requires that 𝐴3𝐵3<0 for 𝐶3 to be negative. Hence, it can be said that there is a unique positive root say 𝜌0 of (3.11). Denoting 𝜔0=𝜌01/2, it follows that the characteristic equation (3.2) has a pair of purely imaginary roots of the form ±𝑖𝜔0. Eliminating sin(𝜔𝑇) from (3.6) and (3.7), we get𝐴cos(𝜔𝑇)=1𝜔2𝐴3𝐵3𝐵1𝜔2+𝜔3𝐴2𝜔𝐵2𝜔𝐵3𝐵1𝜔22+𝐵22𝜔2.(3.12) Then 𝑇𝑛 corresponding to 𝜔0 is given by𝑇𝑛=1𝜔0𝐴arccos1𝜔02𝐴3𝐵3𝐵1𝜔02+𝜔03𝐴2𝜔0𝐵2𝜔0𝐵3𝐵1𝜔022+𝐵22𝜔02+2𝑛𝜋𝜔0.(3.13) Since 𝐸3 is stable for 𝑇=0, it implies from Freedman and Rao [47] that 𝐸3 remains stable for 𝑇<𝑇0.

3.3. Estimation of the Length of Delay to Preserve Stability

Let us consider the linearized system (3.1). Taking the Laplace transform of this system, we get𝑠𝑎1+2𝑎1𝑥1+𝑥2𝑈1(𝑠)=𝑥1𝑈2(𝑠)+𝑢1(0),𝑠+𝑎3𝑈2(𝑠)=𝑎2𝑥3𝑒𝑠𝑇𝑈2(𝑠)+𝑎2𝑥3𝑒𝑠𝑇𝐾1(𝑠)+𝑎2𝑥2𝑒𝑠𝑇𝑈3(𝑠)+𝑎2𝑥2𝑒𝑠𝑇𝐾2(𝑠)+𝑢2(𝑎0),𝑠+4𝑎3𝑎2𝑈3𝑎(𝑠)=5𝑎3𝑎2𝑈2(𝑠)+𝑢3(0),(3.14) where𝑈𝑖𝑢(𝑠)=𝐿𝑖,𝐾(𝜏)1(𝑠)=0𝑇𝑒𝑠𝜏𝑢2𝐾(𝜏)𝑑𝜏,2(𝑠)=0𝑇𝑒𝑠𝜏𝑢3(𝜏)𝑑𝜏.(3.15) Following lines of Erbe et al. [48] and using the Nyquist criterion (Freedman and Rao [47]), it can be shown that the sufficient conditions for the local asymptotic stability of 𝐸3=(𝑥1,𝑥2,𝑥3) are given byIm𝐻𝑖𝜂0>0,(3.16)Re𝐻𝑖𝜂0=0,(3.17) where 𝐻(𝑠)=𝑠3+𝐴1𝑠2+𝐴2𝑠+𝐴3+𝑒𝜆𝑇(𝐵1𝑠2+𝐵2𝑠+𝐵3) and 𝜂0 is the smallest positive root of (3.17).

Inequality (3.16) and (3.17) can alternatively be written as𝐴2𝜂0𝜂30>𝐵2𝜂0𝜂cos0𝑇+𝐵3𝜂sin0𝑇𝐵1𝜂20𝜂sin0𝑇,(3.18)𝐴3𝐴1𝜂20=𝐵1𝜂20𝜂cos0𝑇𝐵3𝜂cos0𝑇𝐵2𝜂0𝜂sin0𝑇.(3.19) Now if (3.18) and (3.19) are satisfied simultaneously, they are sufficient conditions to guarantee stability. These are now used to get an estimate to the length of the time delay. The aim is to find an upper bound 𝜂+ to 𝜂0, independent of 𝑇, from (3.19) and then to estimate 𝑇 so that (3.18) holds true for all values of 𝜂such that 0𝜂𝜂+ and hence, in particular, at 𝜂=𝜂0.

Equation (3.19) is rewritten as𝐴1𝜂20=𝐴3𝐵1𝜂20𝜂cos0𝑇+𝐵3𝜂cos0𝑇+𝐵2𝜂0𝜂sin0𝑇.(3.20) Maximizing the right-hand side of (3.20)

subject to||𝜂sin0𝑇||||𝜂1,cos0𝑇||1,(3.21) we obtain||𝐴1||𝜂20||𝐴3||+||𝐵3||+||𝐵1||𝜂20+||𝐵2||𝜂0.(3.22) Hence if𝜂+=12||𝐴1||||𝐵1||||𝐵2||+𝐵22||𝐴+41||||𝐵1||||𝐴3||+||𝐵3||,(3.23) then clearly from (3.22) we have 𝜂0𝜂+.

From (3.18), we obtain𝜂20<𝐴2+𝐵2𝜂cos0𝑇+𝐵1𝜂0𝜂sin0𝑇𝐵3𝜂sin0𝑇𝜂0.(3.24) Since 𝐸3=(𝑥1,𝑥2,𝑥3) is locally asymptotically stable for 𝑇=0, the inequality (3.24) will continue to hold for sufficiently small 𝑇>0. Using (3.20), (3.24) can be rearranged as𝐵3𝐵1𝜂20𝐴1𝐵2𝜂cos0𝑇+𝐵12𝐴1𝐵1𝜂0+𝐴1𝐵3𝜂0𝜂×sin0𝑇<𝐴1𝐴2𝐴3𝐵3+𝐵1𝜂20+𝐴1𝐵2.(3.25) Using the bound𝐵3𝐵1𝜂20𝐴1𝐵2𝜂cos0𝑇𝐵1=3𝐵1𝜂20𝐴1𝐵22sin2𝜂0𝑇212||𝐵1𝜂2++𝐴1𝐵2𝐵3||𝜂2+𝑇2,𝐵2𝐴1𝐵1𝜂0+𝐴1𝐵3𝜂0𝜂sin0𝑇||𝐵2𝐴1𝐵1||𝜂2++||𝐴1||||𝐵3||𝑇,(3.26) we obtain from (3.24)𝐿1𝑇2+𝐿2𝑇<𝐿3,(3.27) where𝐿1=12||𝐵1𝜂2++𝐴1𝐵2𝐵3||𝜂2+,𝐿2=||𝐵2𝐴1𝐵1||𝜂2++||𝐴1||||𝐵3||,𝐿3=𝐴1𝐴2𝐴3𝐵3+𝐵1𝜂2++𝐴1𝐵2.(3.28) Hence, if𝑇+=12𝐿1𝐿2+𝐿22+4𝐿1𝐿3,(3.29) then for 0𝑇𝑇+ the Nyquist criterion holds true and 𝑇+ estimates the maximum length of the delay preserving the stability.

4. Numerical Simulation

The purpose of this section is to illustrate dynamics of the model (2.4) numerically with variation in the delayed responses of the immune system and relate it with the stability results of the model equilibrium mentioned in Section 3. For this purpose, we consider the following set of parameters 𝑎1=2.5,𝑎2=4.5,𝑎3=0.6,𝑎4=3.5,𝑎5=2, and 𝑎6=0.1, which satisfies the biologically feasible conditions (2.6). It may be mentioned that the same set of parameter values has been used by El-Gohary [32] in numerical calculations. For this set, the positive equilibrium of the model (2.4) is 𝐸3=(0.8720,1.4667,0.1333). It can be seen that for these parameter values the coefficients of the cubic (3.11) are such that 𝐶1>0and 𝐶3<0. It guarantees that the cubic (3.11) has a unique positive root 𝜌0=1.34838, which then provides the unique positive root of (3.9) as 𝜔0=𝜌0=1.1612. Using this value of 𝜔0, it turns out from (3.13) that 𝑇0=0.3427, and the stability result of Section 3 yields that the positive equilibrium 𝐸3 of model (2.4) is stable for 𝑇 such that 0𝑇<0.3427. For 𝑇=0.34, Figure 1 illustrates the approach of the trajectory of the model (2.4) to the equilibrium 𝐸3.

891095.fig.001
Figure 1: Three-dimensional phase portrait depicting stable dynamics of the model (2.4) for 𝑇=0.34. The trajectory spirals and approaches to 𝐸3 starting from the initial point 𝑥1(0)=1,𝑥2(0)=1,𝑥3(0)=1.

Indeed the existence of 𝜔0 as unique root of (3.9) implies that there exists a pair of pure imaginary eigen values 𝜆that satisfies the characteristic equation (3.2) corresponding to the delay value𝑇0. It thus follows that as 𝑇increases from zero and crosses𝑇0=0.3427 a Hopf bifurcation occurs meaning thereby an initiation of periodic solution(s). One such periodic solution (limit cycle) of the model (2.4) is shown to exist for 𝑇=0.4 in Figure 2.

891095.fig.002
Figure 2: Limit cycle solution of the model (2.4) for 𝑇=0.4 starting from the initial point𝑥1(0)=1,𝑥2(0)=1,𝑥3(0)=1.

It has been well known from ecological model results (McDonald [49], Cushing [50], May [51]), especially for models with prey-predator-type interactions, that small delays in general enhance stability where as large delays in the growth response of the species may cause instability. In order to check if such possibility of instability of equilibrium occurs for this model also, we integrated model (2.4) numerically for large values of delay𝑇. It has been quite interesting to note that for large values of delay𝑇, model (2.4) showed irregular pattern in time series for each cell population. The fact that these irregular patterns are indeed chaotic in nature giving rise to chaotic attractors is confirmed by the sensitivity of the solutions to initial conditions. We present here only two illustrations of chaotic attractors for 𝑇=16 and 26 in Figures 3 and 4, respectively.

fig3
Figure 3: (a) Time series solution of the model (2.4) for 𝑇=16. (b) Corresponding chaotic attractor in the phase space. (c) Sensitivity of solution to initial conditions: 𝑥2=1.0,𝑥3=1.0 for both types of curves but 𝑥1=1.0 for solid curve (i.e., on the attractor) and 𝑥1=0.99 for dashed curve.
fig4
Figure 4: (a) Time series solution of the model (2.4) for 𝑇=26. (b) Corresponding chaotic attractor in the phase space. (c) Sensitivity of solution to initial conditions: 𝑥2=1.0,𝑥3=1.0 for both types of curves but 𝑥1=1.0 for solid curve and 𝑥1=0.99 for dashed curve.

5. Discussion and Conclusions

The response of the tumor diseases to treatment depends upon many factors including the severity of the tumor, the application of the treatment, and most importantly the patient’s immune response. Tumor cells are characterized by a vast number of genetic and epigenetic events leading to the appearance of specific antigens called neoantigens triggering antitumor activity by the immune system (El-Gohary [52], d’Onofrio [53]). Though this paper does not deal directly with any external treatment of the tumor, of course it focuses on the indirect treatment aspects of the disease by looking into the role of the immune system if it does not get triggered immediately but shows delayed responses. With this in mind, we modify the model of El-Gohary [32] to incorporate time delayed responses of the immune system through the growth mechanism of the hunting cells. It is assumed that hunting cells do not respond to killing of tumor cells as soon as they get signal from resting cells but they get activated after a constant time delay𝑇. This assumption yields our main model (2.4) as a system of delay differential equations. As has been mentioned in the introduction, the model of this paper represents a link between the super-macro-scale (in terms of ordinary differential equations) and the lower cellular scale (in terms of delay).

The dynamics and the stability results of the model show three main patterns of solutions: (i) stable equilibrium, (ii) limit cycle solution and (iii) chaotic attractor. More specifically, it is found that when hunting cells are either all time-alert (𝑇0) or alert enough (0<𝑇<𝑇0), all three cell populations approach to equilibrium values and the tumor can be said to be nonmalignant. For averagely alert hunting cells (i.e., when 𝑇𝑇0 or slightly greater than𝑇0), all the three cell populations may coexist in a limit cycle or periodic solution. In this case, the tumor can be termed as mildly malignant. The existence of periodic solutions is relevant in cancer models. It implies that the tumor levels may oscillate around a fixed point even in the absence of any treatment. Such a phenomenon, known as Jeff’s phenomenon or self-regression of tumor [54], has been observed clinically. When the hunting cells play too lethargic in their response to killing of tumor cells (i.e., when 𝑇 is large enough), all the three cell populations may grow in an irregular fashion with time leading to chaotic attractors. This is indeed the case when the tumor can be said to be malignant and it is the case where a serious treatment strategy is required because of continuously changing density of tumor cells all the time.

It is well known from ecological-model results that large delays cause instability of equilibriums. Thus one can say that the results of the present model are on the known lines but we feel that instability in the form of chaotic attractors in cancer modeling is quite an interesting observation of this study linking super-macro-scale to lower cellular scale. The allowable time delay for activation of the immune system and the estimation of the length of delay to preserve stability may be the two important parameters that may help decide the mode of action for controlling the disease.

Acknowledgment

The author (M. Saleem) acknowledges the financial support of the UGC, New Delhi, for this work under the Project no. 37-483/2009(SR).

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
  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, pp. 73–79, 2010.
  4. R. de Boer, P. Hoogeweg, H. Dullens, R. de Weger, and W. Den Otter, “Macrophage T lymphocyte interactions in the anti-tumor immune response: a mathematical model,” Journal of Immunology, vol. 134, no. 4, pp. 2748–2758, 1985.
  5. B. Goldstein, J. R. Faeder, and W. S. Hlavacek, “Mathematical and computational models of immune-receptor signalling,” Nature Reviews Immunology, vol. 4, no. 6, pp. 445–456, 2004. View at Scopus
  6. 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 Publisher · View at Google Scholar · View at Scopus
  7. N. Kronik, Y. Kogan, V. Vainstein, and Z. Agur, “Improving alloreactive CTL immunotherapy for malignant gliomas using a simulation model of their interactive dynamics,” Cancer Immunology, Immunotherapy, vol. 57, no. 3, pp. 425–439, 2008. View at Publisher · View at Google Scholar · View at Scopus
  8. J. Aroesty, T. Lincoln, N. Shapiro, and G. Boccia, “Tumor growth and chemotherapy: mathematical methods, computer simulations, and experimental foundations,” Mathematical Biosciences, vol. 17, no. 3-4, pp. 243–300, 1973. View at Publisher · View at Google Scholar · View at Scopus
  9. M. Eisen, Mathematical Models in Cell Biology and Cancer Chemotherapy, vol. 30 of Lecture Notes in Biomathematics, Springer, Berlin, Germany, 1979.
  10. H. Knolle, Cell Kinetic Modelling and the Chemotherapy of Cancer, vol. 75 of Lecture Notes in Biomathematics, Springer, Berlin, Germany, 1988.
  11. J. M. Murray, “Optimal control for a cancer chemotherapy problem with general growth and loss functions,” Mathematical Biosciences, vol. 98, no. 2, pp. 273–287, 1990. View at Publisher · View at Google Scholar
  12. J. A. Adam, “The dynamics of growth-factor-modified immune response to cancer growth: one-dimensional models,” Mathematical and Computer Modelling, vol. 17, no. 3, pp. 83–106, 1993. View at Publisher · View at Google Scholar
  13. J. Adam and J. C. Panetta, “A simple mathematical model and alternative paradigm for certain chemotherapeutic regimens,” Mathematical and Computer Modelling, vol. 22, no. 8, pp. 49–60, 1995. View at Publisher · View at Google Scholar
  14. M. R. Owen and J. A. Sherratt, “Modelling the macrophage invasion of tumors: effects on growth and composition,” IMA Journal of Mathemathics Applied in Medicine and Biology, vol. 15, no. 2, pp. 165–185, 1998. View at Scopus
  15. L. G. De Pillis and A. Radunskaya, “A mathematical tumor model with immune resistance and drug therapy: an optimal control approach,” Journal of Theoretical Medicine, vol. 3, no. 2, pp. 79–100, 2001. View at Scopus
  16. D. Dingli, M. D. Cascino, K. Josić, S. J. Russell, and Z. Bajzer, “Mathematical modeling of cancer radiovirotherapy,” Mathematical Biosciences, vol. 199, no. 1, pp. 55–78, 2006. View at Publisher · View at Google Scholar
  17. S. A. Menchón, R. A. Ramos, and C. A. Condat, “Modeling subspecies and the tumor-immune system interaction: steps toward understanding therapy,” Physica A, vol. 386, no. 2, pp. 713–719, 2007. View at Publisher · View at Google Scholar · View at Scopus
  18. N. Bellomo and L. Preziosi, “Modelling and mathematical problems related to tumor evolution and its interaction with the immune system,” Mathematical and Computer Modelling, vol. 32, no. 3-4, pp. 413–452, 2000. View at Publisher · View at Google Scholar
  19. R. P. Araujo and D. L. S. McElwain, “A history of the study of solid tumour growth: the contribution of mathematical modelling,” Bulletin of Mathematical Biology, vol. 66, no. 5, pp. 1039–1091, 2004. View at Publisher · View at Google Scholar
  20. J. D. Nagy, “The ecology and evolutionary biology of cancer: a review of mathematical models of necrosis and tumor cell diversity,” Mathematical Biosciences and Engineering, vol. 2, no. 2, pp. 381–418, 2005.
  21. H. M. Byrne, T. Alarcon, M. R. Owen, S. D. Webb, and P. K. Maini, “Modelling aspects of cancer dynamics: a review,” Philosophical Transactions of the Royal Society of London. Series A, vol. 364, no. 1843, pp. 1563–1578, 2006. View at Publisher · View at Google Scholar
  22. F. Castiglione and B. Piccoli, “Cancer immunotherapy, mathematical modeling and optimal control,” Journal of Theoretical Biology, vol. 247, no. 4, pp. 723–732, 2007. View at Publisher · View at Google Scholar
  23. M. L. Martins, S. C. Ferreira Jr., and M. J. Vilela, “Multiscale models for the growth of avascular tumors,” Physics of Life Reviews, vol. 4, no. 2, pp. 128–156, 2007. View at Publisher · View at Google Scholar · View at Scopus
  24. T. Roose, S. J. Chapman, and P. K. Maini, “Mathematical models of avascular tumor growth,” SIAM Review, vol. 49, no. 2, pp. 179–208, 2007. View at Publisher · View at Google Scholar
  25. M. A. J. Chaplain, “Modelling aspects of cancer growth: insight from mathematical and numerical analysis and computational simulation,” in Multiscale Problems in the Life Sciences, vol. 1940 of Lecture Notes in Math., pp. 147–200, Springer, Berlin, Germany, 2008.
  26. 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
  27. 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
  28. 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
  29. 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 Scopus
  30. 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 Scopus
  31. 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 · View at Google Scholar
  32. 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 · View at Google Scholar
  33. 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
  34. 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
  35. 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
  36. J. Miękisz, J. Poleszczuk, M. Bodnar, and U. Foryś, “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
  37. 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
  38. 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.
  39. 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.
  40. 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
  41. N. Burić, M. Mudrinic, and N. Vasović, “Time delay in a basic model of the immune response,” Chaos, solitons and fractals, vol. 12, no. 3, pp. 483–489, 2001. View at Publisher · View at Google Scholar · View at Scopus
  42. U. Foryś and M. Kolev, “Time delays in proliferation and apoptosis for solid avascular tumor,” Prep. institute of applied mathematics and mechanics, rw02-10 (110), Warsaw University, 2002.
  43. 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.
  44. R. Yafia, “Dynamics analysis and limit cycle in a delayed model for tumor growth with quiescence,” Nonlinear Analysis: Modelling and Control, vol. 11, no. 1, pp. 95–110, 2006.
  45. 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
  46. 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
  47. H. I. Freedman and V. S. H. Rao, “The trade-off between mutual interference and time lags in predator-prey systems,” Bulletin of Mathematical Biology, vol. 45, no. 6, pp. 991–1004, 1983. View at Publisher · View at Google Scholar
  48. L. H. Erbe, H. I. Freedman, and V. Sree Hari Rao, “Three-species food-chain models with mutual interference and time delays,” Mathematical Biosciences, vol. 80, no. 1, pp. 57–80, 1986. View at Publisher · View at Google Scholar
  49. N. MacDonald, Time Lags in Biological Models, vol. 27 of Lecture Notes in Biomathematics, Springer, Berlin, Germany, 1978.
  50. J. M. Cushing, Integrodifferential Equations and Delay Models in Population Dynamics, vol. 20 of Lecture Notes in Biomathematics, Springer, Berlin, Germany, 1977.
  51. R. M. May, Stability and Complexity in Model Ecosystems, Princeton University Press, Princeton, NJ, USA, 1973.
  52. A. A. El-Gohary, “Optimal control of the genital herpes epidemic,” Chaos, Solitons and Fractals, vol. 12, no. 10, pp. 1817–1822, 2001. View at Publisher · View at Google Scholar · View at Scopus
  53. A. d'Onofrio, “Tumor evasion from immune control: strategies of a MISS to become a MASS,” Chaos, Solitons and Fractals, vol. 31, no. 2, pp. 261–268, 2007. View at Publisher · View at Google Scholar
  54. R. H. Thomlinson, “Measurement and management of carcinoma of the breast,” Clinical Radiology, vol. 33, no. 5, pp. 481–493, 1982. View at Scopus