1. Introduction
Population dynamics has been extensively studied by researchers in biomathematics, particularly the predator–prey models.
In [
1], the authors consider the following two prey one predator model:
where
represent the population densities of the two preys and of the predator, respectively.
In the previous model, for the interaction between the first prey and the predator, they considered a Holling type II functional response where the handling time of predator for the second prey is also involved, whereas for the interaction between the second prey and the predator, they considered a Lotka–Volterra functional response. It is also assumed that there is no intraspecific interaction in the second prey population and its growth is exponential; as a consequence, there is a huge availability of second prey in the absence of a predator, and there is no searching time for the second prey population. They found necessary and sufficient conditions for existence and stability of the nontrivial equilibrium
(see [
1]).
In order to recover a more complex behavior, we consider a modification of the model that takes the indirect effects of predations into account.
The role played by indirect effects in population dynamics has been investigated in the last several decades (see [
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12]). In the case of predation, it has been pointed out (see [
13]) that predator can alter the morphology (see [
6]) or the behavior of the preys. In particular, the preys, in order to avoid contacts with predators, may reduce their normal activity or may stay hidden most of the time. Many kinds of indirect effects have been described in the literature (see, for example, [
8] for a detailed discussion); an interesting example (see [
14]) is the case of refuge indirect effect.
A model was proposed in order to take into account indirect interactions in a plankton community (see [
15] and the references quoted therein). They analyzed the effects of predator
Daphnia over two groups of phytoplankton of different morphology (see [
7,
9]), having phosphorous as a resource (see [
5,
9]). In this case, the predator prefers to predate the smaller size prey group and the other one take advantages of it.
The model was analytically studied in [
15] using persistence theory (see [
16,
17]) and in [
18] using bifurcation theory. Both studies suggest the importance of indirect effects of predation in order to describe cases of coexistence in real life. In [
19], seasonal indirect effects were considered to show the possibility of chaotic motion, whereas in [
20,
21], the authors considered the stochastic version of the model. Regarding the model (
1), since there is a higher availability of the second prey, it is natural to suppose that the predator prefer to predate the second prey and the first one take advantages of it. The easiest way to model this situation consists of adding the indirect effect term
in the second equation and the term
in the first one, with the parameter
describing the intensity of indirect effects. The systems become
We will consider initial conditions
,
,
and we assume all the parameters are positive and with the following meaning:
r and
k are the intrinsic growth rate and carrying capacity of the first prey, respectively;
and
are the intrinsic growth rate and predation rate of the second prey, respectively;
a is the half saturation value of the predator;
b is the maximum growth rate of the predator;
c is the maximum rate of predation for first prey item;
m is the death rate of the predator in the absence of prey;
is the quotient of the handling time of the predator per second prey item and the handling time of the predator per first prey item;
is the quotient of the capture rate of the second prey and the capture rate of the first prey;
is the efficiency with which the second prey consumed by the predator gets converted into predator biomass (see [
1] for more details).
We note that the system is not of Kolmogorov type; indeed, the first equation cannot be written in the form
Such systems can be regarded as semi-Kolmogorov systems using a terminology introduced in [
19].
For system (
2), we perform the stability analysis of equilibria and we analyze the existence of limit cycles by Hopf bifurcation, as they play an important role in the qualitative theory of differential systems. The study of limit cycles was initiated by Poincaré [
22] and motivated by the famous 16th Hilbert problem [
23,
24,
25] and by the fact that the behavior of many natural phenomena has been modelized by limit cycles, as, for instance, the famous limit cycle of van der Pol [
26].
In the last several decades, the existence of limit cycles for systems with biological meaning, such as the Lotka–Volterra or Kolmogorov system, have been studied through Hopf and zero-Hopf bifurcation (see, for example, [
27,
28] for recent results).
The rest of the paper is organised as follows: in
Section 2, we present a preliminary analysis of the features of the model. In
Section 3 we provide a study of existence and stability of equilibria. In
Section 4, we discuss the problem of persistence of the three species. In
Section 5, we present a case of Hopf bifurcation. Finally,
Section 6 contains some conclusive remarks.
2. Analysis of the System (2) on the Invariant Planes
First, we will show that the dynamics of the system, considering positive initial conditions, is contained in the first octant.
Theorem 1. The set is positively invariant for system (2). Proof. At first, we must note that the planes and are invariant. On the plane , we have ; then, solutions do not leave the positive octant, that is, the set is positively invariant. □
We analyze the dynamics on the boundary of . In order to do that, we first study the dynamics on the coordinate axes and then on the planes and .
The three axes are invariant for the dynamics; in particular, any solution with initial conditions on the
x-axis tends to the equilibrium
, any solution with initial conditions on the
z-axis tends to the equilibrium
, and any solution with initial conditions on the
y-axis verifies that
tends to infinity when
t tends to infinity. See
Figure 1.
Now, we present some considerations about the dynamics on the invariant planes. On the invariant plane,
the system is
and for this system, solutions are unbounded except that on the positive
x-axis. The equilibrium points are
and
. The eigenvalues of
are
r and
, which are both positive, so this equilibrium is an unstable node. The eigenvalues of
are
and
, so the equilibrium point is a saddle.
On the plane
, the system is
and the equilibria are
,
. If
there exists a further equilibrium point with coordinates
where
The eigenvalues of are r and , so it is a saddle. The eigenvalues of are and . The first one is negative and the second changes it sign when . When , it is a stable node, and when , it loses its stability; it becomes a saddle and the third equilibrium appears.
The
x-nullclines are
and
. If
, then
is positive if
and negative if
. The local phase portrait in this case is given in
Figure 2.
When
, the equilibrium
appears, as is it shown in
Figure 3. The eigenvalues of the Jacobian matrix of system (
4) at
are
where
Note that
by the existence conditions of the equilibrium point. These eigenvalues are complex if
, and in that case, they have a positive real part if
in which case,
is unstable. If
the real part of the eigenvalues is negative, so
is asymptotically stable. In the case with
, as the determinant of the Jacobian matrix is positive, it is not possible that the eigenvalues have different sign. Then, if
is positive, both eigenvalues are positive and
is unstable. In other case, from the conditions
and
we would get that
, which is a contradiction. In the same way, we obtain that if
, then both eigenvalues are negative and
is asymptotically stable.
We give here some results about the possible existence of periodic orbits surrounding the equilibrium point in the plane .
The equilibrium point is a Hopf equilibrium if and only if and , it is, when . Note that this occurs only for . In general, when a differential system in has an equilibrium with eigenvalues , it can exhibit a Hopf bifurcation, that is, a local bifurcation in which the equilibrium point loses stability as a pair of complex conjugate eigenvalues of the linearization around the equilibrium point, cross the imaginary axis of the complex plane. To show that this bifurcation takes place, it is necessary to compute the first Lyapunov coefficient of the differential system at the equilibrium. When , the equilibrium is a weak focus of the differential system restricted to the central surface of , associated to the pair of complex eigenvalues, which cross the imaginary axis, and the limit cycle that emerges from is stable. In this case, we say that the Hopf bifurcation is supercritical.
Theorem 2. The equilibrium of system (4) undergoes a supercritical Hopf bifurcation at . For , the system has a unique and stable limit cycle bifurcating from the equilibrium point . Proof. We use the results presented on Chapter 3 of [
29] for computing the first Lyapunov coefficient
at the equilibrium
. At first, to simplify calculation, we introduce in system (
4) a new time variable
by
, obtaining the polynomial system:
This system has the positive equilibrium
which is the same as
with the notation introduced in
Section 2. The Jacobian matrix at this equilibrium is
and it has eigenvalues
, where
We get
for
which is positive because as we have said before, a neccesary condition for Hopf bifurcation is
. Moreover,
Therefore, at
, the equilibrium point
has a pair of pure imaginary eigenvalues
and the system has a Hopf bifurcation. The equilibrium is stable for
and unstable for
. In order to analyze this Hopf bifurcation, we will apply Theorem 3.3 in [
29], so we must prove if the genericity conditions are satisfied. We check that the transversality condition is satisfied as
where
represents the derivative with respect to
m, and the sign is determined because
.
To check the second condition, we must compute the first Lyapunov coefficient. We fix the value
, and then, the equilibrium
has the expression
We translate
to the origin of coordinates obtainig the system
which can be represented as
where
and the multilinear functions
B and
C are given by
We need to find two eigenvectors
of the matrix
A verifying
as for example
Now, we compute
and the first Lyapynov coefficient
which is negative for any values of the parameters, and so the second condition of the theorem we are applying is satisfied, and we can conclude that a unique and stable limit cycle bifurcates from the equilibrium point
through a Hopf Bifurcation for
. □
Now, we include some numerical experiments. We fix the parameters as follows:
In this case,
and the eigenvalues of
are
In
Figure 4, we represent the case
, in which the equilibrium is locally asymptotically stable. In
Figure 5, we represent the case
, in which the equilibrium loses stability and a limit cycle arises due to Hopf bifurcation.
We conclude this section by considering a case in which there are no periodic orbits in :
Theorem 3. then system (4) does not admit periodic orbits in the set . Proof. In order to prove the nonexistence of periodic orbits, we use the Bendixson–Dulac theorem that states that if there exists a function
such that the term
does not change sign in a simply connected set
, then there are no periodic orbits on
.
We consider then function
; then:
We observe that, since
for
, there are no periodic orbits in the set
and for the same reason, there are no periodic orbits crossing the half line
. As a consequence we will restrict to the case
for which we obtain
Then, in if , and we conclude that there are no periodic orbits in the whole set . □
Remark 1. We observe that, for , the condition of Theorem 3 on the nonexistence of periodic orbits and then is compatible with the results of bifurcation analysis.
3. Existence and Stability Analysis of Equilibria
The first step for studying the dynamics of the system (
2) is to find all the equilibrium points and analyzing their stability.
Theorem 4. System (2) has the following boundary equilibria on : for any values of the parameters,
for any values of the parameters,
If , the equilibrium with
Proof. From direct calculation. □
We also analyze the existence of nontrivial positive equilibria for system (
2).
Theorem 5. The system (2) has at least one positive equilibrium if and only if where is a solution of the equation with the coefficients defined below.
Proof. By direct calculation, we obtain the equilibrium
, with
and
is the solution of the equation
We must require that
verifies
so that
is positive. This is always verified if
Otherwise, this is satisfied when
The expressions of the coefficients in Equation (
21) are given by
We apply Descartes rule of signs in Equation (
21) to study the existence of positive roots. The coefficient of degree zero,
, is always negative and the coefficient of degree four,
, is always positive. This means that there always exist a real positive and a real negative zero of the polynomial. The other two zeroes can be complex or real with the same sign.
The other coefficients of the polynomial can be positive, negative or zero, and combining all the possible signs we obtain that:
So that there exist at least one non-trivial interior equilibrium
of the system (
2) if condition (
19) is satisfied. □
Corollary 1. A sufficient condition for the system (2) has at least one positive equilibrium is . Remark 2. The positive equilibria, if they exist, are all on the plane .
Remark 3. Numerically, we have found different cases acording to the existence of positive solutions for the Equation (21) and the existence of positive equilibria . There exist systems for which there is only one positive solution of Equation (21), and this solution is such that verifies condition (19). As an example the system with parameters: There are systems for which there is only one positive solution of Equation (21), and this solution does not verify condition (19), so there are no positive equilibria for system (2). As an example the system with parameters: There are systems for which there are three positive solutions of Equation (21), and only for one of this solutions condition (19) holds, so there are one positive equilibrium for system (2). As an example the system with parameters: Also there are systems for which there are three positive solutions of Equation (21), and none of them verifies condition (19), so system (2) has none positive equilibria. As an example the system with parameters:
We have not found, numerically, any conditions for which Equation (21) has three positive solutions and the three of them verify condition (19) but we are not able to exclude this case analytically. We will analyze the local stability of the equilibria. To do this we consider the Jacobian functional of the vector field
where
and where we have set for simplicity
.
Theorem 6. The stability of the boundary equilibria is the following:
The equilibrium point is always a saddle, the z-axis is the stable manifold and it has an unstable manifold of dimension two.
The equilibrium point is a saddle with a stable manifold of dimension one if and with a stable manifold of dimension two if .
The equilibrium point is unstable if or if and one of the following statements holds:
and condition (7) holds, and .
The equilibrium point is asymptotically stable if and one of the following statements holds:
and condition (8) holds, and .
where the coefficients are the ones given in (6).
Proof. The local stability analysis of the equilibria and is the same as in the case without indirect effects.
For the equilibrium
the Jacobian matrix is
so there are two positive eigenvalues and one negative eigenvalue so that
is a saddle with a stable manifold of dimension one, which is the
z-axis, and an unstable manifold of dimension two.
For the equilibrium
we have
is a saddle, with two positive eigenvalues and one negative eigenvalue if
and with one positive eigenvalue and two negative eigenvalues if
When the value of
surpasses the value of
, a new equilibrium
appears, and the equilibrium
changes from having an unstable manifold of dimension two to an unstable manifold of dimension one. The functional Jacobian of
is
with eigenvalues
with the coefficients
given in (
6).
The sign of eigenvalues
has been analyzed in
Section 2. We note that the expression
can be written as
, so we conclude that the eigenvalue
is positive if
and it is negative in the other case. Combining the different possibilities for the three eigenvalues we obtain the conditions for the stability of
. □
Remark 4. We recall that in the case , the equilibrium point has two non-zero eigenvalues in the plane , as stated in Section 2, and the third eigenvalue is zero. The direction on the flow in the y-direction is determined by the z coordinate. Note that in this case , so is positive if and is negative if . Remark 5. We have seen that if (22) is verified then there exists at least one positive equilibrium . We can prove that, at least in this case, the first eigenvalue of is positive, that is is unstable. We recall that is positive if In this case we observe that , in fact which is always verified.
When all boundary equilibrium points are unstable we expect that there exist solutions that are attracted by an interior equilibrium or a limit cycle.
Theorem 7. The positive equilibrium is locally asymptotically stable if and only if and , where the constants are defined below.
Proof. The characteristic polynomial of
is given by the following expression:
where
it is,
We use Hurwitz criterion to study the stability of the equilibrium point , and we conclude that it is asymptotically stable if and only if and also the expression is positive. □
Remark 6. We observe that if then is negative, that is . This means that at least one eigenvalue has negative real part.
4. Some Remarks about Coexistence of the Three Species
The problem of the coexistence of the three species can be reformulated in mathematical terms by finding the conditions for which the positive solutions starting in the interior of
do not approach the boundary of the set,
as
. These ideas can be made rigorous in the context of persistence theory (see [
16]). There are several definitions that are used by mathematicians depending on the context. If we consider a nonlinear system of the form
then we have persistence (see [
30]) if
while we have permanence (see [
31]) if there exists
, independent of
, such that
Finally, we have uniform persistence (see [
32]) if there exists
, independent of
, such that
In many cases, the favorite choice for analysis is uniform persistence, since in real cases, requiring that
is not sufficient. In fact, a small stochastic or non-autonomous perturbations may lead solutions converges to the boundary. For this reason, in general, it is important to require that
. We recall the definition introduced in [
17]:
Definition 1. The system (29) is uniformly persistent if there exists such that In order to prove uniform persistence (see Theorem 8.17 page 188, hypothesis (H) page 185 and Theorem 5.2 page 126 in [
16]), we would need to prove the following conditions:
Hypothesis 1 (H1). There exists a compact attractor of bounded set.
Hypothesis 2 (H2). The invariant sets of are weakly repelling.
Hypothesis 3 (H3). The invariant sets of are acyclic.
Unfortunately, in the case of system (
29), we cannot prove the existence of a compact attractor of bounded set of
and as a consequence we are not able to obtain uniform persistence of the system. This lack of dissipativity can be guessed by considering the unbounded solutions on the invariant plane
. Moreover, we recall that the divergence of the vector field
F associated to the system
is related to the evolution of the three dimensional volumes elements under the flow of the system. We observe that it has a complex expression, in particular, there are no values of the parameters for which it has a negative sign (which means contraction), on a region of
. Moreover, while the first term and the second term are negative for high values of
x and
z respectively, the third term is positive for high value of
y.
Then, we are able only to prove conditions (H2)–(H3) which only ensures that the invariant sets of does not attracts positive solutions. For condition (H2) we first recall the definition:
Definition 2. A set is called is called weakly repelling if there is no solution of system (29) starting at , with , such that as . Theorem 8. We suppose that hypothesis of Theorem 3 is verified. Then the invariant set of are weakly repelling in any of the following cases
- (1)
;
- (2)
and .
Proof. In order to prove the theorem, it is sufficient to show that the stable manifolds of the invariant sets of are contained all contained in .
The equilibria
and
have their stable manifolds on
for any value of the parameters. Then if
there are no further equilibria. Otherwise, if
exists, a sufficient condition that ensures that its stable manifold is in
consists in requiring that
that is, its first eigenvalues is positive. To conclude the proof we have to exclude the existence of other invariant set contained in
. We use Theorem 3 that guarantees the non existence of periodic orbits in the plane
, that is the only invariant sets of
are the equilibria. □
Now we pass to check condition (H3).
Definition 3. Let . Then A is chained to B in , and we write if there exists a total trajectory with such that Definition 4. A finite collection of subsets of is called cyclic if after possibly renumbering or in for some . Otherwise it is called acyclic.
Thanks to this notion we are able to exclude the case in which orbits are attracted by an heteroclinic cycle that connects the equilibria such as in the well known case of [
33].
Theorem 9. If hypothesis of Theorem 3 is verified, then the invariant sets of are acyclic.
Proof. By hypothesis the only invariant sets of
are equilibria. The analysis of the previous sections is sufficient to exclude the existence of cycle between equilibria. In
Figure 6 below we represent the case in which there are only two boundary equilibria (
and
) while in
Figure 7 we represent the case in which we have also
. In the latter situation we distinguish two cases: both eigenvalues with positive and both with negative real part respectively. □
In conclusion, we are not able to prove uniform persistence, however the above results guarantee a sort of weak persistence of the three species.
5. Hopf Bifurcation
In this section, we analyze the possible existence of a limit cycle by Hopf bifurcation for the positive equilibrium . We recall that Hopf bifurcation occurs when a pair of complex conjugate eigenvalues of the Jacobian matrix of an equilibrium crosses the imaginary axis. In this case a limit cycle arises and its stability character can be obtained by the analysis of the first Lyapunov coefficient. If it is negative the cycle is stable and the bifurcation is called supercritical; otherwise, it is unstable and the bifurcation is called subcritical.
Because of the complexity of the system and the high number of parameters we are not able to perform a general bifurcation analysis. In this section, we simplify this task by fixing the value of parameters and using m as bifurcation parameter.
With the above choice of the parameters we have
and as a consequence
. In this case the equilibrium is:
where
and
is the solution of the equation
with
Moreover the functional Jacobian at
is
where
The characteristic polynomial is
where
The Hurwitz matrix of the characteristic polynomial is given by
If , we always have a negative (real) eigenvalue, while if , that is , we ensure that the other two eigenvalues are complex conjugate. If (resp. ) then we have at least two eigenvalues with negative (resp. positive) real part.
From Hurwitz-Routh criterion we obtain that a necessary condition for Hopf Bifurcation in this case becomes
that is
We have numerically solved the previous equation by using the software Matlab and we obtained the critical value
For this value we obtain , , and .
For
the eigenvalues of the Jacobian matrix
are
We will see below that as m passes trough the value the real part of the eigenvalues change sign from negative to positive.
Then a cycle appears due to Hopf Bifurcation of the equilibrium
. We are not able to exactly compute the first Lyapunov exponent of the system, however simulations and the sign of the term
suggests that the cycle is stable and the equilibrium looses stability. In
Figure 8 below we represents the solutions for
, as we expect, they converge to a periodic solution of period
.
In order to illustrate the results of this section we present several numerical simulations. We fix the parameters as above and initial data as follows.
In a first numerical experiment we fix
and as we expect solutions converges to the equilibrium (see
Figure 9 below)
In this case the eigenvalues of
are
and have negative real parts.
In a second numerical experiment, we set
. The positive equilibrium point
is unstable, the eigenvalues of
are
We observe that the eigenvalues
now have positive real part and a Hopf bifurcation occurs at
. Solutions converge to a limit cycle as shown in
Figure 10.
6. Conclusions
In this paper, we have considered a model describing the dynamics of an ecological system with two prey species and a predator species, which is a modification of the model proposed in [
1]. In particular, due to the high availability of one of the two prey populations, we supposed that the predator prefer to predate the more available prey population whereas the other one take advantages of it. In previous articles in the literature (see Introduction for details), it has been pointed out the importance of indirect effects in order to describe real cases of coexistence.
We have performed the stability analysis of equilibria and we have made a detailed analysis of the system on the invariant planes, including the study of the existence of Hopf bifurcation at the equilibrium point . We have proved that through this bifurcation a stable limit cycle appears.
Regarding the existence of positive equilibria, the expression (
21) and the verification of one of the conditions (
22) or (
23) give us all the positive equilibrium points. The expressions of the equilibria as a function of the parameters are too complicated and not easy to handle. Because of this, we have not been able to determine, in general, for which conditions appear none, one or three equilibria. We have obtained sufficient conditions for the existence of at least one positive equilibrium (see Corollary 1). Also, for fixed values of the parameters it is easy to compute the positive equilibria and maybe it would be possible for certain subfamilies on which less parameters are considered.We have found values of the parameters for which there exist one equilibrium point, and others for which there are not any positive equilibria, but numerically, we have not found values for which three positive equilibria exist, and therefore we think that probably this situation is not feasible, although we have not been able to prove it.
We have shown that Hopf bifurcation can occur also at the positive equilibrium
and as a consequence, coexistence of the three species via the existence of an attracting limit cycle is possible, by taking into account indirect effects of predation. It is worth mentioning that in [
1] Hopf bifurcation is obtained only by considering a version of the model with time delay.
Furthermore, we have included a detailed discussion about the problem of persistence of the system.
Throughout the paper, several numerical simulations are provided in order to illustrate the theoretical results.
Due to the complexity of the model, we have not been able to perform a complete bifurcation analysis, then this point remains as an open problem and it would be interesting to study it in the future. A further analysis of predator prey models incorporating indirect effects can be done, for example, considering time delay or non autonomous (seasonal) terms.
A final comment for further research: the proposed model in this paper could be easily implemented by using analog devices ([
34]). Therefore, it could be used to perform experiments in this way.