# Hopf Bifurcation, Cascade of Period-Doubling, Chaos, and the Possibility of Cure in a 3D Cancer Model

Authors: Marluci Cristina Galindo; Cristiane Nespoli; Marcelo Messias
Journal: Abstract and Applied Analysis
 (2015)

---

## 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.


---

## Body

## 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 t will be denoted here by x(t), y(t), and z(t), 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]),(1)dxdt=x(1-x)-axy-bxz,dydt=cy(1-y)-dxy,dzdt=exzx+f-gxz-hz,in which the parameters a, b, c, d, e, f, g, h are all positive.The equations of system (1) describe the variation rate of tumour cells, healthy cells, and effector cells with the time t. In the first equation, a logistic growth of tumour cells is considered, the term axy represents the negative effects that the healthy cells have on tumour cells, and bxz 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 c. The term dxy represents the negative effects that tumour cells exert on healthy cells. In the last equation, exz/(x+f) 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 gxz represents the negative effects that tumour cells exert on the effector cells of the immune system and, finally, h 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: a=1.0, b=2.5, c=0.6, d=1.5, e=4.5, f=1.0, g=0.2, and h=0.5. 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 a and b, 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 b>0 fixed and varying the parameter a, the system presents a supercritical Hopf bifurcation at an equilibrium point for the critical value a≅0.717455, 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 a 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 x(t) and z(t)) vanish and only the healthy cells remain positive and their amount tends to the carrying capacity y(t)=1. The chaotic dynamics occurs when the parameter a varies from the critical value a≅0.717455 to a=1.0 and the dynamics is not affected by the variation of parameter b, 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 a.The parametera 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) as(2)dxdt=x(1-x)-axy+baz.In this way, for a smaller than 1.0 and b≠0 fixed, the negative effect of the immune system cells onto the tumour cells is more significant and it decreases as a approaches the critical value a=1.0, 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 Section2, 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

Consideringx˙=y˙=z˙=0 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]): (3)A=0,0,0,B=0,1,0,C=1,0,0,E=x¯,0,z¯,E1=x1¯,0,z1¯,F=ca-1ad-c,d-cad-c,0,G=x2¯,y2¯,z2¯,G1=(x3¯,y3¯,z3¯).These equilibria are called biologically admissible if their three coordinates are greater than or equal to zero. The equilibrium pointsA, B, and C satisfy such condition.For the equilibrium pointF we have the following possibilities:(i)
ifa=1, the equilibrium F coincides with the equilibrium B;(ii)
ifc=d, F coincides with the equilibrium C;(iii)
ifa>1, ad>c, and d>c, F is biologically admissible;(iv)
ifa<1, ad<c, and d<c, F is also admissible;(v)
if none of the previous inequalities are satisfied,F is not biologically admissible.As the equilibriaE, E1 and G, G1 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, c=0.6, d=1.5, e=4.5, f=1, g=0.2, and h=0.5, and letting the parameters a and b vary. For these parameter values we obtain (4)E=0.132503,0,0.867497b,E1=18.8675,0,-17.867497b,F=2(a-1)5a-2,35a-2,0,G=0.132503,0.668742,0.1(8.675-6.687a)b,G1=18.8675,-46.1687,0.1(-178.675+461.687a)b.Note that the equilibria A, B, C remain the same and E1 and G1 are not biologically admissible; therefore it will not be considered here. Let us analyse the conditions on the equilibria E, F, and G. (i)
E is biologically admissible, since we are considering b>0;(ii)
the equilibriumF exists if a≠0.4 and coincides with the equilibrium B for a=1, as seen previously. If a>1, F is biologically admissible, which does not occur if a<1;(iii)
G is biologically admissible if a<1.2972, since b>0.


### 2.1. Linear Stability of the Equilibrium Points

