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 [33–39]. 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: 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,:maximum carrying capacity of tumor cells,:rate of killing of tumor cells by hunting cells,:conversion rate of the resting cells to tumor cells,:natural death rate of hunting cells,:growth rate of resting cells,:natural death rate of resting cells.
Using nondimensional variables and parameters as El-Gohary [32] obtained the following nondimensional form of model (2.1):
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 in (2.3) by . Thus our model takes the form It is obvious that model (2.4) would have same equilibrium solutions as model (2.3), whose equilibriums have been reported as in El-Gohary [32], under the biologically feasible conditions as 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., )
The local stability results for equilibriums of model (2.4) for have been reported in El-Gohary [32] as shown in (Table 1).
3.2. Stability with Delay (i.e., )
Assuming small deviations around the equilibrium such that , the linearized system of the model (2.4) becomes In the case of a positive delay , the characteristic equation for this system can be written as where such that with Now substituting (where is positive) in (3.2) and separating the real and imaginary parts, we obtain the following system of transcendental equations: Squaring and adding (3.6) and (3.7), we get which can be simplified to where Equation (3.9) can be written as a cubic with .
For parameter values such that is positive, the simplest assumption that (3.11) will have a unique positive root is . Since is positive, it requires that for to be negative. Hence, it can be said that there is a unique positive root say of (3.11). Denoting , it follows that the characteristic equation (3.2) has a pair of purely imaginary roots of the form . Eliminating from (3.6) and (3.7), we get Then corresponding to is given by Since is stable for , it implies from Freedman and Rao [47] that remains stable for .
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 where 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 are given by where and is the smallest positive root of (3.17).
Inequality (3.16) and (3.17) can alternatively be written as 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 , independent of , from (3.19) and then to estimate so that (3.18) holds true for all values of such that and hence, in particular, at .
Equation (3.19) is rewritten as Maximizing the right-hand side of (3.20)
subject to we obtain Hence if then clearly from (3.22) we have .
From (3.18), we obtain Since is locally asymptotically stable for , the inequality (3.24) will continue to hold for sufficiently small . Using (3.20), (3.24) can be rearranged as Using the bound we obtain from (3.24) where Hence, if then for 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 , and , 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 . It can be seen that for these parameter values the coefficients of the cubic (3.11) are such that and . It guarantees that the cubic (3.11) has a unique positive root , which then provides the unique positive root of (3.9) as . Using this value of , it turns out from (3.13) that , and the stability result of Section 3 yields that the positive equilibrium of model (2.4) is stable for such that . For , Figure 1 illustrates the approach of the trajectory of the model (2.4) to the equilibrium .
Indeed the existence of 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. It thus follows that as increases from zero and crosses 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 in Figure 2.
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 and 26 in Figures 3 and 4, respectively.
(a)
(b)
(c)
(a)
(b)
(c)
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 ( or alert enough (, 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 or slightly greater than), 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).