Periodic Orbits for a Three-Dimensional Biological Differential Systems

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.

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) where 𝑢_𝑔, 𝑢_𝑐, 𝑎_𝑐, 𝑎_𝑔, 𝑒, 𝑚∈ℝ are positive parameters and where (𝑍,𝐶,𝐺,𝑃)⊂ℝ⁴₊.

In the previous system 𝑍 represents the density of population of a predator (Daphnia or Zooplankton) that predates the preys (Phytoplankton) 𝐶 and 𝐺 that are of different size, in particular 𝐶 being of a smaller size than 𝐺. The variable 𝑃 represents the amount of resources (Phosphorous) for the preys 𝐶 and 𝐺. The system admits a constant of motion 𝐼(𝑡)=𝑍+𝐶+𝐺+𝑃, and then we can reduce the degree of freedom of the problem by fixing a value of the first integral, 𝐼(𝑡)=𝐼₀. The system can be rewritten in the following way: 𝑍^′=𝑍−𝑒+𝑢_𝑐𝐶+𝑢_𝑔𝐺,𝐶^′=𝐶𝑎_𝑐𝐼₀−𝑎_𝑐+𝑢_𝑐𝑍−𝑎_𝑐𝐶−𝑎_𝑐𝐺𝐺−𝑚𝐶𝑍,_′=𝐺𝑎_𝑔𝐼₀−𝑎_𝑔+𝑢_𝑔𝑍−𝑎_𝑔𝐶−𝑎_𝑔𝐺+𝑚𝐶𝑍. (2) 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., 𝑚=0) of the terms that describe indirect effects, the system (2) does not admits coexistence of the three populations, while for 𝑚>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 𝑚 being a bifurcation parameter for coexistence. In particular under the hypotheses 𝐼₀𝑒>max𝑢_𝑐,𝑒𝑢_𝑔,𝑀=𝑎_𝑐𝑢_𝑔−𝑎_𝑔𝑢_𝑐−𝑎_𝑔𝑚>0, (3)(4) the system is uniformly persistent. Under the previous hypotheses the system admits an attractor 𝒜 that lies at positive distance from 𝜕ℝ³₊ and inside the set Δ=(𝑍,𝐶,𝐺)∈ℝ³₊∶𝑍+𝐶+𝐺≤𝐼₀. (5) In particular the boundary fixed point 𝑃₀=(0,0,0),𝑃₂=0,𝐼₀,0,𝑃₃=0,0,𝐼₀,𝑃₄=𝑎_𝑔𝑎_𝑔+𝑢_𝑔𝐼₀−𝑒𝑢_𝑔𝑒,0,𝑢_𝑔 (6) and the segment of fixed point 𝑆∶=(𝑍,𝐶,𝐺)∈ℝ³∶𝑍=0,𝐶+𝐺=𝐼₀ (7) are all instable and their stable and central manifolds are all on 𝜕𝑅³₊. Under the hypotheses (3)-(4) the set Δ always contains an interior fixed point 𝑄=(𝑍_𝑄,𝐶_𝑄,𝐺_𝑄) as follows: 𝐺_𝑄=𝑒/𝑢_𝑐𝑀/𝑎_𝑐𝑚+𝑢_𝑔/𝑢_𝑐,𝐶_𝑄=𝑀𝑎_𝑐𝑚𝑒/𝑢_𝑐𝑀/𝑎_𝑐𝑚+𝑢_𝑔/𝑢_𝑐,𝑍_𝑄=𝑎_𝑐𝑚_𝑐𝐼₀−𝑀𝑚𝑎_𝑐+1𝑒/𝑢_𝑐𝑀/𝑎_𝑐𝑚+𝑢_𝑔/𝑢_𝑐, (8) whose stability character has not been investigated in [14]. The matrix 𝐴(𝑄) of the linearized system at 𝑄 is 0𝑢_𝑐𝑍_𝑄𝑢_𝑔𝑍_𝑄−𝑚_𝑐𝑀𝑚𝑎_𝑐𝑀𝐷−𝑚𝑀𝐷−𝑚𝐷−𝑎_𝑔𝑚_𝑐𝑎_𝑐𝐷𝑚𝑚_𝑐𝑎_𝑐𝐼₀−𝑚_𝑔𝑚𝐷−𝑀𝑎_𝑐𝑍_𝑄−𝑎_𝑔𝐷, (9) where 𝑚_𝑐=𝑚+𝑎_𝑐+𝑢_𝑐 and 𝑚_𝑔=𝑚+𝑎_𝑔+𝑢_𝑔 and 𝐷=𝐺_𝑄. By numerical experiments (see [14]) it was pointed out that the attractor 𝒜 should take the following forms: (1)the interior fixed point 𝑄; (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 𝐼₀. In general cases we can use the following remark to facilitate our analysis.

Remark 1. The trace of the matrix 𝐴(𝑄) satisfies the following inequality (𝐴(𝑄𝑀Tr))=−𝑎_𝑐𝑍_𝑄−𝐺_𝑄𝑎_𝑔+𝑀𝑚<0; (10) 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 Section 2 we present general condition to have Hopf bifurcation, while in Section 3 we give an example fixing all the parameters except for 𝐼₀. 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 matrix 𝐴(𝑄) as follows: 𝑥³+𝛼𝑥²+𝛽𝑥+𝛾=0, (11) where [𝐴(𝑄=𝑀𝛼=−Tr)]𝑚_𝑐𝐼₀𝑀−𝐷𝑎_𝑐𝑚−𝑚+1^𝑐𝑎_𝑔𝑀−𝑚_𝑐𝑚,(𝐴(𝑄=𝛽=𝑍))𝑀𝐷𝑚𝐼₀−𝑀𝑎_𝑐𝑚𝐷×𝑢+1_𝑐+𝑀𝑚_𝑐+𝑢_𝑔𝑎_𝑔𝑚𝑀+𝑚𝑎_𝑐𝑚_𝑐,𝛾=−Det[𝐴(𝑄)]=𝑒𝑀𝑎_𝑐𝑚_𝑐𝐼₀−𝑀𝑎_𝑐𝑚+1𝐷². (12) A necessary condition for Hopf bifurcation is that each of the previous coefficients has to be positive. It easy to verify that 𝛼, 𝛽, 𝛾>0 if 𝐼₀>𝑀𝑎_𝑐𝑚+1𝐷, (13) and the last inequality follows from hypothesis (3).

From Hurwitz-Routh criterion we need to study the sign of the term 𝑀𝛼𝛽−𝛾=𝑚_𝑐𝐵𝐼₀−𝑀𝑎_𝑐𝑚+1𝐷. (14) The sign of the previous term depends only on the factor 𝐵; in details we have 𝐵=𝐼₀𝑚−𝐷𝐴−^𝑐𝑎_𝑔𝑀−𝑚_𝑐𝑚𝐹−𝑒𝑎_𝑐𝐼₀=−𝐴𝐷𝐹−𝑒𝑎_𝑐𝐼₀−𝐼^∗₀, (15) where 𝑀𝐴=𝑎_𝑐𝑚+1,𝐹=𝑀𝐷𝑚𝑢_𝑐+𝑀𝑚_𝑐+𝑢_𝑔𝑎_𝑔𝑚𝑀+𝑚𝑎_𝑐𝑚_𝑐,𝐼^∗₀𝐹𝑚=𝐴𝐷−𝐷^𝑐𝑎_𝑔/𝑀+𝑚_𝑐/𝑚𝐹−𝑒𝑎_𝑐. (16) From the hypothesis (3) we already have 𝐼₀>𝐴𝐷; then if 𝐹−𝑒𝑎_𝑐>0, (17) we have that 𝛼𝛽−𝛾>0 for any 𝐼₀>𝐴𝐷, and this implies that the real parts of the eigenvalues are all negative and so the fixed point 𝑄 is asymptotically stable. Moreover, all the persistence solutions converge to the fixed point 𝑄.

On the contrary if 𝐹−𝑒𝑎_𝑐<0, (18) we have that 𝛼𝛽−𝛾>0 for any 𝐼₀∈(1,𝐼^∗₀). In this case the eigenvalues have all negative real parts and the fixed point 𝑄 is asymptotically stable. If 𝐼₀>𝐼^∗₀, then 𝛼𝛽−𝛾 becomes negative and the fixed point 𝑄 loses stability. Moreover a periodic orbit appears thanks to Hopf bifurcation theorem (with 𝐼^∗₀ 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 𝐴(𝑄) has one negative and two purely imaginary eigenvalues if and only if 𝐼₀=𝐼^∗₀ 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 Section 2. We prove the existence of a limit cycle by considering 𝐼₀ as a bifurcation parameter and fixing the values of the others parameters.

We set 𝑢_𝑐=𝑢_𝑔=𝑒=𝑎_𝑔=𝑚=1,𝑎_𝑐=3. (19) We will choose 𝐼₀ in order to satisfy hypothesis (3) while hypothesis (4) is already satisfied as 𝐼₀>1,𝑀=1>0. (20) The fixed point 𝑄 has the following coordinates: 𝑍_𝑄,𝐶_𝑄,𝐺_𝑄=35𝐼₀,1−14,34, (21) and it then exists for 𝐼₀>1 which is verified by hypothesis (3). The Jacobian matrix at point 𝑄 is the following: 𝐴(𝑄)=035𝐼₀3−15𝐼₀−5−14−34−34−5435𝐼₀−94−15𝐼₀+114. (22) The characteristic polynomial of the matrix 𝐴(𝑄) takes the form 𝜆³+𝛼𝜆²+𝛽𝜆+𝛾=0, (23) where 1𝛼=5𝐼₀+132>0,𝛽=21𝐼10⁰3−1>0,𝛾=5𝐼₀−1²>0, (24) 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 9𝛼𝛽−𝛾=−𝐼50⁰−1𝐼₀−372, (25) which is positive for 𝐼₀∈(1,37/2) and negative for 𝐼₀>37/2. Then Hopf bifurcation occurs at 𝐼^𝐻₀=372, (26) and this gives the existence of a family of periodic orbits.

In particular if 𝐼₀ satisfies hypothesis (3) and 𝐼₀<𝐼^𝐻₀, we have three eigenvalues with negative real part and then the fixed point is stable and 𝒜=𝑄 (see Figure 1). In the case in which 𝐼₀>𝐼^𝐻₀ 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 𝐼₀=𝐼^𝐻₀ the eigenvalues are 𝜆₁=−5,𝜆_2,37√=±𝑖32. (27) The period of the limit cycle is 𝑇=2𝜋/𝜔, where 𝜔=7√3/2 as represented in Figure 3.

Figure 1 

The solution of the system with 𝐼0<𝐼𝐻0. The attractor is the fixed point 𝑄.

Figure 2 

The solution of the system with 𝐼0>𝐼𝐻0. The attractor is a periodic orbit.

(a)

(b)

(c)

Figure 3 

The period of the limit cycle for 𝐼0=𝐼𝐻0, 𝑇=4𝜋/7√3≃1.03.

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: 𝑍^′(),𝐶=𝑍−1+𝐶+𝐺_′=𝐶3𝐼₀,𝐺−5𝑍−3𝐶−3𝐺_′=𝐺𝐼₀−2𝑍−𝐶−𝐺+𝐶𝑍. (28) We consider the following change of coordinates: 𝑍=𝑥₁+35𝐼^𝐻₀−1,𝐶=𝑥₂+14,𝐺=𝑥₃+34, (29) and we rewrite the system 𝑥^′₁=𝑥₁+212𝑥₂+𝑥₃,𝑥^′₂=−𝑥₂+145𝑥₁+3𝑥₂+3𝑥₃,𝑥^′₃=𝑥₃+34−72−2𝑥₁−𝑥₂−𝑥₃+𝑥₂+14𝑥₁+212. (30) in matrix form 𝑥^′₁𝑥^′₂𝑥^′₃(𝑄)𝑥=𝐴𝑥¹₂𝑥₃+𝑥₁𝑥₂+𝑥₃−𝑥₂5𝑥₁+3𝑥₂+3𝑥₃𝑥₁𝑥₂−2𝑥₁𝑥₃−𝑥₂𝑥₃−𝑥²₃, (31) where the matrix 𝐴(𝑄) is calculated at 𝐼₀=𝐼^𝐻₀. Let 𝑣₁, 𝑣₂, 𝑣₃ be the eigenvectors of 𝐴(𝑄). Consider the following transformations: 𝑣^′₁=𝑣₁=−84355,−,𝑣11,1_′2=𝑣₂+𝑣₃2=12,12,𝑣,1_′3=𝑣₂−𝑣₃=−32𝑖2√3,12√3,0 (32) and the change of variable 𝑋₁𝑋₂𝑋₃𝑥=𝑃𝑥¹₂𝑥₃, (33) where 𝑃 is the matrix with vectors 𝑣^′₁, 𝑣^′₂, 𝑣^′₃ as columns. Then the system becomes 𝑋̇₁𝑋̇₂𝑋̇₃=𝑃⁻¹𝐴(𝑄)𝑃𝑋₁𝑋₂𝑋₃+𝑃⁻¹𝑥₁𝑥₂+𝑥₃−𝑥₂5𝑥₁+3𝑥₂+3𝑥₃𝑥₁𝑥₂−2𝑥₁𝑥₃−𝑥₂𝑥₃−𝑥²₃; (34) that is 𝑋̇₁𝑋̇₂𝑋̇₃=7√−50000327√0−320𝑋₁𝑋₂𝑋₃+𝑓₁𝑓₂𝑓₃, (35) where 𝑓₁=55𝑋₁+9𝑋₂−5𝑋₃×√244036−55+853𝑋₁−495+223√3𝑋₂+3−73+23√3𝑋₃,𝑓₂=55𝑋₁+495𝑋₂+219𝑋₃×5√24403622+513𝑋₁+990−669√3𝑋₂+−2032+207√3𝑋₃,𝑓₃=1√244036−25988+1873𝑋²₁+3−131404+36795√3𝑋²₂+260021−102821√3𝑋₂𝑋₃+315314+17687√3𝑋²₃−10𝑋₁−5928+2981√3𝑋₂+3211−6157√3𝑋₃. (36) We consider a 2-dimensional center manifold of the type 𝑋₁=ℎ𝑋₂,𝑋₃=𝑎𝑋²₂+𝑏𝑋₂𝑋₃+𝑐𝑋²₃. (37) Using the first equation of the system, we get 𝑋̇₁=−5𝑋²₂√𝑎+99495+2233244036−5𝑋²₃√𝑐+165−73+233244036+5𝑋₂𝑋₃√−𝑏+1118513+504,244036 (38) while using (37) and the equations of the systems, we get 𝑋̇₁=7√3−𝑏2𝑋²₂+𝑏2𝑋²₃+(𝑎−𝑐)𝑋₂𝑋₃. (39) Equaling the coefficients, we get √−5𝑎+99495+2233√244036=−73𝑏2,√−5𝑐+165−73+233√244036=73𝑏2,5√−𝑏+1118513+504√244036=73(𝑎−𝑐), (40) with solutions √𝑎=−9958035+476063,√83948384𝑏=49535345+1386√3,√41974192𝑐=−33211175+1268783.83948384 (41) The system restricted on the central manifold takes the form 𝑋̇₂𝑋̇₃=02133√160−38152133√160𝑋38150𝑋²₃+𝑔₂𝑋₂,𝑋₃𝑔₃𝑋₂,𝑋₃, (42) where 𝑔₂(𝑋₂,𝑋₃)=𝑓₂(ℎ(𝑋₂,𝑋₃),𝑋₂,𝑋₃), 𝑔₃(𝑋₂,𝑋₃)=𝑓₃(ℎ(𝑋₂,𝑋₃),𝑋₂,𝑋₃). As remarked in the previous section the periodic orbit should persist and should be stable for 𝐼₀≥𝐼^𝐻₀.

The first Lyapunov coefficient [19, page 99] can be calculated as 𝑙₁(0)=−32340518937237+21544969548257√3501917482026688≃−0.0099149, (43) 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 𝑚 is replaced by the time-dependent term: 𝑚|sin𝑡|. (44) 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 𝑍^′=𝑍−𝑒+𝑢_𝑐𝐶+𝑢_𝑔𝐺,𝐶^′=𝐶𝑎_𝑐𝐼₀−𝑎_𝑐+𝑢_𝑐𝑍−𝑎_𝑐𝐶−𝑎_𝑐𝐺||𝐺−𝑚sin𝑡𝐶𝑍,_′=𝐺𝑎_𝑔𝐼₀−𝑎_𝑔+𝑢_𝑔𝑍−𝑎_𝑔𝐶−𝑎_𝑔𝐺+𝑚|sin𝑡|𝐶𝑍. (45) 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 but 𝐼₀ as follows: 𝑢_𝑔=𝑢_𝑐=𝑚=𝑒=𝑎_𝑔=1,𝑎_𝑐=3. (46) In Figure 4 we represent the solution for 𝐼₀=7, 14, 28 while in Figure 5 for 𝐼₀=35, 82, 177. We note that as 𝐼₀ increases, the period of the solution increases, and in Figure 6 we represent the solutions 𝑍, 𝐶, 𝐺 as function of time for 𝐼₀=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 𝑍(𝑡) for 𝐼₀=7,35,177.

(a)

(b)

(c)

Figure 4 

The solutions for 𝐼0=7, 14, 28, respectively.

(a)

(b)

(c)

Figure 5 

The solution for 𝐼0=35, 82, 177, respectively.

(a)

(b)

(c)

Figure 6 

The solutions 𝑍, 𝐶, 𝐺 as function of time 𝐼0=177.

(a)

(b)

(c)

Figure 7 

The period of the function 𝑍 for 𝐼0=7, 35, 177.

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) then the indirect effects terms disappear and so indirect effects do not change the total amount of the preys 𝐶+𝐺 but only the proportion between them. A simple generalization of the model consists in considering different values of the parameter 𝑚, that is, 𝑚₁ and 𝑚₂, respectively, in the second and the third equations of the system 𝑍^′=𝑍−𝑒+𝑢_𝑐𝐶+𝑢_𝑔𝐺,𝐶^′=𝐶𝑎_𝑐𝐼₀−𝑎_𝑐+𝑢_𝑐𝑍−𝑎_𝑐𝐶−𝑎_𝑐𝐺−𝑚₁𝐺𝐶𝑍,_′=𝐺𝑎_𝑔𝐼₀−𝑎_𝑔+𝑢_𝑔𝑍−𝑎_𝑔𝐶−𝑎_𝑔𝐺+𝑚₂𝐶𝑍. (48) The dynamical behavior of the system (48) is very similar to the case 𝑚₁=𝑚₂ (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: 𝐼₀>𝑒𝑢_𝑐,𝑒𝑢_𝑔,𝑀₁=𝑎_𝑐𝑢_𝑔−𝑎_𝑔𝑢_𝑐−𝑎_𝑔𝑚₁>0. (49)

We note that the previous conditions do not depend on the parameter 𝑚₂. The differences between the two models appear in the coordinate of the fixed point 𝑄𝑄𝑄𝑄=(𝑍,𝐶,𝐺) as follows: 𝐺𝑄=𝑒/𝑢_𝑐𝑀₁/𝑎_𝑐𝑚₂+𝑢_𝑔/𝑢_𝑐,𝐶𝑄=𝑀₁𝑎_𝑐𝑚₂𝑒/𝑢_𝑐𝑀₁/𝑎_𝑐𝑚₂+𝑢_𝑔/𝑢_𝑐,𝑍𝑄=𝑎_𝑐𝑚¹_𝑐𝐼₀−𝑀₁𝑚₂𝑎_𝑐+1𝑒/𝑢_𝑐𝑀₁/𝑎_𝑐𝑚₂+𝑢_𝑔/𝑢_𝑐, (50) where 𝑚¹_𝑐=𝑎_𝑐+𝑢_𝑐+𝑚₁. In particular, as for the original system, if the hypotheses of the previous theorem are satisfied, the interior fixed point 𝑄 always exists.

The parameter 𝑚₂ plays a role only in the stability of the point 𝑄 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 𝐴(𝑄) has one negative and two purely imaginary eigenvalues if and only if 𝐼₀=𝐼^∗₀̃ and 𝐹−𝑒𝑎_𝑐<0, (51) where 𝑀𝐹=¹𝑚₂𝑒/𝑢_𝑐𝑀₁/𝑎_𝑐𝑚₂+𝑢_𝑔/𝑢_𝑐×𝑢_𝑐+𝑀₁𝑚¹_𝑐+𝑢_𝑔𝑎_𝑔𝑚₂𝑀₁+𝑚₂𝑎_𝑐𝑚¹_𝑐,𝑄,𝑀𝐷=𝐺𝐴=¹𝑎_𝑐𝑚₂,𝐼̃+1_∗0=𝐴𝐷𝐹𝑚𝐷−^1𝑐𝑎_𝑔/𝑀₁+𝑚¹_𝑐/𝑚₂𝐹−𝑒𝑎_𝑐. (52)

We note that inequality (51) depends on 𝑚₁ and 𝑚₂.

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 𝑚₁|sin𝑡| and 𝑚₂|sin𝑡|, respectively, in the second and third equations of the system. If 𝑚₁>𝑚₂, the numerical experiments showed a behavior of the system similar to that described in the previous section (i.e., the case 𝑚₁=𝑚₂). If 𝑚₂>𝑚₁, we have that 𝐶̇+𝐺̇=𝑎_𝑐𝐶+𝑎_𝐺𝐺𝐼₀𝑢−𝑍−𝐶−𝐺−𝑍_𝑐𝐶+𝑢_𝑔𝐺+𝑚₂−𝑚₁|sin𝑡|, (53) 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: 𝑎_𝑐=5,𝑎_𝑔=110,𝑚₁=45,𝑚₂=32𝑢_𝑐=810,𝑢_𝑔=2𝑍(0)(0)(0)=210,𝑒=1,=𝐶=𝐺3. (54) For 𝐼₀≥1.3 the attractor is a periodic orbit, and in Figure 8 we represent the case 𝐼₀=30. The most interesting case is when 𝐼₀=40 where a new structure appears; see Figures 9 and 10.

Figure 8 

The attractor for 𝐼0=30.

Figure 9 

The attractor for 𝐼0=40.

(a)

(b)

(c)

Figure 10 

The solutions 𝑍, 𝐶, 𝐺 as function of time for 𝐼0=40.

For higher value of 𝐼₀ 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 𝑍, 𝐶, 𝐺 for 𝐼₀=50.

Figure 11 

The attractor for 𝐼0=50.

(a)

(b)

(c)

Figure 12 

The solutions 𝑍, 𝐶, 𝐺 as function of time for 𝐼0=50.

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

Acknowledgment

This work was partially supported by the research group EMAP (Grupo de investigaciòn en matemàatica y estadistica aplicada) of Pontificia Universidad Javeriana of Cali, Colombia.