# Periodic Orbits for a Three-Dimensional Biological Differential Systems

Authors: Renato Colucci; Daniel Nuñez
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We study the existence of periodic orbit for a differential system describing the effects of indirect predation over two preys. Besides discussing a generalized version of the model, we present some remarks and numerical experiments for the nonautonomous version of the two models.


---

## Body

## 1. Introduction

The role of indirect effects in biology has been extensively studied in the last decades (see [1–7]); for example, in the case of predation (see [8]), the predator can alter the morphology (see [9]) or the behavior of the preys. The preys, in order to reduce the possibility of contacts with the predators, could modify their normal conduct by reducing their activity or by hiding themselves for long time. There are many types of indirect effects (see [3] for a detailed discussion). Another interesting case is the refuge indirect effect (see [10] for an example); anyway it is of great interest trying to describe the indirect interactions in population dynamics. In [11] a model was proposed including indirect effects, modeling the effects of predator Daphnia over two groups of Phytoplankton of different morphology (see [12, 13]), having Phosphorous as resource (see [13] or [7]). The system takes the following form
(1)Z′=Z(-e+ucC+ugG),C′=C(acP-ucZ)-mCZ,G′=G(agP-ugZ)+mCZ,P′=P(-agG-acC)+eZ,
where ug, uc, ac, ag, e, m∈ℝ are positive parameters and where (Z,C,G,P)⊂ℝ+4.In the previous systemZ represents the density of population of a predator (Daphnia or Zooplankton) that predates the preys (Phytoplankton) C and G that are of different size, in particular C being of a smaller size than G. The variable P represents the amount of resources (Phosphorous) for the preys C and G. The system admits a constant of motion I(t)=Z+C+G+P, and then we can reduce the degree of freedom of the problem by fixing a value of the first integral, I(t)=I0. The system can be rewritten in the following way:
(2)Z′=Z(-e+ucC+ugG),C′=C[acI0-(ac+uc)Z-acC-acG]-mCZ,G′=G[agI0-(ag+ug)Z-agC-agG]+mCZ.
The dynamics of system (2) have been studied in [14] by using Persistence Theory (see [15, 16]) and in particular an acyclicity approach. It has been shown that in absence (i.e., m=0) of the terms that describe indirect effects, the system (2) does not admits coexistence of the three populations, while for m>0 and for an appropriate choice of the parameters, the system admits coexistence. The paper [14] suggests the importance of indirect effects in describing cases of coexistence with the parameter m being a bifurcation parameter for coexistence. In particular under the hypotheses
(3)I0>max{euc,eug},(4)M=acug-aguc-agm>0,
the system is uniformly persistent. Under the previous hypotheses the system admits an attractor 𝒜 that lies at positive distance from ∂ℝ+3 and inside the set
(5)Δ={(Z,C,G)∈ℝ+3:Z+C+G≤I0}.
In particular the boundary fixed point
(6)P0=(0,0,0),P2=(0,I0,0),P3=(0,0,I0),P4=(agag+ug[I0-eug],0,eug)
and the segment of fixed point
(7)S:={(Z,C,G)∈ℝ3:Z=0,C+G=I0}
are all instable and their stable and central manifolds are all on ∂R+3. Under the hypotheses (3)-(4) the set Δ always contains an interior fixed point Q=(ZQ,CQ,GQ) as follows:
(8)GQ={e/ucM/acm+ug/uc},CQ=Macm{e/ucM/acm+ug/uc},ZQ=acmc{I0-(Mmac+1)(e/ucM/acm+ug/uc)},
whose stability character has not been investigated in [14]. The matrix A(Q) of the linearized system at Q is
(9)(0ucZQugZQ-mcMmacD-MmD-MmD-agmcacDmmcac{I0-mgmD}-MacZQ-agD),
where mc=m+ac+uc and mg=m+ag+ug and D=GQ. By numerical experiments (see [14]) it was pointed out that the attractor 𝒜 should take the following forms: (1)
the interior fixed pointQ;(2)
a periodic orbit;(3)
a bidimensional (possibly chaotic) set.In this work we are concerned with the proof of existence of periodic orbits (see [17, 18] for examples in the classical case), and we will use Hopf bifurcation theory in order to prove existence of periodic orbits. Since the system depends on 7 parameters, the conditions to have Hopf bifurcation are quite complicated and for sake of simplicity we also present a particular case in which all the parameters are fixed but I0. In general cases we can use the following remark to facilitate our analysis.Remark 1.
The trace of the matrixA(Q) satisfies the following inequality
(10)Tr(A(Q))=-MacZQ-GQ(ag+Mm)<0;
then at least one eigenvalue has negative real part and as a consequence fold-Hopf bifurcation is not possible.The rest of the paper is organized as follows. In Section2 we present general condition to have Hopf bifurcation, while in Section 3 we give an example fixing all the parameters except for I0. In Section 4 we consider a nonautonomous version of the model and present several numerical experiments in order to study the structure of the attractor. In the last case we discuss a generalized version of the model, both autonomous and nonautonomous, and present some numerical experiments.