The Jacobian matrix of system (1) at a point (x,y,z) is given by(5)J=1-2x-ay-bz-ax-bx-1.5y0.6-1.2y-1.5x00.1z(43-2x2-4x)x+1200.1(38x-2x2-5)x+1.It will be used bellow to study the local stability of the admissible equilibria.(i)
EquilibriumA=(0,0,0): the Jacobian matrix J applied at A has the eigenvalues(6)λ1=1,λ2=0.6,λ3=-0.5.As λ1>0 and λ2>0, A 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)
EquilibriumB=(0,1,0): the eigenvalues of the matrix J at this point are(7)λ1=1-a,λ2=-0.6,λ3=-0.5.Thus, the stability of B depends on the value of λ1, whereas the other two eigenvalues are negative. If λ1>0, that is, if a<1, then the equilibrium B is a saddle point, therefore unstable; if a>1, B is asymptotically stable (a stable node); and if a=1, B is nonhyperbolic, since λ1=0. The equilibrium point B itself represents the absence of cancer and immune system cells and the permanence of the healthy cells, which tend to its carrying capacity y(t)=1. In this way, from the biological point of view the instability of B (case a<1) 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 B 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 B (case a>1).(iii)
EquilibriumC=(1,0,0): the Jacobian matrix J at C has the eigenvalues(8)λ1=-1,λ2=-0.9,λ3=1.55.As λ3>0, C 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)
EquilibriumE=0.132503,0,0.867497/b: the Jacobian matrix J calculated at E has the eigenvalues(9)λ1=0.401245,λ2=-0.6625+0.613i,λ3=-0.6625-0.613i,from which follows that E is a saddle-focus. As λ1>0, the equilibrium E 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)
EquilibriumF=2(a-1)/5a-2,3/(5a-2),0: the Jacobian matrix J at this point has two eigenvalues of the form(10)λ1=α+β(7a-3)(5a-2),λ2=α-β(7a-3)(5a-2)with (11)α=-7.5a2+4.5a-0.6,β=5525a4-11850a3+9025a2-2880a+324.The third eigenvalue is given by (12)λ3=0.1(124-416a+247a2)(8-34a+35a2).Recalling that F is biologically admissible if a>1, we have λ1>0. Hence F is a saddle point. The instability of F 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)
EquilibriumG=0.132503,0.668742,0.1(8.675-6.687a)/b: the Jacobian matrix J calculated at the point G has the following characteristic polynomial(13)λ3+a1λ2+a2λ+a3=0,where a1=0.533748, a2=0.433476-0.426092a, and a3=0.152597-0.117636a. Hence the linear stability of G depends on the value of parameter a.We observe that the equilibrium pointG 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 polynomial(14)p(x)=x3+ax2+bx+c,with real coefficients. (a)
Then (Routh-Hurwitz Criterion)p is stable (i.e., all its roots have negative real part) if and only if(15)a,b,c>0,ab-c>0.(b)
Ifc>0 and ab-c<0, then p has two roots with positive real part and one negative real root.(c)
Ifa,b,c>0 and ab-c=0, then p 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 G as a candidate for the occurrence of a Hopf bifurcation in system (1). We will prove that this bifurcation indeed occurs at the point G 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 l1.Consider the differential equation(16)x′=f(x,ζ),where x∈R3 and ζ∈R3 are, respectively, vectors representing phase variables and control parameters. Assume that f is of class C∞ in R3×R3. Suppose that (16) has an equilibrium point x=x0 at ζ=ζ0 and, denoting the variable x-x0 also by x, write(17)F(x)=f(x,ζ0)as (18)Fx=Ax+12Bx,x+16Cx,x,x+124Dx,x,x,x+1120E(x,x,x,x,x)+1720Kx,x,x,x,x,x+15040Lx,x,x,x,x,x,x+Ox8,where A=fx(0,ζ0) and, for i=1,2,3, (19)Bi(x,y)=∑j,k=13‍∂2Fi(ξ)∂ξj∂ξkξ=0xjyk,Cix,y,z=∑j,k,l=13‍∂3Fiξ∂ξj∂ξk∂ξlξ=0xjykzl,and so on for Di, Ei, Ki, and Li.Suppose that(x0,ζ0) is an equilibrium point of (16) where the Jacobian matrix A has a pair of purely imaginary eigenvalues λ2,3=±iω0, ω0>0, and admits no other eigenvalue with zero real part. Let Tc be the generalized eigenspace of A corresponding to λ2,3. By this it is meant the largest subspace invariant by A on which the eigenvalues are λ2,3.Letp,q∈C3 be vectors such that(20)Aq=iω0q,A⊤p=-iω0p,p,q=∑i=13‍p¯iqi=1,where A⊤ is the transpose of the matrix A. Any vector y∈Tc can be represented as y=wq+w¯q¯, where w=〈p,y〉∈C. The two dimensional center manifold associated with the eigenvalues λ2,3=±iω0 can be parameterized by the variables w and w¯ by means of an immersion of the form x=H(w,w¯), where H:C2→R3 has a Taylor expansion of the form(21)H(w,w¯)=wq+w¯q¯+∑2≤j+k≤7‍1j!k!hjkwjw¯k+O(w8),with hjk∈C3 and hjk=h¯kj. Substituting this expression into (16) we obtain the following differential equation:(22)Hww′+Hw¯w¯′=F(H(w,w¯)),where F is given by (17). The complex vectors hij are obtained solving the system of linear equations defined by the coefficients of (22), taking into account the coefficients of F, so that system (22), on the chart w for a central manifold, is written as follows:(23)w′=iω0w+12G21ww2+112G32ww4+1144G43ww6+Ow8,with Gjk∈C.The first Lyapunov coefficientl1 is defined by(24)l1=12ReG21,where G21=〈p,H21〉 and H21=C(q,q,q¯)+B(q¯,h20)+2B(q,h11).AHopf point(x0,ζ0) of system (16) is an equilibrium point where the Jacobian matrix A has a pair of purely imaginary eigenvalues λ2,3=±iω0, ω0>0 and the other eigenvalue λ1≠0. 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 calledtransversal if the parameter dependent complex eigenvalues cross the imaginary axis with nonzero derivative. In a neighborhood of a transversal Hopf point with l1≠0 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: (25)w′=(η+iω)w+l1ww2,where w∈C and η, ω, and l1 are real functions having derivatives of arbitrary higher order, which are continuations of 0, ω0, and the first Lyapunov coefficient at the Hopf point. See [8] for details. When l1<0 (l1>0) 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 systemx˙=f(x,μ) with x∈R3 and μ∈R can be obtained. This algorithm has been used to make the calculations of the next section.


