A Nutrient-Prey-Predator Model: Stability and Bifurcations

# A Nutrient-Prey-Predator Model: Stability and Bifurcations

Mary Ballyk Keywords: Chemostat, Hopf bifurcation, coexistence equilibrium    Ibrahim Jawarneh

Department of Mathematical Sciences
New Mexico State University
Las Cruces, NM 88003, USA
Mathematics Subject Classification: Primary 37G10; Secondary 34C23 92D25 34A34
Ross Staffeldt
###### Abstract

In this paper we consider a model of a nutrient-prey-predator system in a chemostat with general functional responses, using the input concentration of nutrient as the bifurcation parameter. We study the changes in the existence of isolated equilibria and in their stability, as well as the global dynamics, as the nutrient concentration varies. The bifurcations of the system are analytically verified and we identify conditions under which an equilibrium undergoes a Hopf bifurcation and a limit cycle appears. Numerical simulations for specific functional responses illustrate the general results.

## 1 Introduction

We consider a mathematical model of two-species predator-prey interaction in the chemostat under nutrient limitation. With the exception of one nutrient, all nutrients that the prey species requires are supplied to the growth vessel from the feed vessel in ample supply. The predator species grows exclusively on the prey. With the concentration of the limiting nutrient, the concentration of prey (say, phytoplankton), and the concentration of predator (say, zooplankton), we consider the following model:

 dN/dt =(μ−N)D−Pf1(N) dP/dt =γ1Pf1(N)−D1P−Zf2(P) (1) dZ/dt =γ2Zf2(P)−D2Z

for initial conditions , , and .

The concentration of the growth-limiting nutrient in the feed vessel is denoted , and will be the bifurcation parameter in our analysis. is the input rate from the feed vessel to the growth vessel as well as the washout rate from the growth vessel to the receptacle, so that constant volume is maintained. The parameters and are the removal rates of and , respectively, from the growth vessel, incorporating the washout rate and the intrinsic death rates of and . Our analysis does not necessarily require that and are positive; however, and should be positive. The yield coefficient gives the amount of prey biomass produced per unit of nutrient consumed, while gives the amount of predator biomass produced per unit of prey biomass consumed.

The function represents the per capita consumption rate of nutrient by the prey populations as a function of the concentration of available nutrient; similarly, the function represents the per capita consumption rate of the prey by the predator as a function of available prey. These functions are assumed to satisfy , , 2, for all , and for all . We further assume that and are bounded. To avoid the case of washout due to an inadequate resource, we assume that , and to avoid the case of an inadequate prey, we assume that . Define and to be the unique numbers satisfying

 (2)

The number is the break-even concentration of nutrient, at which the growth and removal of phytoplankton balance; is similarly interpreted. The number plays a central role in our investigation, so we assume that and are also defined. From the perspective of , due to the boundedness assumptions on and , and are perturbations of .

###### Lemma 1.1.

From a functional point of view, and are right inverses of and , respectively. Accordingly, on their respective domains and are as differentiable as and . We note for later use

 λ′P(D1)=(γ1⋅f′1(λP(D1)))−1andλ′Z(D2)=(γ2⋅f′2(λP(D2)))−1.\qed (3)

Kuang and Li  studied this system with general functional responses and distinct values of , , and . However, they fixed the input nutrient concentration, whereas we have this as a parameter. With the hypothesis that , they provide numerical criteria for the stability of a coexistence equilibrium, and prove that a cycle exists when the equilibrium is unstable [7, Theorem 3.2]. When the hypothesis does not hold, they provide numerical evidence that stability of the coexistence equilibrium breaks down and a cycle appears. A similar model was studied in  with functional response of Holling type I and of Holling type II, demonstrating the existence of a Hopf bifurcation in response to varying nutrient concentration. These results inspire our work. Our goal is to prove analytically that the system undergoes a Hopf bifurcation without restricting the forms of the uptake functions or the values of the removal rates and .