## 2. Condition for Hopf Bifurcation

We consider the characteristic polynomial of the matrixA(Q) as follows:
(11)x3+αx2+βx+γ=0,
where
(12)α=-Tr[A(Q)]=Mmc{I0-D[(Macm+1)-mcagM-mcm]},β=Z(A(Q))=MDm{I0-(Macm+1)D}×{uc+Mmc+ugagmM+macmc},γ=-Det[A(Q)]=eMacmc{I0-(Macm+1)D}2.
A necessary condition for Hopf bifurcation is that each of the previous coefficients has to be positive. It easy to verify that α, β, γ>0 if
(13)I0>(Macm+1)D,
and the last inequality follows from hypothesis (3).From Hurwitz-Routh criterion we need to study the sign of the term(14)αβ-γ=MmcB{I0-(Macm+1)D}.
The sign of the previous term depends only on the factor B; in details we have
(15)B={I0-D[A-mcagM-mcm]}F-eac{I0-AD}=(F-eac)(I0-I0*),
where
(16)A=(Macm+1),F=MDm{uc+Mmc+ugagmM+macmc},I0*=AD-DF(mcag/M+mc/m)F-eac.
From the hypothesis (3) we already have I0>AD; then if
(17)F-eac>0,
we have that αβ-γ>0 for any I0>AD, and this implies that the real parts of the eigenvalues are all negative and so the fixed point Q is asymptotically stable. Moreover, all the persistence solutions converge to the fixed point Q.On the contrary if(18)F-eac<0,
we have that αβ-γ>0 for any I0∈(1,I0*). In this case the eigenvalues have all negative real parts and the fixed point Q is asymptotically stable. If I0>I0*, then αβ-γ becomes negative and the fixed point Q loses stability. Moreover a periodic orbit appears thanks to Hopf bifurcation theorem (with I0* as the bifurcation value). Recalling that all the boundary fixed points are instable and repulsive (see [14]), we can conclude that the attractor consists of the periodic orbit generated by Hopf bifurcation. Since the system is uniformly persistent under hypotheses (3)-(4), we expect that the periodic orbit is persistent. We summarize as follows.Proposition 2.
Suppose that (3) and (4) hold. Then the matrix A(Q) has one negative and two purely imaginary eigenvalues if and only if I0=I0* and condition (18) holds.Thanks to Hopf bifurcation theorem, a periodic orbit appears. In order to illustrate this result, in the next section we will consider an example of existence of limit cycle.


## 3. An Example of Existence of a Limit Cycle