## 2.1. Linear Stability of the Equilibrium Points

The Jacobian matrix of system (1) at a point (x,y,z) is given by(5)J=1-2x-ay-bz-ax-bx-1.5y0.6-1.2y-1.5x00.1z(43-2x2-4x)x+1200.1(38x-2x2-5)x+1.It will be used bellow to study the local stability of the admissible equilibria.(i)
EquilibriumA=(0,0,0): the Jacobian matrix J applied at A has the eigenvalues(6)λ1=1,λ2=0.6,λ3=-0.5.As λ1>0 and λ2>0, A 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)
EquilibriumB=(0,1,0): the eigenvalues of the matrix J at this point are(7)λ1=1-a,λ2=-0.6,λ3=-0.5.Thus, the stability of B depends on the value of λ1, whereas the other two eigenvalues are negative. If λ1>0, that is, if a<1, then the equilibrium B is a saddle point, therefore unstable; if a>1, B is asymptotically stable (a stable node); and if a=1, B is nonhyperbolic, since λ1=0. The equilibrium point B itself represents the absence of cancer and immune system cells and the permanence of the healthy cells, which tend to its carrying capacity y(t)=1. In this way, from the biological point of view the instability of B (case a<1) 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 B 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 B (case a>1).(iii)
EquilibriumC=(1,0,0): the Jacobian matrix J at C has the eigenvalues(8)λ1=-1,λ2=-0.9,λ3=1.55.As λ3>0, C 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)
EquilibriumE=0.132503,0,0.867497/b: the Jacobian matrix J calculated at E has the eigenvalues(9)λ1=0.401245,λ2=-0.6625+0.613i,λ3=-0.6625-0.613i,from which follows that E is a saddle-focus. As λ1>0, the equilibrium E 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)
EquilibriumF=2(a-1)/5a-2,3/(5a-2),0: the Jacobian matrix J at this point has two eigenvalues of the form(10)λ1=α+β(7a-3)(5a-2),λ2=α-β(7a-3)(5a-2)with (11)α=-7.5a2+4.5a-0.6,β=5525a4-11850a3+9025a2-2880a+324.The third eigenvalue is given by (12)λ3=0.1(124-416a+247a2)(8-34a+35a2).Recalling that F is biologically admissible if a>1, we have λ1>0. Hence F is a saddle point. The instability of F 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)
EquilibriumG=0.132503,0.668742,0.1(8.675-6.687a)/b: the Jacobian matrix J calculated at the point G has the following characteristic polynomial(13)λ3+a1λ2+a2λ+a3=0,where a1=0.533748, a2=0.433476-0.426092a, and a3=0.152597-0.117636a. Hence the linear stability of G depends on the value of parameter a.We observe that the equilibrium pointG 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 polynomial(14)p(x)=x3+ax2+bx+c,with real coefficients. (a)
Then (Routh-Hurwitz Criterion)p is stable (i.e., all its roots have negative real part) if and only if(15)a,b,c>0,ab-c>0.(b)
Ifc>0 and ab-c<0, then p has two roots with positive real part and one negative real root.(c)
Ifa,b,c>0 and ab-c=0, then p 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 G as a candidate for the occurrence of a Hopf bifurcation in system (1). We will prove that this bifurcation indeed occurs at the point G 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 l1.Consider the differential equation(16)x′=f(x,ζ),where x∈R3 and ζ∈R3 are, respectively, vectors representing phase variables and control parameters. Assume that f is of class C∞ in R3×R3. Suppose that (16) has an equilibrium point x=x0 at ζ=ζ0 and, denoting the variable x-x0 also by x, write(17)F(x)=f(x,ζ0)as (18)Fx=Ax+12Bx,x+16Cx,x,x+124Dx,x,x,x+1120E(x,x,x,x,x)+1720Kx,x,x,x,x,x+15040Lx,x,x,x,x,x,x+Ox8,where A=fx(0,ζ0) and, for i=1,2,3, (19)Bi(x,y)=∑j,k=13‍∂2Fi(ξ)∂ξj∂ξkξ=0xjyk,Cix,y,z=∑j,k,l=13‍∂3Fiξ∂ξj∂ξk∂ξlξ=0xjykzl,and so on for Di, Ei, Ki, and Li.Suppose that(x0,ζ0) is an equilibrium point of (16) where the Jacobian matrix A has a pair of purely imaginary eigenvalues λ2,3=±iω0, ω0>0, and admits no other eigenvalue with zero real part. Let Tc be the generalized eigenspace of A corresponding to λ2,3. By this it is meant the largest subspace invariant by A on which the eigenvalues are λ2,3.Letp,q∈C3 be vectors such that(20)Aq=iω0q,A⊤p=-iω0p,p,q=∑i=13‍p¯iqi=1,where A⊤ is the transpose of the matrix A. Any vector y∈Tc can be represented as y=wq+w¯q¯, where w=〈p,y〉∈C. The two dimensional center manifold associated with the eigenvalues λ2,3=±iω0 can be parameterized by the variables w and w¯ by means of an immersion of the form x=H(w,w¯), where H:C2→R3 has a Taylor expansion of the form(21)H(w,w¯)=wq+w¯q¯+∑2≤j+k≤7‍1j!k!hjkwjw¯k+O(w8),with hjk∈C3 and hjk=h¯kj. Substituting this expression into (16) we obtain the following differential equation:(22)Hww′+Hw¯w¯′=F(H(w,w¯)),where F is given by (17). The complex vectors hij are obtained solving the system of linear equations defined by the coefficients of (22), taking into account the coefficients of F, so that system (22), on the chart w for a central manifold, is written as follows:(23)w′=iω0w+12G21ww2+112G32ww4+1144G43ww6+Ow8,with Gjk∈C.The first Lyapunov coefficientl1 is defined by(24)l1=12ReG21,where G21=〈p,H21〉 and H21=C(q,q,q¯)+B(q¯,h20)+2B(q,h11).AHopf point(x0,ζ0) of system (16) is an equilibrium point where the Jacobian matrix A has a pair of purely imaginary eigenvalues λ2,3=±iω0, ω0>0 and the other eigenvalue λ1≠0. 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 calledtransversal if the parameter dependent complex eigenvalues cross the imaginary axis with nonzero derivative. In a neighborhood of a transversal Hopf point with l1≠0 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: (25)w′=(η+iω)w+l1ww2,where w∈C and η, ω, and l1 are real functions having derivatives of arbitrary higher order, which are continuations of 0, ω0, and the first Lyapunov coefficient at the Hopf point. See [8] for details. When l1<0 (l1>0) 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 systemx˙=f(x,μ) with x∈R3 and μ∈R can be obtained. This algorithm has been used to make the calculations of the next section.