This paper is organized as follows. In section 2, conditions for existence and local stability of predator-free equilibria are obtained for general functional responses and . In section 3, we study stability of a coexistence equilibrium which appears as the parameter increases. We quote a version of the Hopf bifurcation theorem, stating the result in a limited form most useful for us. Accordingly, application of this theorem requires us to develop specific information about the behavior of eigenvalues of the linearizations as the paramter increases. In section 4, we use a smooth change of variables to reach a situation where the Hopf bifurcation theorem in three dimensions applies, enabling us to conclude the existence of cycles. In section 5, our results are illustrated using simulations arising by choosing rate functions of Holling type II and Holling type III forms. In section 6, we explain in detail the approximations and estimates that support our work in the latter part of section 3.

## 2 Steady States and Their Stability

To begin, we establish that the solutions of (1) are nonnegative and bounded. These are minimum requirements for a reasonable model of the chemostat. We then develop conditions for the existence and stability of equilibria. We conclude the section by proving uniform persistence in the sense of  when is sufficiently large and the initial values and are positive.

###### Lemma 2.1.

All solutions , , and of (1) are nonnegative and bounded.

###### Proof.

The plane is invariant for system (1). Therefore, by existence and uniqueness, if then for all . Similarly, since , the plane is invariant, so implies for all .

Suppose . If there exists a with , then there is a least such number, say . Then since . Consequently, there is such that , which is a contradiction to the choice of .

For the boundedness of solutions, set . From system (1), it follows that

 U′(t)=Dμ−DN(t)−γ−11D1P(t)−(γ1γ2)−1D2Z(t)≤Dμ−ˆDU(t),

where . Then

 U(t) ≤(Dμ)/ˆD+(U(0)−(Dμ)/ˆD)e−ˆDt ≤{U(0),if U(0)>(Dμ)/ˆD,(Dμ)/ˆD,if U(0)≤(Dμ)/ˆD.

Since , , and are nonnegative, the boundedness of implies the boundedness of , , and . ∎

There are at most three biologically relevant equilibria of system (1) depending on the value of . The equilibria and the conditions of their existence are summarized in the following theorem.

###### Theorem 2.2.

Let and be as in (2). The equilibria of the system (1) satisfy the following conditions:

1. The washout equilibrium exists for all .

2. The single-species equilibrium exists for all , where

 P(μ,D1)=(μ−λP(D1))Dγ1D1. (4)
3. The coexistence equilibrium

 E2(μ,D1,D2)=(N(μ,D1,D2),λZ(D2),Z(μ,D1,D2))

exists for all , where

 μc1(D1,D2)=λP(D1)+D1λZ(D2)Dγ1 (5)

and , satisfy the simultaneous equations

 (μ−N(μ,D1,D2))D−λZ(D2)f1(N(μ,D1,D2))=0, (6) (7)

Thus, for only the equilibrium exists; for , there are two equilibria ; and, when , there are three equilibria , , and .

###### Proof.

From the -equation, either or (so that ). If , the -equation yields

 0=γ1Pf1(N)−D1P,

which implies either or (so that ).

If and , then the -equation gives , so that . Thus, the washout equilibrium is given by and exists for . Note that as increases, moves along the -axis in . This is the proof of the first part of the theorem.

If and , then in the -equation gives

 0=(μ−λP(D1))D−PD1γ1,%withsolutionP=(μ−λP(D1))Dγ1D1:=P(μ,D1). (8)

Note that for all , and so the single-species equilibrium exists in the positive cone for all . This is the proof of the second part of the theorem.

If , then in the - and -equations, and equations (6) and (7) define and as implicit functions of , , and .

We now determine the critical value of at which first appears in the positive cone. Equation (7) implies

 Z(μ,D1,D2)=γ2D2λZ(D2)(γ1f1(N(μ,D1,D2))−D1). (9)

Since is strictly increasing, if and only if , and for . Substituting in (6), we obtain the equation

 0=(μ−λP(D1))D−λZ(D2)D1γ1with solutionμ=λP(D1)+D1λZ(D2)Dγ1:=μc1(D1,D2).

Thus, is positive when , and the coexistence equilibrium exists in the positive cone for all . This proves part three of the theorem. ∎

In Theorems 2.3, 2.4, and 2.5 we investigate the stability of the equilibria of system (1) by finding the eigenvalues of the associated Jacobian matrices. The Jacobian matrix of the system (1) takes the form

 J=⎡⎢⎣−D−f1(N)0γ1Pf′1(N)γ1f1(N)−D1−Zf′2(P)−f2(P)0γ2Zf′2(P)γ2f2(P)−D2⎤⎥⎦. (10)

We first summarize the stability of in the following theorem. Here, the breakeven concentration of nutrient given in (2) plays a critical role.

###### Theorem 2.3.

The equilibrium point is locally asymptotically stable if and unstable if . When , is globally asymptotically stable with respect to solutions initiating in . That is, the plane is , the stable manifold of .

###### Proof.

The Jacobian at is

 J(E0)=⎡⎢⎣−D−f1(μ)00γ1f1(μ)−D1000−D2⎤⎥⎦,

so that the eigenvalues of are , and . The stability of now follows from (2) and the fact that is strictly increasing: when and when . To see that when , consider the Lyapunov function

 L(N,Z)=(μ−N)22+Z22.

Clearly, and if . The time derivative of at a point on a trajectory of system (1) is

 L′(N,Z)=−D(μ−N)2−D2Z2<0

for . Thus is globally asymptotically stable in the plane . ∎

For , , so that and coalesce (see equation (4)). When , enters the positive cone. We summarize the stability of in the following theorem. Note that the critical value of given in (5) now plays a central role.

###### Theorem 2.4.

The equilibrium point is locally stable if and unstable if . When , is globally asymptotically stable with respect to solutions initiating in . That is, the plane is , the stable manifold of .

###### Proof.

The Jacobian matrix at is

 J(E1(μ,D1))=⎡⎢ ⎢⎣−D−P(μ,D1)f′1(λP(D1))−f1(λP(D1))0γ1P(μ,D1)f′1(λP(D1))0−f2(P(μ,D1)))00γ2f2(P(μ,D1))−D2⎤⎥ ⎥⎦.