In this section we present an example of the results of Section2. We prove the existence of a limit cycle by considering I0 as a bifurcation parameter and fixing the values of the others parameters.We set(19)uc=ug=e=ag=m=1,ac=3.
We will choose I0 in order to satisfy hypothesis (3) while hypothesis (4) is already satisfied as
(20)I0>1,M=1>0.
The fixed point Q has the following coordinates:
(21)(ZQ,CQ,GQ)={35(I0-1),14,34},
and it then exists for I0>1 which is verified by hypothesis (3). The Jacobian matrix at point Q is the following:
(22)A(Q)=(035(I0-1)35(I0-1)-54-34-34-5435(I0-94)-15(I0+114)).
The characteristic polynomial of the matrix A(Q) takes the form
(23)λ3+αλ2+βλ+γ=0,
where
(24)α=15{I0+132}>0,β=2110{I0-1}>0,γ=35{I0-1}2>0,
and where the three inequalities follow from hypothesis (3). We use Hurwitz-Routh criterion to study the sign of the real parts of the eigenvalues. In particular we are interested in the sign of the term
(25)αβ-γ=-950{I0-1}{I0-372},
which is positive for I0∈(1,37/2) and negative for I0>37/2. Then Hopf bifurcation occurs at
(26)I0H=372,
and this gives the existence of a family of periodic orbits.In particular ifI0 satisfies hypothesis (3) and I0<I0H, we have three eigenvalues with negative real part and then the fixed point is stable and 𝒜=Q (see Figure 1). In the case in which I0>I0H we have one eigenvalue with negative real part and two eigenvalues with positive real part. In this case the attractor is a periodic orbit (see Figure 2). For I0=I0H the eigenvalues are
(27)λ1=-5,λ2,3=±i732.
The period of the limit cycle is T=2π/ω, where ω=73/2 as represented in Figure 3.Figure 1
The solution of the system withI0<I0H. The attractor is the fixed point Q.Figure 2
The solution of the system withI0>I0H. The attractor is a periodic orbit.The period of the limit cycle forI0=I0H, T=4π/73≃1.03.
(a)
(b)
(c)We derive the equation of the central manifold in order to investigate the character of the bifurcation. We rewrite the systems with parameters as in (19) as follows:
(28)Z′=Z(-1+C+G),C′=C(3I0-5Z-3C-3G),G′=G(I0-2Z-C-G)+CZ.
We consider the following change of coordinates:
(29)Z=x1+35(I0H-1),C=x2+14,G=x3+34,
and we rewrite the system
(30)x1′=(x1+212)(x2+x3),x2′=-(x2+14)(5x1+3x2+3x3),x3′=(x3+34)(-72-2x1-x2-x3)+(x2+14)(x1+212).
in matrix form
(31)(x1′x2′x3′)=A(Q)(x1x2x3)+(x1(x2+x3)-x2(5x1+3x2+3x3)(x1x2-2x1x3-x2x3-x32)),
where the matrix A(Q) is calculated at I0=I0H. Let v1, v2, v3 be the eigenvectors of A(Q). Consider the following transformations:
(32)v1′=v1=(-8455,-311,1),v2′=v2+v32=(12,12,1),v3′=v2-v32i=(-323,123,0)
and the change of variable
(33)(X1X2X3)=P(x1x2x3),
where P is the matrix with vectors v1′, v2′, v3′ as columns. Then the system becomes
(34)(X˙1X˙2X˙3)=P-1A(Q)P(X1X2X3)+P-1(x1(x2+x3)-x2(5x1+3x2+3x3)(x1x2-2x1x3-x2x3-x32));
that is
(35)(X˙1X˙2X˙3)=(-500007320-7320)(X1X2X3)+(f1f2f3),
where
(36)f1=55(X1+9X2-5X3)244036×{((-55+853)X1-(495+2233)X2+3(-73+233)X3)},f2=55X1+495X2+219X3244036×{5(22+513)X1+(990-6693)X2+(-2032+2073)X3},f3=1244036{-25(988+1873)X12+3(-131404+367953)X22+2(60021-1028213)X2X3+3(15314+176873)X32-10X1((-5928+29813)X2+(3211-61573)X3)}.
We consider a 2-dimensional center manifold of the type
(37)X1=h(X2,X3)=aX22+bX2X3+cX32.
Using the first equation of the system, we get
(38)X˙1=-5X22(a+99(495+2233)244036)-5X32(c+165(-73+233)244036)+5X2X3(-b+11(18513+504)244036),
while using (37) and the equations of the systems, we get
(39)X˙1=73(-b2X22+b2X32+(a-c)X2X3).
Equaling the coefficients, we get
(40)-5(a+99(495+2233)244036)=-73b2,-5(c+165(-73+233)244036)=73b2,5(-b+11(18513+504)244036)=73(a-c),
with solutions
(41)a=-99(58035+476063)83948384,b=4953(5345+13863)41974192,c=-33(211175+1268783)83948384.
The system restricted on the central manifold takes the form
(42)(X˙2X˙3)=(021331603815-213316038150)(X2X3)+(g2(X2,X3)g3(X2,X3)),
where g2(X2,X3)=f2(h(X2,X3),X2,X3), g3(X2,X3)=f3(h(X2,X3),X2,X3). As remarked in the previous section the periodic orbit should persist and should be stable for I0≥I0H.The first Lyapunov coefficient [19, page 99] can be calculated as
(43)l1(0)=-32340518937237+215449695482573501917482026688≃-0.0099149,
and this confirm that the Hopf bifurcation is supercritical.


## 4. The Nonautonomous Case

Since indirect effects can be of seasonal type (see [14]), it is interesting to study the case in which the constant m is replaced by the time-dependent term:
(44)m|sint|.
We consider values of the parameters for which a periodic orbit exists and make the substitution (44) in the indirect effects terms of the system
(45)Z′=Z(-e+ucC+ugG),C′=C[acI0-(ac+uc)Z-acC-acG]-m|sint|CZ,G′=G[agI0-(ag+ug)Z-agC-agG]+m|sint|CZ.
We consider some numerical experiments which show that not only the periodic orbit persists but also sort of period bifurcation phenomena occur.For simplicity we fix all the parameters butI0 as follows:
(46)ug=uc=m=e=ag=1,ac=3.
In Figure 4 we represent the solution for I0=7, 14, 28 while in Figure 5 for I0=35, 82, 177. We note that as I0 increases, the period of the solution increases, and in Figure 6 we represent the solutions Z, C, G as function of time for I0=177 showing different frequencies of oscillations. It would be interesting to analyze this phenomenon in details together with the possibility to have chaotic motion. In Figure 7 we compare the period and oscillations of the function Z(t) for I0=7,35,177.The solutions forI0=7, 14, 28, respectively.
(a)
(b)
(c)The solution forI0=35, 82, 177, respectively.
(a)
(b)
(c)The solutionsZ, C, G as function of time I0=177.
(a)
(b)
(c)The period of the functionZ for I0=7, 35, 177.
(a)
(b)
(c)