## 3. Hopf Bifurcation at the PointG

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: b=2, c=0.6, d=1.5, e=4.5, f=1, g=0.2, and h=0.5. If a<0.717455, then the equilibrium point G is a stable focus of system (1). If a≅0.717455,G is a weak stable focus (or a vague attractor). On the other hand, for a>0.717455, and close to this value, the equilibrium G becomes an unstable focus and a small stable limit cycle arises around G. In short, a supercritical Hopf bifurcation occurs at the equilibrium point G for the critical value a=0.717455.Proof.
Considering the characteristic polynomial of the Jacobian matrixJ calculated at the point G, given by (13), we have (1)
a1>0;(2)
a2>0 if a<1.017;(3)
a3>0 if a<1.297;(4)
a1a2-a3>0 if a<0.717455.
The necessary condition for the pointG being biologically admissible is a<1.297; therefore a3 is always positive. If a<0.717455, then a2>0 and a1a2-a3>0; then, according to the Routh-Hurwitz Criterion, G is a stable focus. On the other hand, if a>0.717455, then a1a2-a3<0 and a3>0. Hence, from item (b) of Lemma 1 it follows that (13) has two roots with positive real part and a negative real root. Thus G is an unstable focus, remaining unstable for a2>0 or a2<0.
Ifa=0.717455, then a1,a2,a3>0 and a1a2-a3=0; 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 G 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 G: (26)l1=-1.854776782·108.
To finish the proof, the transversality condition remains to be checked. Using software MAPLE, we calculate the derivative of the Jacobian matrixA(a) with respect to the parameter a applying the critical value a=0.717455 from which we find (27)dAda(0.717455)=-0.668742-0.1325030000000.
Solving the equation(28)dγda(0.717455)=Rep,dAda(0.717455)q,we obtain (29)dγda(0.717455)=9.18·107≠0.
This ends the proof of Theorem2.The numerical simulations presented in Figures1 and 2 illustrate the results stated in Theorem 2.Figure 1
Phase portrait of system (1) near the equilibrium point G, with initial conditions: (0.1325,0.66,0.19). (a) G is a stable focus for a=0.6. Integration time: [0,300]. (b) G is a weakly stable focus for a=0.71. Integration time: [0,3000].
(a)
(b)Figure 2
Phase portrait of system (1) near the equilibrium G, with initial conditions: (0.1325,0.66,0.19). (a) G is an unstable focus for a=0.8. Integration time: [50,500]. (b) Limit cycle created with the Hopf bifurcation that occurs when a passes through the critical value a=0.717455. Integration time: [450,500].
(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 parametera moves away from the critical value a=0.717455, 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 a tends to the limit value a=1.


## 4. Period-Doubling Bifurcations and Chaotic Dynamics

Fora close to and greater than the critical value a=0.717455, for which the Hopf bifurcation at the equilibrium G occurs, we have a stable limit cycle of period T (see Figure 3); increasing the value of a, we numerically found that this limit cycle splits into a period 2T periodic orbit (see Figure 4). When the value of a increases a bit more, a 4T periodic orbit is observed (see Figure 5); increasing further the value of a, a series of n duplications and a limit cycle of period 2nT are observed, until a chaotic attractor is formed for a≅1 (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).Figure 3
Periodic orbit of periodT created by the Hopf bifurcation that occurs at the equilibrium point G. (a) Limit cycle for system (1), for a=0.9. Integration time: [1500,1600]. Initial conditions: (0.1325,0.6687,0.19), (0.13,0.6,0.2). (b) Periodic behaviour of the coordinate x(t) of the solution in (a). Integration time: [1500,1600]. Initial conditions: (0.1325,0.6687,0.2).
(a)
(b)Figure 4
Period-doubling bifurcation occurred with the periodic solution created in the Hopf bifurcation at the equilibriumG. (a) Periodic orbit of period 2T, for a=0.93. Integration time: [1500,1600]. Initial conditions: (0.1325,0.6687,0.19), (0.13,0.6,0.2). (b) Periodic behaviour of the coordinate x(t) of the solution in (a). Integration time: [1500,1600]. Initial conditions: (0.1325,0.6687,0.2).
(a)
(b)Figure 5
(a) Periodic orbit of period4T, for a=0.948. Integration time: [1500,1600]. Initial conditions: (0.1325,0.6687,0.19), (0.13,0.6,0.2). (b) Periodic behaviour of the coordinate x(t) of the solution in (a). Integration time: [1500,1600]. Initial conditions: (0.1325,0.6687,0.2).
(a)
(b)Figure 6
Chaotic attractor created through the period-doubling bifurcations. Initial conditions:(0.1325,0.6687,0.19), (0.13,0.6,0.2). (a) Chaotic attractor of system (1), with a=0.99. Integration time: [0,1850]. (b) Sensitive dependence on initial conditions. Two solutions very close diverge exponentially in time, for a=0.99. Integration time: [800,1550].
(a)
(b)Figure6(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 Figure7.Figure 7
Bifurcation diagram for system (1). From left to right we see bifurcations of period-2k: as the value of the parameter a increases, successive bifurcations occur, resulting in chaos.


### 4.1. The Elimination of Tumour and Immune System Cells after the Chaotic Dynamics

The chaotic dynamics occurs when the value ofa is approximately 1. In the adimensional system (1) the parameter a, 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 by(30)dTdt~=r1T1-Tk1-a~TH-b~TE,dHdt~=c~H1-Hk2-d~TH,dEdt~=e~TET+f~-g~TE-h~E,in the following way: (31)a=a~k2r1,from which we have (32)r1=a~k2a.Thus, for small values of parameter a we have large values of r1, since k2 is the largest parameter of the original system, because it represents the carrying capacity of healthy cells. Then the value of r1, which represents the growth rate of tumour cells, is large when a≤1. This explains the fact that the tumour-free equilibrium point B=(0,1,0) is unstable, since, as the growth rate of the tumour cells is large, the tumour cells are increasing when a≤1. 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 a≤1. 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 a=1. 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, fora>1, the value of r1 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 Figure8 we present some numerical simulations of the solutions of system (1) with a>1, showing the stability of the equilibrium point B=(0,1,0), illustrating the fact that, for the set of parameters considered here, only healthy cells survive after the battle between tumour cells and immune cells.Figure 8
Phase portrait of system (1) for a=1.2 showing the asymptotic stability of the point B=(0,1,0). Integration time: [0,500]. Initial conditions: (0.01,0.01,0.01), (0.01,0.9,0.01), (1.5,0.01,0.01), (0.01,0.2,0.01), (0.01,1.5,0.01),(0.9,0.01,0.01),(0.01,0.1,0.01),(0.01,0.01,0.1).


## 4.1. The Elimination of Tumour and Immune System Cells after the Chaotic Dynamics

The chaotic dynamics occurs when the value ofa is approximately 1. In the adimensional system (1) the parameter a, 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 by(30)dTdt~=r1T1-Tk1-a~TH-b~TE,dHdt~=c~H1-Hk2-d~TH,dEdt~=e~TET+f~-g~TE-h~E,in the following way: (31)a=a~k2r1,from which we have (32)r1=a~k2a.Thus, for small values of parameter a we have large values of r1, since k2 is the largest parameter of the original system, because it represents the carrying capacity of healthy cells. Then the value of r1, which represents the growth rate of tumour cells, is large when a≤1. This explains the fact that the tumour-free equilibrium point B=(0,1,0) is unstable, since, as the growth rate of the tumour cells is large, the tumour cells are increasing when a≤1. 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 a≤1. 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 a=1. 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, fora>1, the value of r1 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 Figure8 we present some numerical simulations of the solutions of system (1) with a>1, showing the stability of the equilibrium point B=(0,1,0), illustrating the fact that, for the set of parameters considered here, only healthy cells survive after the battle between tumour cells and immune cells.Figure 8
Phase portrait of system (1) for a=1.2 showing the asymptotic stability of the point B=(0,1,0). Integration time: [0,500]. Initial conditions: (0.01,0.01,0.01), (0.01,0.9,0.01), (1.5,0.01,0.01), (0.01,0.2,0.01), (0.01,1.5,0.01),(0.9,0.01,0.01),(0.01,0.1,0.01),(0.01,0.01,0.1).


## 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 parametera, then we observe that after the chaotic dynamics (i.e., for a>1) 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.



---

*Source: 354918-2015-03-22.xml*