The determinant of the upper lefthand -by- submatrix is positive and its trace is negative, so its eigenvalues have negative real parts. The third eigenvalue is . If , so that , then , and is locally stable. Similarly, if , so that , then and is unstable with one dimension of instability.

To see that when , consider the Lyapunov function introduced by Hsu in 

 L(N,P)=∫NλP(D1)f1(n)−f1(λP(D1))f1(n)dn+1γ1(P−P(μ,D1)−P(μ,D1)ln(P/P(μ,D1))).

Notice that ,

 ∂L∂N=f1(N)−f1(λP(D1))f1(N)=0and∂L∂P=1γ1(1−P(μ,D1)/P)=0

precisely when . Moreover,

 ∂2L∂N2(λP(D1),P(μ,D1)) =f′1(λP(D1))f1(λP(D1))>0 and ∂2L∂P2(λP(D1),P(μ,D1)) =1γ1P(μ,D1)>0.

Therefore, is the only critical point of and it is a local minimum, so that for all .

Now we compute the time derivative of at a point along a trajectory of system (1). Noting from (2) and (4) that

we have

If , then and , so

 μ−Nμ−λP(D1)−f1(N)f1(λP(D1))>0.

Thus, when .

If , then . When ,

 μ−Nμ−λP(D1)<1

while for , then . In either case, , and so when .

Finally, if and only if . By LaSalle’s extension theorem , any trajectory of system (1) in the plane for which approaches the largest invariant set in the line , and this is simply . Therefore, is globally asymptotically stable in the plane . ∎

When , and coalesce. With , we have (see equations (5) and (4)). Also , so that (see the discussion around (9)). Thus,

 E2(μc1(D1,D2),D1,D2)=(λP(D1),λZ(D2),0)=E1(μc1(D1,D2),D1).

Said another way, as increases through , enters the positive cone by passing through .

With , the Jacobian matrix at takes the form

 (11)

