Abstract
We study a cancer model given by a three-dimensional system of ordinary differential equations, depending on eight parameters, which describe the interaction among healthy cells, tumour cells, and effector cells of immune system. The model was previously studied in the literature and was shown to have a chaotic attractor. In this paper we study how such a chaotic attractor is formed. More precisely, by varying one of the parameters, we prove that a supercritical Hopf bifurcation occurs, leading to the creation of a stable limit cycle. Then studying the continuation of this limit cycle we numerically found a cascade of period-doubling bifurcations which leads to the formation of the mentioned chaotic attractor. Moreover, analyzing the model dynamics from a biological point of view, we notice the possibility of both the tumour cells and the immune system cells to vanish and only the healthy cells survive, suggesting the possibility of cure, since the interactions with the immune system can eliminate tumour cells.
1. Introduction
Mathematical models describing cancer tumour growth have been extensively studied in the literature in order to understand the mechanism of disease and to predict its future behaviour; see [1] and references therein as a brief review on the subject. One class of these cancer models, well known as Lotka-Volterra systems, is based on nonlinear mathematical systems of interacting species.
The competition for nutrients, the effects of growth factors, and the actions of the immune system stand out among possible interactions.
The immune system has the function of recognizing internal and external threats to the organism, reacting to eliminate, neutralize, or tolerate such threats. The recognition of tumour cells by the immune system can happen in distinct and complementary ways. According to Chammas et al. [2], the lymphocytes of type B participate in the immune response against tumours by producing specific antibodies against tumour antigens. These antibodies bind with specific antigens on tumour cells which facilitate their recognition and destruction by NK cells and phagocytosis by macrophages. The study of different antibodies found in patients with cancer tumour has been useful for determining tumour-associated antigens and pointed out some information for the development of therapeutic antibodies with clinical application. More specifically, several researches seeking the elimination of tumour cells through the immune cells have been developed [3–5]. Despite the increasing number of researches in this area, in [1] the authors pointed out the existence of a gap in the understanding of tumour growth and invasion, requiring an intense and closer collaboration in investigations by different scientists, such as mathematicians, biologists, and experimentalists.
Aiming to contribute to the understanding of these types of model we perform a bifurcation analysis of a cancer tumour growth model of three competing cell populations (tumour, healthy, and immune cells), which was proposed in [6]. The populations of these cells in a given instant of time will be denoted here by , , and , respectively. Analogous to what is done in competition models in population ecology, the mathematical model for the interaction of these cells is given by the following three-dimensional adimensional system of ordinary differential equations (for details on the derivation of the system, see the nice paper [6]),in which the parameters , , , , , , , are all positive.
The equations of system (1) describe the variation rate of tumour cells, healthy cells, and effector cells with the time . In the first equation, a logistic growth of tumour cells is considered, the term represents the negative effects that the healthy cells have on tumour cells, and represents the negative effects that the effector cells of the immune system exert on tumour cells. In the second equation, the healthy cells also grow logistically, with growth rate . The term represents the negative effects that tumour cells exert on healthy cells. In the last equation, describes the stimulation of the immune system by tumour cells through specific antigens. Therefore, the stimulus of the immune system depends directly on the number of tumour cells. The term represents the negative effects that tumour cells exert on the effector cells of the immune system and, finally, refers to the natural death rate of the effector cells. It is important to notice that system (1) is a slight modification of a model of cancer without treatment, proposed in [7], where a first phase space analysis of the model was presented.
In [6], by using the Shilnikov’s theorem [8] and by the calculation of Lyapunov exponents and Lyapunov dimension, the authors have shown that system (1) presents a chaotic attractor for the following fixed set of parameters: , , , , , , , and . In a second paper concerning the dynamical analysis of system (1) [9], the authors described in detail the topological structure of this chaotic attractor. However, in none of these papers the authors show how this attractor is formed nor do they relate the chaotic behaviour of system (1) with biological aspects of the cancer evolution.
In this paper we analyse the local stability of the equilibria of system (1); then, varying two of the parameters involved, namely, the parameters and , and taking for the other parameters the same values considered in [6], we perform a bifurcation analysis of system (1). Through this analysis we prove that, for fixed and varying the parameter , the system presents a supercritical Hopf bifurcation at an equilibrium point for the critical value , which leads to the creation of a small stable limit cycle. Furthermore, numerically studying the continuation of this stable limit cycle by increasing the parameter from the Hopf point, we numerically found the occurrence of a cascade of period-doubling bifurcations which leads to the formation of the chaotic attractor shown to exist in [6]. Trying to analyse this dynamical behaviour from the biological point of view, we observed that, after the occurrence of chaotic dynamics, both the tumour cells and the immune system cells (represented by the coordinates and ) vanish and only the healthy cells remain positive and their amount tends to the carrying capacity . The chaotic dynamics occurs when the parameter varies from the critical value to and the dynamics is not affected by the variation of parameter , which must only be taken different from zero, as will be shown in the calculations ahead. In this way, the bifurcation analysis presented here is performed focusing on the variation of parameter .
The parameter considered in the mentioned bifurcation analysis represents the interaction between the cancer cells and the healthy cells. However, this parameter can also be related to the action of the immune system cells against cancer cells. In fact, we can rewrite the first equation of system (1) asIn this way, for smaller than and fixed, the negative effect of the immune system cells onto the tumour cells is more significant and it decreases as approaches the critical value , when the chaotic dynamics occurs. Although being a very complex matter, this and another tentative analysis relating the dynamical behaviour and the biological meaning of the variables of system (1) will be further considered later on, in Sections 2, 4, and 5.
We believe that the bifurcation analysis and the biological considerations presented here complement the results and the analysis of system (1) performed in [6, 7, 9].
The paper is organized as follows. In this introductory section, we present the studied model, describing the variables and parameters involved, and state the main results obtained. In Section 2, we make a linear analysis of system (1): calculate the equilibrium points, check whether these equilibria are biologically admissible (i.e., if all their coordinates are nonnegative), analyse their local stability, and give an outline of the Hopf bifurcation theorem. In Section 3 we prove the main result of the paper, the occurrence of a supercritical Hopf bifurcation in one of the equilibrium points of system (1), which gives rise to a small stable limit cycle; we also present some numerical simulations which corroborate with the occurrence of this bifurcation. In Section 4 we numerically show that a cascade of period-doubling bifurcations occurs with the limit cycle created at the Hopf bifurcation, leading to the creation of a chaotic attractor (exactly that one shown to exist in [6]). Finally, Section 5 presents some concluding remarks, including a possible interpretation of the dynamical properties of system (1) from a biological point of view.
2. Linear Analysis
Considering in system (1), we obtain the following equilibrium points (for details on the calculation and on the local analysis of the equilibrium points of system (1), see [6]):
These equilibria are called biologically admissible if their three coordinates are greater than or equal to zero. The equilibrium points , , and satisfy such condition.
For the equilibrium point we have the following possibilities:(i)if , the equilibrium coincides with the equilibrium ;(ii)if , coincides with the equilibrium ;(iii)if , , and , is biologically admissible;(iv)if , , and , is also admissible;(v)if none of the previous inequalities are satisfied, is not biologically admissible.
As the equilibria , and , may present negative coordinates, aiming to simplify the analysis of system (1), we will present a study of the feasibility of these equilibria using the values of parameters extracted from [6], that is, , , , , , and , and letting the parameters and vary. For these parameter values we obtain Note that the equilibria , , remain the same and and are not biologically admissible; therefore it will not be considered here. Let us analyse the conditions on the equilibria , , and . (i) is biologically admissible, since we are considering ;(ii)the equilibrium exists if and coincides with the equilibrium for , as seen previously. If , is biologically admissible, which does not occur if ;(iii) is biologically admissible if , since .
2.1. Linear Stability of the Equilibrium Points
The Jacobian matrix of system (1) at a point is given byIt will be used bellow to study the local stability of the admissible equilibria.(i)Equilibrium : the Jacobian matrix applied at has the eigenvaluesAs and , is (unstable) saddle point, which is coherent from the biological point of view, since this point represents the simultaneous extinction of normal cells, immune cells, and tumour cells.(ii)Equilibrium : the eigenvalues of the matrix at this point areThus, the stability of depends on the value of , whereas the other two eigenvalues are negative. If , that is, if , then the equilibrium is a saddle point, therefore unstable; if , is asymptotically stable (a stable node); and if , is nonhyperbolic, since . The equilibrium point itself represents the absence of cancer and immune system cells and the permanence of the healthy cells, which tend to its carrying capacity . In this way, from the biological point of view the instability of (case ) represents the existence of cancer cells and, consequently, the existence of immune system cells that is the persistence of the disease; on the other hand, the stability of represents the absence of disease, which could be possible through a given treatment or by the action of the immune system. As the model considered does not admit treatment and includes the action of immune system cells, we expect that these cells somehow combat the tumour cells, thus enabling the cure, which justifies the possible stability of the equilibrium point (case ).(iii)Equilibrium : the Jacobian matrix at has the eigenvaluesAs , is a saddle point; therefore it is unstable, which is expected from the biological standpoint, since this equilibrium point represents the extinction of healthy and immune cells, remaining only tumour cells.(iv)Equilibrium : the Jacobian matrix calculated at has the eigenvaluesfrom which follows that is a saddle-focus. As , the equilibrium is unstable, which could be possible from the biological point of view, since this equilibrium represents the extinction of healthy cells and the persistence of tumour cells and immune system cells.(v)Equilibrium : the Jacobian matrix at this point has two eigenvalues of the formwith The third eigenvalue is given by Recalling that is biologically admissible if , we have . Hence is a saddle point. The instability of represents the biological fact that, for the set of parameters adopted, tumour cells and healthy cells cannot coexist without the presence of immune cells; that is, the tumour cells actually activate the immune system cells.(vi)Equilibrium : the Jacobian matrix calculated at the point has the following characteristic polynomialwhere , , and . Hence the linear stability of depends on the value of parameter .
We observe that the equilibrium point plays an important role in our analysis, since it represents the coexistence of the three types of cells considered in the model. Thus a detailed analysis of its stability will be presented here. The following result will be used in such analysis (for a proof, see [10]).
Lemma 1. Consider the cubic polynomialwith real coefficients. (a)Then (Routh-Hurwitz Criterion) is stable (i.e., all its roots have negative real part) if and only if(b)If and , then has two roots with positive real part and one negative real root.(c)If and , then has two complex conjugate roots with zero real part and one negative real root.
Applying the previous lemma to the polynomial (13), we can obtain the equilibrium as a candidate for the occurrence of a Hopf bifurcation in system (1). We will prove that this bifurcation indeed occurs at the point in the next section. Before doing this, let us remember some results about the Hopf bifurcation theory.
2.2. Outline of the Hopf Bifurcation Theory and the Projection Method
This section is a brief review of the projection method described in [8] for the calculation of the first Lyapunov coefficient associated with the well-known Hopf bifurcation, denoted by .
Consider the differential equationwhere and are, respectively, vectors representing phase variables and control parameters. Assume that is of class in . Suppose that (16) has an equilibrium point at and, denoting the variable also by , writeas where and, for , and so on for , , , and .
Suppose that is an equilibrium point of (16) where the Jacobian matrix has a pair of purely imaginary eigenvalues , , and admits no other eigenvalue with zero real part. Let be the generalized eigenspace of corresponding to . By this it is meant the largest subspace invariant by on which the eigenvalues are .
Let be vectors such thatwhere is the transpose of the matrix . Any vector can be represented as , where . The two dimensional center manifold associated with the eigenvalues can be parameterized by the variables and by means of an immersion of the form , where has a Taylor expansion of the formwith and . Substituting this expression into (16) we obtain the following differential equation:where is given by (17). The complex vectors are obtained solving the system of linear equations defined by the coefficients of (22), taking into account the coefficients of , so that system (22), on the chart for a central manifold, is written as follows:with .
The first Lyapunov coefficient is defined bywhere and .
A Hopf point of system (16) is an equilibrium point where the Jacobian matrix has a pair of purely imaginary eigenvalues , and the other eigenvalue . From the Center Manifold Theorem, at a Hopf point a two-dimensional center manifold is well-defined, and it is invariant under the flow generated by (16) and can be continued with arbitrary high class of differentiability to nearby parameter values (see [8]).
A Hopf point is called transversal if the parameter dependent complex eigenvalues cross the imaginary axis with nonzero derivative. In a neighborhood of a transversal Hopf point with the dynamic behaviour of the system (16), reduced to the family of parameter-dependent continuations of the center manifold, is orbitally topologically equivalent to the following complex normal form: where and , , and are real functions having derivatives of arbitrary higher order, which are continuations of , , and the first Lyapunov coefficient at the Hopf point. See [8] for details. When () one family of stable (unstable) periodic orbits can be found on this family of manifolds, shrinking to an equilibrium point at the Hopf point.
Based on these results, a numerical algorithm for calculating the first Lyapunov coefficient of a system with and can be obtained. This algorithm has been used to make the calculations of the next section.
3. Hopf Bifurcation at the Point
By using the Hopf bifurcation theorem stated in [8] and the projection method presented at the same reference and synthesized in the previous section we prove the following result.
Theorem 2. Consider system (1) with the following parameter values: , , , , , , and . If , then the equilibrium point is a stable focus of system (1). If is a weak stable focus (or a vague attractor). On the other hand, for , and close to this value, the equilibrium becomes an unstable focus and a small stable limit cycle arises around . In short, a supercritical Hopf bifurcation occurs at the equilibrium point for the critical value .
Proof. Considering the characteristic polynomial of the Jacobian matrix calculated at the point , given by (13), we have (1);(2) if ;(3) if ;(4) if .
The necessary condition for the point being biologically admissible is ; therefore is always positive. If , then and ; then, according to the Routh-Hurwitz Criterion, is a stable focus. On the other hand, if , then and . Hence, from item (b) of Lemma 1 it follows that (13) has two roots with positive real part and a negative real root. Thus is an unstable focus, remaining unstable for or .
If , then and ; hence from item (c) of Lemma 1, the polynomial (13) has two complex conjugate roots with zero real part and a negative real root, which is one of the conditions for the occurrence of the Hopf bifurcation. In order to guarantee that the Hopf bifurcation actually occurs at the equilibrium and determine the stability of the limit cycle originated from it, we have to calculate the first Lyapunov coefficient of system (1) associated with this point. Using the projection method presented in the previous section, we obtained the value for the first Lyapunov coefficient of system (1) at the point : To finish the proof, the transversality condition remains to be checked. Using software MAPLE, we calculate the derivative of the Jacobian matrix with respect to the parameter applying the critical value from which we find Solving the equation we obtain This ends the proof of Theorem 2.
The numerical simulations presented in Figures 1 and 2 illustrate the results stated in Theorem 2.
(a)
(b)
(a)
(b)
The existence of a stable limit cycle implies the occurrence of a periodic behaviour of system (1); that is, from the biological point of view, healthy cells, tumour cells, and immune system cells can coexist.
In the next section, by studying the continuation of this limit cycle when the parameter moves away from the critical value , we numerically find a cascade of period-doubling bifurcations initiated with the limit cycle created in the Hopf bifurcation and leading to the creation of a chaotic attractor, when tends to the limit value .
4. Period-Doubling Bifurcations and Chaotic Dynamics
For close to and greater than the critical value , for which the Hopf bifurcation at the equilibrium occurs, we have a stable limit cycle of period (see Figure 3); increasing the value of , we numerically found that this limit cycle splits into a period periodic orbit (see Figure 4). When the value of increases a bit more, a periodic orbit is observed (see Figure 5); increasing further the value of , a series of duplications and a limit cycle of period are observed, until a chaotic attractor is formed for (see Figure 6(a)). This type of bifurcation is called a period-doubling or “flip bifurcation.’’ System (1) presents therefore an infinite cascade of flip bifurcations with a finite accumulation point; the associated dynamics at this limit point is chaotic; that is, the cascade of period-doubling bifurcations is a continuous transition from periodic solutions to chaos (see Figures 3–6).
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
Figure 6(b) shows the sensitivity to initial conditions of two solutions of system (1) with very close initial conditions, confirming the chaotic behaviour of the system, for these parameter values.
A bifurcation diagram confirming the occurrence of the period-doubling cascade mentioned above, culminating in chaos, is shown in Figure 7.
4.1. The Elimination of Tumour and Immune System Cells after the Chaotic Dynamics
The chaotic dynamics occurs when the value of is approximately . In the adimensional system (1) the parameter , which represents the interaction between tumour cells and healthy cells, can be obtained in terms of the parameters in the original model (see [6]), given byin the following way: from which we have Thus, for small values of parameter we have large values of , since is the largest parameter of the original system, because it represents the carrying capacity of healthy cells. Then the value of , which represents the growth rate of tumour cells, is large when . This explains the fact that the tumour-free equilibrium point is unstable, since, as the growth rate of the tumour cells is large, the tumour cells are increasing when . On the other hand, the effector cells of the immune system also grow because their growth is proportional to the growth of tumour cells; then the number of effector cells is also large when . Thus, we assume that this growth of tumour cells and effector cells, along with the interactions of healthy cells, results in a dispute whose apex occurs for . This analysis suggests that the apex of the battle between cancer cells, healthy cells, and immune system cells is represented by the chaotic behaviour of the system.
In turn, for , the value of becomes smaller, so the number of tumour cells decreases, which implies in a decreasing in the number of effector cells of the immune system, since its variation rate is directly proportional to the number of tumour cells (see the last equation in (30)). In this way, it is possible that the interactions with the immune system eliminate the tumour cells. Thus, if there are no more tumour cells, there will be no immune cells, leaving only healthy cells, which makes the tumour-free equilibrium point stable, allowing a possible cure of the disease for the set of parameters considered in this analysis.
In Figure 8 we present some numerical simulations of the solutions of system (1) with , showing the stability of the equilibrium point , illustrating the fact that, for the set of parameters considered here, only healthy cells survive after the battle between tumour cells and immune cells.
5. Concluding Remarks
In this work, using the concepts of qualitative theory, bifurcations, and chaos in dynamical systems, we analyze a three-dimensional cancer model, describing the interactions among healthy cells, tumour cells, and immune system cells. The analyzed model was proposed in [6], where the existence of chaotic dynamics for a fixed set of parameter values was proved. Here we performed a bifurcation analysis by varying one of the parameters involved in the system, showing that a Hopf bifurcation occurs, giving rise to a stable limit cycle, which represents the possibility of coexistence of tumour cells, immune cells, and healthy cells. Analyzing the continuation of this limit cycle numerically, we find a cascade of period-doubling bifurcations, leading to the creation of a chaotic attractor, which is exactly the one shown to exist in [6].
In a first glance one could argue that, from the biological point of view, the chaotic dynamics could be related to the process of uncontrolled growth of tumour cells. However, we surprisingly observed that, for the set of parameter values considered the period-doubling bifurcations which occur to the cancer cells lead to the growth of effector cells, since their growth is proportional to the growth of tumour cells. In this way, it seems that, for the particular values of the parameter considered here, the growth of both types of cells occurs simultaneously and is represented mathematically by the period-doubling cascade, which leads to the creation of a chaotic attractor. Furthermore, supposing the increasing values of the parameter , then we observe that after the chaotic dynamics (i.e., for ) both tumour cells and immune system cells vanish and the healthy cells tend to their carrying capacity. It indicates that in the studied model the chaotic dynamics could be considered as a kind of mechanism related to a possible cure of the cancer, through the elimination of tumour cells, at least for the set of parameter values considered. This situation is obviously very optimistic but it is mathematically possible in the studied model. Of course, the model could be studied with other parameter values, in order to detect other dynamical behaviour, which could represent biological aspects different from the one presented here. Indeed, after the acceptance of this paper for publication, we knew about the paper by Lettelier et al. [11] concerning the cancer model studied in this note, which can be of interest of the readers.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was partially supported by PROPG/UNESP. The third author is supported by CNPq-Brazil, under the Project 308315/2012-0 and by FAPESP-Brazil under the Process 2012/18413-7. The authors would like to thank the referees for their valuable comments and suggestions, which helped them to improve the presentation of the paper.