## 5. A Generalized Model

In this section we consider a slighter modification of the model, both in the autonomous and nonautonomous cases. If we sum the second and the third equations of the system, we get(47)C˙+G˙=(acC+aGG)(I0-Z-C-G)-Z(ucC+ugG);
then the indirect effects terms disappear and so indirect effects do not change the total amount of the preys C+G but only the proportion between them. A simple generalization of the model consists in considering different values of the parameter m, that is, m1 and m2, respectively, in the second and the third equations of the system
(48)Z′=Z(-e+ucC+ugG),C′=C[acI0-(ac+uc)Z-acC-acG]-m1CZ,G′=G[agI0-(ag+ug)Z-agC-agG]+m2CZ.
The dynamical behavior of the system (48) is very similar to the case m1=m2 (see [14] for details). In particular it has the same boundary fixed points with the same stability character. Then we can state the following without giving the details of computation.Theorem 3.
The system (48) admits coexistence for any positive solution if the following two conditions are satisfied:
(49)I0>{euc,eug},M1=acug-aguc-agm1>0.We note that the previous conditions do not depend on the parameterm2. The differences between the two models appear in the coordinate of the fixed point Q~=(ZQ~,CQ~,GQ~) as follows:
(50)GQ~={e/ucM1/acm2+ug/uc},CQ~=M1acm2{e/ucM1/acm2+ug/uc},ZQ~=acmc1{I0-(M1m2ac+1)(e/ucM1/acm2+ug/uc)},
where mc1=ac+uc+m1. In particular, as for the original system, if the hypotheses of the previous theorem are satisfied, the interior fixed point Q always exists.The parameterm2 plays a role only in the stability of the point Q~ and as a consequence in the structure of the attractor. By the same method of Section 2 it is possible to derive a condition similar to (18) for Hopf bifurcation.Proposition 4.
Suppose that (49) hold. Then the matrix A(Q~) has one negative and two purely imaginary eigenvalues if and only if I0=I~0* and
(51)F~-eac<0,
where
(52)F~=M1m2(e/ucM1/acm2+ug/uc)×{uc+M1mc1+ugagm2M1+m2acmc1},D~=GQ~,A~=(M1acm2+1),I~0*=A~D~-D~F~(mc1ag/M1+mc1/m2)F~-eac.We note that inequality (51) depends on m1 and m2.From the previous discussion we remark that the behavior of the system (48) is qualitative, the same of that of the system (2).Numerical simulations suggest that the nonautonomous version of system (48) presents peculiar features. We put m1|sint| and m2|sint|, respectively, in the second and third equations of the system. If m1>m2, the numerical experiments showed a behavior of the system similar to that described in the previous section (i.e., the case m1=m2). If m2>m1, we have that
(53)C˙+G˙=(acC+aGG)(I0-Z-C-G)-Z(ucC+ugG)+(m2-m1)|sint|,
then, since the last term is nonnegative, we expect a more complicated behavior of the solutions.We present some simulations to described the dynamics of this case. We consider the following values of the parameters:(54)ac=5,ag=110,m1=45,m2=32uc=810,ug=210,e=1,Z(0)=C(0)=G(0)=23.
For I0≥1.3 the attractor is a periodic orbit, and in Figure 8 we represent the case I0=30. The most interesting case is when I0=40 where a new structure appears; see Figures 9 and 10.Figure 8
The attractor forI0=30.Figure 9
The attractor forI0=40.The solutionsZ, C, G as function of time for I0=40.
(a)
(b)
(c)For higher value ofI0 we obtain an attractor similar to what we have found in the previous section but with a more complicated geometry. In Figures 11 and 12 we represent, respectively, the attractor and the function Z, C, G for I0=50.Figure 11
The attractor forI0=50.The solutionsZ, C, G as function of time for I0=50.
(a)
(b)
(c)The numerical experiments suggest that indirect effects encourage an oscillatory behavior of the system and, as a consequence, they are useful to describe real cases of coexistence (see [20]). It could be interesting to investigate in details the nonautonomous case, in particular the phenomenon of period doubling and the existence of chaotic motion (see [20–22] for examples for the classical case).



---

*Source: 465183-2013-06-25.xml*