The eigenvalues of satisfy

 x3+a1x2+a2x+a3=0,

where

 a1(μ,D1,D2) −γ1f1(N(μ,D1,D2))+D1+D, (12) a2(μ,D1,D2) =λZ(D2)Z(μ,D1,D2)f′1(N(μ,D1,D2))f′2(λZ(D2)) +D2Z(μ,D1,D2)f′2(λZ(D2))+DZ(μ,D1,D2)f′2(λZ(D2)) +D1λZ(D2)f′1(N(μ,D1,D2))−Dγ1f1(N(μ,D1,D2))+DD1, (13) and a3(μ,D1,D2) =D2Z(μ,D1,D2)f′2(λZ(D2))(D+λZ(D2)f′1(N(μ,D1,D2))). (14)
###### Theorem 2.5.

The coexistence equilibrium point

 E2(μ,D1,D2)=E2(N(μ,D1,D2),λZ(D2),Z(μ,D1,D2))

is asymptotically stable if and only if and .

###### Proof.

Since is positive, this follows from the Routh-Hurwitz criterion. ∎

We conclude this section with a significant strengthening of Lemma 2.1. We use the concept of uniform persistence introduced in .

###### Theorem 2.6.

Assume that . Then system (1) is uniformly persistent with respect to all solutions satisfying and .

###### Proof.

Recall from Lemma 2.1 that all solutions of system (1) for which and are positive and bounded.

We first show that . If and , then . But this is impossible, for then it follows from the -equation that as and this, in turn, contradicts the fact that is bounded.

Now, suppose while . Then there exists a sequence of local minima of satisfying as . Thus,

1. , since is a local minimum, and

2. as , since .

From the -equation we have

 N′(τn)=μD−(N(τn)D+P(τn)f1(N(τn))),

so that

Rearranging and using the facts that as , is continuous and , we get

 0=limτn→∞∣∣N(τn)D+P(τn)f1(N(τn))∣∣≥μD>0,

Choose . Then is a nonempty, compact invariant set with respect to system (1). We claim and are not in . Suppose . Since , is an unstable hyperbolic equilibrium point. By theorem 2.3, is globally asymptotically stable with respect to solutions initiating in the plane . Since , . By the Butler-McGehee lemma (see lemma A1, ), there exists , so that . For such an initial condition, the governing system is , and . But then becomes unbounded as . This is a contradiction to the compactness of , and so .

Suppose is in . Since , is an unstable hyperbolic equilibrium point. By theorem 2.4 is globally asymptotically stable with respect to solutions initiating in the plane . Since , . By the Butler-McGehee lemma, there exists , so that . If , then , a contradiction. Thus , and this implies is unbounded as , contradicting the compactness of . Therefore, .

Suppose the system (1) is not persistent. Then there exists such that or . Then , which implies either or , neither of which can be true. Thus, and , and it follows from the main result of  that system (1) is uniformly persistent. ∎

## 3 Stability of the coexistence equilibrium

In this section we study the coexistence equilibrium point

 E2(μ,D1,D2)=(N(μ,D1,D2),λZ(D2),Z(μ,D1,D2))

first when and varies and then after relaxing the assumption . In particular, we study the eigenvalues of the Jacobians at the coexistence equilibria as increases. To prepare for the study of the evolution of the equilibrium , we observe the following consequence of the implicit function theorem [5, p.122].

###### Lemma 3.1.

For , , and there is a local parameterization of the locus of coexistence equilibria

 E2(μ,D1,D2)=(N(μ,D1,D2),λZ(D2),Z(μ,D1,D2))

defined on an interval containing and a disc around that is smooth to the smaller of the degree of smoothness of and the degree of smoothness of .

###### Proof.

Set in the - and - equations of system (1). With from (2), let

 G1(N,Z,μ,D1,D2)=D(μ−N)−f1(N)λZ(D2)

and

 G2(N,Z,μ,D1,D2)=γ1f1(N)λZ(D2)−D1λZ(D2)−(D2/γ2)Z,

and define by