Open access peer-reviewed chapter

Stability and Hopf Bifurcation Analysis of a Simple Nutrient- Prey-Predator Model with Intratrophic Predation in Chemostat

By Zabidin Salleh and Liyana Abd Rahim

Submitted: April 8th 2017Reviewed: October 11th 2017Published: December 20th 2017

DOI: 10.5772/intechopen.71624

Downloaded: 272

Abstract

A 3-dimensional nutrient-prey-predator model with intratrophic predation is proposed and studied. Some elementary properties such as invariance of nonnegativity, boundedness and dissipativity of the system are presented. The purpose of this chapter is to study the existence and stability of equilibria along with the effects of intratrophic predation towards the positions and stability of those equilibria of the system. We also investigate the occurrence of Hopf bifurcation. In the case when there is no presence of predator organisms, intratrophic predation may not give impact on the stability of equilibria of the system. We also analysed global stability of the equilibrium point. A suitable Lyapunov function is defined for global stability analysis and some results of persistence analysis are presented for the existence of positive interior equilibrium point. Besides that, Hopf bifurcation analysis of the system are demonstrated.

Keywords

  • chemostat
  • intratrophic predation
  • local and global stability
  • Hopf bifurcation

1. Introduction

Chemostat, a piece of laboratory apparatus is frequently used in mathematical ecology. This device carries an important role for ecological studies because the mathematics are tractable, and the relevant parameters are readily measurable. Chemostat also can be used for modelling microorganism systems such as lake and wastewater treatment. For detailed information regarding chemostat device, see [1].

There are abundant of research journals and papers for analysing chemostat models. One of them is Li and Kuang [2], considered a simple food chain with distinct removal rate, which in this case, conservation law has failed. Therefore, to overcome the problem, they constructed a Lyapunov function in global stability analysis of predator-free steady state. Local and global stability of other steady states were shown in the paper along with persistence analysis of the system.

Another paper regarding chemostat model is by El-Sheikh and Mahrouf [3] that presented a 4-dimensional food chain in a chemostat with removal rates. They studied local and global stability of equilibria along with elementary properties including boundedness of solutions, invariance of nonnegativity, dissipativity and persistence analysis. Hopf bifurcation theory was applied.

Recently, several analysis on chemostat models are carried out, for example, by Hamra and Yadi [4] and Yang et al. [5]. In the work of [4], they studied a chemostat model with constant recycle sludge concentration. Number of parameters are reduced by considering a dimensionless model. Next, they successfully proved the existence of a positive uniform attractor for the model with different removal rates by using dissipative theory. Hence, they used methods of singular perturbation theory in order to investigate the asymptotic behaviour of the chemostat model under small pertubations. Thus, it is shown that in the case of two species in competition, the positive unique equilibrium point is globally asymptotically stable.

In the research of Yang et al. [5], piecewise chemostat models which involve control strategy with threshold window are proposed and analysed. They investigated the qualitative analysis such as existence and stability of equilibrium points of the system and it is proved that the regular equilibria and pseudo-equilibrium cannot coexist. The global stability analysis of both the regular equilibria and pseudo-equilibrium have been studied using qualitative analysis techniques of non-smooth Filippov dynamic systems. Furthermore, the bifurcation analysis of the system is investigated with theoretical and numerical techniques.

Moreover, one of the interesting topic is the research involving intratrophic predation. Intratrophic predation is a situation where members of a trophic group consume other members of the same trophic group (for the purpose of mathematical modelling). Several previous studies based on intratrophic predation have been discovered. One of them is by Pitchford and Brindley [6]. They studied a predator-prey model with intratrophic predation and successfully shown that the model had desirable and plausible features. Furthermore, to investigate the effect of intratrophic predation towards the position and stability of equilibrium points of a model, they had developed a simple asymptotic method.

For Hopf bifurcation analysis, a study by Mada Sanjaya et al. [7] has been carried out. They introduced an ecological model of three species food chain with Holling Type-III functional response. Two equilibrium points are obtained. They proved that the system has periodic solutions around those equilibrium points. Not only that, they also investigated the dynamical behaviours of the system and found that it was sensitive when the parameter values varied.

Tee and Salleh [8] investigated Hopf bifurcation of a nonlinear modified Lorenz system using normal form theory that was the same technique used in Hassard et al. [9]. Then, the dynamics on centre manifold of the system was presented as it will be applied in the technique of normal form theory. Another research based on Hopf bifurcation analysis is carried out by [10]. They considered a three-species food chain models with group defence. It is proved that the model without delay undergone Hopf bifurcation by using the carrying capacity of the environment as the bifurcation parameter. In the analysis of Hopf bifurcation in a delay model, they used a computer code BIFDD to determine the stability of bifurcation solutions.

Furthermore, a simple 3-dimensional food chain model in chemostat with variable yield for prey population and constant yield for predator population is proposed by Rana et al. [11]. In this model, the prey consumes the nutrient and the predator consumes the prey but the predator does not consume the nutrient. The functional response functions are assumed in Michaelis-Menton type. The stability of equilibrium points, the existence of limit cycles, the Hopf bifurcation and the positive invariant set for the system are discussed by qualitative analysis of differential equations. Finally, numerical simulations are carried out in support of the theoretical results.

Our work is a modification of the models in Li and Kuang [2] and El-Sheikh and Mahrouf [3]. It should be emphasised that this work is different from [2, 3]. The modified model contains parameter b, known as the measure of intensity of intratrophic predation which is not in their models. The parameter band term1+D1x1+x+D1byare added to differential equations xand yin the interest of studying the behaviour of the modified model. By this motivation, we analysed the stability and Hopf bifurcation of the model with intratrophic predation, as intratrophic predation analysis is rarely considered in mathematical models of populations biologically [6].

The purpose of this chapter is to study the existence and stability of equilibria along with the effects of intratrophic predation towards the positions and stability of those equilibria of the system. We shall also investigate the occurrence of Hopf bifurcation towards the system. Hence, we introduced a simple nutrient-prey-predator model in chemostat with intratrophic predation. Some notations regarding the model will be explained. Not only that, several elementary properties such as nonnegativity, boundedness and dissipativity along with definitions and several results will be presented. The local and global stability together with existence of the equilibrium points are shown. For global stability analysis, a suitable Lyapunov function is defined. Lastly, we applied Hopf bifurcation theorems (see [9]) in the analysis of Hopf bifurcation.

2. The model

In this work, we consider a nutrient-prey-predator model with one prey organism and one predator organism in chemostat with intratrophic predation. Biologically, the predator organisms will feed upon the prey organisms, while the prey organisms will consume the nutrient in the chemostat. Precisely, the modified model of [2, 3] is

{s(t)=(s0s)D01β1f1(s)x,x(t)=( f1(s)D1)x1β2f2(x)(y+D1xyD2+x+by),y(t)=f2(x)(y+D1xyD2+x+by)y,E1

where s0=s0>0,x0=x0>0,y0=y0>0,Di=1;i=0,1,2,βj1;j=1,2and 0b1.This leads to the following assumptions on functional responsefn,n=1,2:

ifn:R+R+andfn,n=1,2arecontinuously differential equations,iifn0=0,iiif1s>0,f2x>0foralls,x>0.E2

Generally, functional response is a common component in predator-prey system. The term ‘functional response’ was first stated by Solomon [12] which defined the relationship between the rate of predation, i.e., the number of prey organisms consumed per predator organism in time, twith the density of prey organisms.

By referring to system (1), srepresents the concentration of nutrient at time twhile s0is the input nutrient concentration. The variables xand yrepresent the concentration of prey and predator at time t, respectively. Parameters β1and β2denote the yield constants, D0is the washout rate of the chemostat while D1and D2are the removal rates of xand y, respectively. Removal rate is the sum of washout rate and death rate θi, i.e., Di=D0+θi,i=1,2.f1sand f2xdenote the specific growth rate of prey xand predator y, respectively, while bis the measure of intensity of intratrophic predation in predator organisms y.These observations are based on numerical simulations. We can rescale system (1) by reducing the number of parameters using standard change of variables such as.

s¯=ss0,x¯=xD0β1s0,y¯=yD0D1β1β2s0,
f1s=f1~s0s¯D0,f2x=f2~D0β1s0x¯D0,b=b¯β2,
t¯=Dit,i=0,1,D2=D0β1s0.

To make our works more convenient, we just drop the bars and tildes. So after rescaling, the system (1) becomes

s=1sf1sx,x=f1sD11xf2xy1+D1x1+x+D1by,y=f2xy1+D1x1+x+D1byy,E3

where

s0>0,x0>0,y0>0andD1>0.E4

Now, the variables are non-dimensional and the discussion is in R+3=sxy:s>0x>0y>0.

2.1. Elementary properties of system (3)

In this section, we present nonnegativity, boundedness and dissipativity of the system (3) with respect to a region inR+3. Firstly, we consider dissipativity of the system (3).

[13]. A system with differential equations x=fxis defined to be dissipative if there exists a bounded subset Γof R3, such that there is a time t0for any x0R3which depends on x0and Γso that the solution of the system ϕtx0Γfor tt0.

Let Hbe the region

H=sxyR+3:1Pmaxqs+x+y1Pmin+q,

where Pmax=max1f1s+f1sD11,Pmin=min1f1s+f1sD11and qis a positive constant.

Then,

  1. Ηis positively invariant.

  2. All nonnegative solutions of (3) with initial values in R+3are uniformly bounded and they eventually attracted into region Η.

  3. The system is dissipative.

We must show that the solutions of system (3) are nonnegative and bounded, so that the system becomes biologically meaningful.

Proof of Theorem 1.

(i): First, we must show that 1Pmaxqs+x+y1Pmin+q.Let us define

Pmax=max1f1s+f1sD11 and Pmin=min1f1s+f1sD11.

By adding those differential equations s,xand yin (3), we shall get

s+x+y=1s+f1s+f1sD11x+y.

Hence, this leads to

1Pmaxs+x+y1s+f1s+f1sD11x+y1Pmins+x+y.

Thus, solving the above inequality gives

1Pmaxqs+x+y1Pmin+q,

as qis the positive constant

(ii) and (iii): Let 0<s0<1and consider.

s=1sf1sx<1s.

Then st<1et1s0for 0<s0<1,and hence limtsupst<1for 0<s0<1.

Now, consider the xequation;

x=f1sD11xf2xy1+D1x1+x+D1by<f1s1x<α1x,

where α1=maxsΗ f1s1.We assume that α1=maxsΗ f1s1<0.Hence, xt<x0eα1t,α1<0,and thus, limtsupxt<x0,for x0>0.

Similarly, consider the yequation,

y=f2xy1+D1x1+x+D1byy<α2y,

where α2=maxxΗf2x+f2xD1x1+x+D1by.We shall assume that maxxΗf2x+f2xD1x1+x+D1by<0.Then, yt<y0eα2twhere α2<0.Hence, limtsup yt<y0for y0>0.Thus, system (3) is proved to be uniformly bounded and dissipative, following Definition 1.

2.2. Existence of equilibrium points

According to Robinson [14], a point is termed ‘equilibrium point’ because of the forces are in equilibrium and the mass did not move. Therefore, in order to find equilibrium points of the system (3); E1,E2andE3, we equate the differential equations s,xand yto zero and solve the resulting equations simultaneously. The possible equilibrium points are as follows;

  1. E1100where no predator organism yand prey xexist. By equating system (3) to zero, we will get s=1from the equation1sf1sx=0. Then, from equation x=0,we get x=0when y=0.

  2. E2ζs1ζsD10. s=ζsis the unique solution of f1sD11=0.Let y=0,then Eq. (3) give 1sf1sx=0and f1sD11x=0.When x 0,f1sD11=0. Therefore, f1s=D1.Next, we find x.From equation 1sf1sx=0, x=1sf1s.By substituting f1s=D1and s=ζs, we get x=1ζsD1.Thus, E2is the one possible equilibrium point that consists only prey organisms and not predator organisms.

  3. E3sxydenotes the positive interior equilibrium point with s,x,y>0. sis a unique solution of 1sf1sx=0.From this, we have x=1sf1s. Next, let the equationy=0. Sincey0,f2x1+D1x1+x+D1by=1. We solve for y,and get y=2x+1f2xx1bD11f2x. Thus, the equilibrium point E3sxyis s1sf1s2x+1f2xx1bD11f2x, where x=1sf1s.

    To show the existence of E2,we let the system (3) be restricted to Rsx+as

s=1sf1sx,x=f1sD11x,E5

where s0>0and x0>0. Thus, Lemma 1 below shows the existence of non-trivial equilibrium point E2.

Suppose that a point ζsx̂exists in Rsx+such that ζs+D1x̂1=0as time, tapproaches .Then, the non-trivial equilibrium point E2ζs1ζsD10exists.

Proof. We will get two surfaces by equating system (5) to zero;

1sf1sx=0,
f1sD11x=0.

Then,

1sx=f1sandD1=f1s.

Thus, 1s=D1x, i.e., x=1sD1. By takings=ζs, we will have x=x̂=1ζsD1and satisfyingζs+D1x̂1=0. Hence, the equilibrium point E2ζs1ζsD10exists.

2.3. Stability analysis of equilibrium points E1,E2and E3

In this section, we analyse the stability of equilibrium points as it plays an important part in ordinary differential equations with their applications. As we cannot easily identify the positions of equilibrium point in applications of dynamical system, but only approximately, so the equilibrium point must be in stable state to be biologically meaningful [15]. Several definitions and theorem from Ref. [14] are stated to make us understand clearly.

[14] A fixed point xis Lyapunov stable or L-stable, provided that any solution ϕtx0stays near xfor all future time t0if the initial condition x0starts near enough to x.Specifically, a fixed point xis L-stable, provided that for any ϵ>0,there is δ>0, such that if x0x<δ,then ϕtx0x<ϵfor all t0.

Another form of stability is asymptotically stable. This is stated in the following Definition 3. below. Before going further to understand this concept, we need the definition of ω-limit set as follows.

[14] A point kis an ω-limit point of the trajectory of x0, if ϕtx0keeps coming near kas t(i.e., there is a sequence of times tj, with tjas jsuch that ϕtjx0converges to k). Certainly, if ϕtx0x0as t, then xis the only ω-limit point of x0. There can be more than one point that is an ω-limit point of x0. The set of all ω-limit points of x0is denoted by ωx0and is called the omega limit set of x0.

[14] A fixed point xis weakly asymptotically stable, if there exists δ1>0such that ωx0=xfor all x0x<δ1(i.e., ϕtx0x0as time tfor allx0x<δ1). Therefore, a fixed point is weakly asymptotically stable, if the stable manifold contains all points in a neighbourhood of the fixed point (i.e., all points are sufficiently close). A fixed point xis asymptotically stable if it is both L-stable and weakly asymptotically stable. An asymptotically stable fixed point is also called sink.

Moreover, Theorem 2 below clearly shows that if a fixed point xis hyperbolic (the real parts of all eigenvalues of the Jacobian matrix at x, i.e., DFxare nonzero), then the stability type of the fixed point for the nonlinear system is the same as for the linearized system at that fixed point.

[14] Let ẋ=Fxbe a differential equation in nvariables, with a hyperbolic fixed pointx. Suppose thal F,Fixjxand 2Fixjxkxare all continuous. Then, the stability type of the fixed point for the nonlinear system is the same as for the linearized system at that fixed point.

  1. If the real parts of all the eigenvalues of DFxare negative, then the fixed point is asymptotically stable for the nonlinear equation (i.e., if the origin of the linearized system is asymptotically stable, then xis asymptotically stable for the nonlinear equation). In this case, the basin of attraction WSxis an open set that contains some solid ball about the fixed point.

  2. If there is at least one eigenvalue of DFxhas positive real part, then the fixed point xis unstable. (For the linearized system, the fixed point can be a saddle, unstable node, unstable focus, etc.)

  3. If one of the eigenvalues of DFxhas zero real part, then the situation is more complicated. In particular, for n=2,if the fixed point is an elliptic center (eigenvalues ±βi) or one eigenvalue is zero of multiplicity one, then the linearized system does not determine the stability type of the fixed point.

Now by utilising the Theorem 2, we are going to analyse the stability of the fixed points E1,E2and E3.

(i) Stability analysis of E1.

Now, we will discuss the stability type of the equilibrium pointE1100. The Jacobian matrix of the system (3) is

J=f1sx1f1sf1sxD1f1sD1f2xyD11+x+D1byD1x1+x+D1by2f2xyD1x1+x+D1by+110f2xyD11+x+D1byD1x1+x+D1by2+f2xy1+D1x1+x+D1by0f2xD12bxy1+x+D1by2f2x1+D1x1+x+D1byf2x1+D1x1+x+D1byf2xD12bxy1+x+D1by21.

Then,

JE1=1f1100f11D110001.

As Jacobian matrix JE1above is an upper triangular, therefore the diagonal is the value of all of its eigenvalues; λ1=1,λ2=1and λ3=f11D11.If all the eigenvalues of JE1are negatives,then this leads to the following result.

If the eigenvalue λ3such that

λ3=f11D11<0,

then the trivial equilibrium point E1100is locally asymptotically stable.

(ii) Stability analysis of E2.

Next, we analyse the local stability of the system (3) that restricted to the neighbourhood of the equilibrium pointE2ζs1ζsD10. The Jacobian matrix atE2ζs1ζsD10is given as

JE2=f1ζs1ζsD11f1ζs0f1ζs1ζsD12f1ζsD11f21ζsD11+1ζs1+1ζSD100f21ζsD11+1ζs1+1ζSD11.

We can see from both Jacobian matrices JE1and JE2that there are no parameter binvolved when the predator organisms yis zero. This means that intratrophic predation does not affect the local stability and existence of E1and E2when no predator organisms involved. The characteristic equation is λ3+c1λ2+c2λ+c3=0where

c1=1D1D1ζs+1f1ζs+3D1f1ζs+f1ζsD1D1f21ζsD12ζsf1ζs3D1ζs2D12f21ζsD1+ζs2f1ζsf1ζsD1+ζsf1ζs+3D12D1ζsf1ζs+D1ζsf21ζsD1+D12ζsf21ζsD1),
c2=1D1D1ζS+12f1ζs3D12f1ζs2D1f1ζs+f21ζsD1f1ζs+2D1f21ζsD1+4ζsf1ζs+3D1ζs+4D12f21ζsD12ζs2f1ζs+2D1f1ζsf1ζsf21ζsD12ζsf1ζs3D122D1f1ζsf21ζsD1+ζsf1ζsf21ζsD1+2D1f1ζsf21ζsD1+2D1ζsf1ζs2ζsf1ζsf21ζsD12D1ζsf21ζsD1+ζs2f1ζsf21ζsD12D12ζsf21ζsD13D1ζsf1ζsf21ζsD1+D1ζs2f1ζsf21ζsD1+D1ζsf1ζsf21ζsD1,
c3=1D1D1ζs+1f1ζsD1f1ζsD1f1ζs+f1ζsf21ζsD1+D1f21ζsD1+2ζsf1ζs+D1ζsζs2f1ζs+f1ζsD1f1ζsf21ζsD1ζsf1ζsD122D1f1ζsf21ζsD1+ζsf1ζsf21ζsD1+2D1f1ζsf21ζsD1+D1ζsf1ζs2ζsf1ζsf21ζsD1D1ζsf21ζsD1+ζs2f1ζsf21ζsD1D12ζsf21ζsD13D1ζsf1ζsf21ζsD1+D1ζs2f1ζsf21ζsD1+D1ζsf1ζsf21ζsD1.

Then, by using MATLAB R2015b, we get the eigenvalues of JE2which are λ4=1, λ5=f1ζsζs1+f1ζsD1D1and λ6=f21ζsD11+2D1ζsD1ζs1D1ζs+1where D1ζs+10.We summarise the above discussion in the following theorems.

Suppose the assumptions such as

f1ζsζs1+f1ζsD1D1<0and f21ζsD11+2D1ζsD1ζs1D1ζs+1<0

are satisfied, then the equilibrium point, E2ζs1ζsD10is locally asymptotically stable.

If

f1ζsζs1+f1ζsD1D1<0and f21ζsD11+2D1ζsD1ζs1D1ζs+1<0,

then the equilibrium point E2ζs1ζsD10is a hyperbolic saddle and is repels locally in y-direction. Particularly, the dimension of the stable manifold WsE2and the unstable manifold WuE2are given by

dimWsE2=2and dimWuE2=1.

Proof. The results follow from inspections of the eigenvalues of the matrixJE2and Theorem 2 (see [16, 18]).

[17] The flow Fwill be called uniformly persistent if there exists ε0>0such that for all xE0, limtdπxtEε0.

The following theorem shows the existence of the equilibrium point E2using persistence analysis.

Assume that

  1. Lemma 1 being holds,

  2. E2is a unique hyperbolic saddle point in Rsxy+and repels locally in y-direction (as in Theorem 5),

  3. no existence of periodic orbits in the planes of Rsxy+.

Then

limt+infst>k,limt+infxt>k,limt+infyt>k,

where k>0.

Particularly, the system (3) exhibits uniform persistence and the equilibrium point E2ζs1ζsD10exists.

Proof. The result follows from the Definition 4, which defines uniform persistence by Butler et al. [17] and Nani and Freedman [18].

(iii) Stability analysis of E3sxy.

The Jacobian matrix at E3is

JE3=f1sx1f1sf1sxD1f1sD1f2xyD11+x+D1byD1x1+x+D1by2f2xyD1x1+x+D1by+110f2xyD11+x+D1byD1x1+x+D1by2+f2xyD1x1+x+D1by+10f2xD12bxy1+x+D1by2f2xD1x1+x+D1by+1f2xD1x1+x+D1by+1f2xD12bxy1+x+D1by21.

The eigenvalues of JE3are resulted from the characteristic equation below

λ3+c4λ2+c5λ+c6=0,E6

where

c4=f1sx+Q+R+VZ+3f1sD1,
c5=f1sxQ+R+VZ+2+f1sZR2D1+2Q+2R+2V2Z+3,
c6=f1sxQ+R+VZ+1+f1sZR1D1+Q+R+VZ+1,

and

Q=f2xyD11+x+D1byD1x1+x+D1by2,
R=f2xD12bxy1+x+D1by2,
V=f2xyD1x1+x+D1by+1,
Z=f2xD1x1+x+D1by+1.

From (6), we obtain the eigenvalues; λ7=1,and λ8, λ9are

λ8=12D1D1f1sx+Q+R+VZ+2f1s+D1Q2+2QR+2QV2QZ2Qf1sf1sx2Q+2R+2V2Zf1sx+R2+2RV2RZ+V22VZ+Z2+f1s2D1R2D1V2D1Zf1sx+f1s2+f1sD1x,
λ9=12D1D1f1sx+Q+R+VZ+2f1sD1Q2+2QR+2QV2QZ2Qf1sf1sx2Q+2R+2V2Zf1sx+R2+2RV2RZ+V22VZ+Z2+f1s2D1R2D1V2D1Zf1sx+f1s2+f1sD1x.

These results lead to the following theorem.

Suppose that λ8<0and λ9<0, then the equilibrium point E3is locally asymptotically stable.

2.4. Global stability and uniform persistence analysis

Now we will analyse the global stability of the equilibrium point E2ζs1ζsD10and the existence of positive interior equilibrium point E3sxy. For global stability analysis, we use the similar technique as in [3, 18].

(i) Global stability analysis of E2ζs1ζsD10

Consider system (3) restricted to Rsx+as in (5). Let Nbe a neighbourhood of equilibrium point E2ζs1ζsD10in Rsx+. To analyse the global stability of the equilibrium point E2, a suitable Lyapunov function L=12n1sŝ2+n2xx̂2is used, where ŝand x̂denote the components of E2ŝx̂, i.e., ŝ=ζsand x̂=1ζsD1, while n1and n2are positive constants. Note that Lis a positive definite function with respect to E2in Rsx+and a Lyapunov function for system (5) in N. By differentiating Lwith respect to time t,we get

dLdt=n1sŝs+n2xx̂x,E7

where x̂=1ŝD1. From (5), we have

1=ŝ+x̂f1ŝandD1=f1ŝ.

Hence (7) can be written as

L=n1sŝŝ+x̂f1ŝsf1sx+n2xx̂f1sf1ŝ1x=n1sŝ2+n1sŝf1ŝx̂f1sx+n2xx̂xf1sf1ŝ1=n11sŝ2+12n12sŝxx̂+12n21sŝxx̂+n22xx̂2,

where

n11=n1<0,
n12=n21=n1f1ŝx̂f1sxxx̂,
n22=n2xf1sf1ŝxx̂f1ŝ.

Clearly that Lcan be written as L=XTNX,which Tdenotes the transpose and the matrix Nis particularly a real, symmetric 2×2matrix, where Xand Ncan be represented by

X=v1v2=sŝxx̂ and N=n1112n1212n21n22.

Thus, it leads to the following theorem.

The equilibrium point E2is global asymptotically stable with respect to solution trajectories are initiated from intRsx+if the assumptions n22<0and detN>0are satisfied.

Proof. By using the Frobenius Theorem in ([18], Lemma 6.2), we can see that n22and detNare the leading principal minors of the matrix N.It is shown that matrix Nis negative definite if

n22<0, and detN=detn1112n1212n21n22>0.

This completes the proof of the theorem.

(ii) Existence of positive interior equilibrium point E3

In this subsection, we present some results of persistence analysis, including uniform persistence and state the necessary conditions for the existence of positive equilibrium point E3. The following lemma from Ref. [19] is applied to obtain the persistence results.

[19] Let Gbe an isolated hyperbolic equilibrium point in the omega limit set, ωXof the orbit OX.Then either ωX=Gor there exist points J+and Jin ωXwith J+WsGand JWuG, where WsGand WuGdenote stable and unstable manifolds of G.

Assume that

  1. E2ζs1ζsD10is a hyperbolic saddle point and is repels locally in ydirection (as in Theorem 5),

  2. system (3) is dissipative and all solutions with initial values in Rsxy+are uniformly bounded and attracted into region Η(as in Theorem 1), and

  3. equilibrium point E2ζs1ζsD10is globally asymptotically stable with respect to Rsx+.

Then, the system (3) is uniformly persistence.

Proof. This proof strictly depends on Lemma 2. Suppose Ηis the region as stated in Theorem 1. It showed that region Ηis positive invariant set and any solutions of system (3) emanating at a point in R+3is uniformly bounded. Despites that, the only compact invariant set on R+3is E2ζs1ζsD10. Let P=E3sxybelongs to the interior of R+3, i.e., PintR+3. We shall show that there is no points JiR+3where R+3belongs to ωP, the omega limit set of P.Now, we prove that the equilibrium point E2ζs1ζsD10ωP.Assume that E2ωPis true. Then, there is a point J1+WsE2\E2such that J1+ωPby Lemma 2. But WsE2R+3\E2is empty which is a contradiction for the positive invariance property of ΗR+3.Thus, the equilibrium point E2is not in the omega limit set of P; E2ωP.Next, we shall show R+3ωP=.Assume that R+3ωP, and let JR+3and JωP.Then the closure of the orbit of the point J, i.e., clOJmust either contains E2or unbounded. This is a contradiction, and hence it is proved that R+3ωP=.We deduce that if E2ζs1ζsD10is unstable, then, for stable manifold WsE2;WsE2intR+3=, and for unstable manifold WuE2;WuE2intR+3. Therefore, the result of uniform persistence follows since the omega limit set of P, ωPmust be in intR+3.This completes the proof.

Global stability of equilibrium point E2ζs1ζsD10with respect to Rsx+indicates that the boundary flow is isolated and a cyclic with respect to regionΗ. Thus, the system (3) undergoes uniform persistence and implies that a positive interior equilibrium point E3sxyexists (see [20]).

2.5. Hopf bifurcation

In this section, we investigate Hopf bifurcation on the system (3) with a bifurcation real parameter, σ.Particularly, σis selected in such a way that the growth rate function f2is a function of xand σ.Therefore, system (3) takes of the form

s=1sf1sx,x=f1sD11xf2xσy1+D1x1+x+D1by,y=f2xσy1+D1x1+x+D1byy,E8

where s0>0,x0>0,y0>0.Next, we do linearization on the system (8). First, let

S=sh1;X=xh2;Y=yh3;s=S+h1;x=X+h2;y=Y+h3;

where h1h2h3is the non-trivial equilibrium point. Then, we obtain the following differential equations

S=1Sh1f1S+h1X+h2,X=f1S+h1D11X+h2f2X+h2σY+h31+D1X+h21+X+h2+D1bY+h3,Y=f2X+h2σY+h31+D1X+h21+X+h2+D1bY+h3Y+h3.E9

The Jacobian matrix of system (8) is given by

Jσ=f1sx1f1sf1sxD1f1sD1f2xσyD11+x+D1byD1x1+x+D1by2f2xσyD1x1+x+D1by+110f2xσyD11+x+D1byD1x1+x+D1by2+f2xσy1+D1x1+x+D1by0f2xσD12bxy1+x+D1by2D1x1+x+D1by+1f2xσ1+D1x1+x+D1byD12bxy1+x+D1by21.

Thus, the Jacobian matrix of system (8) about E2ζs1ζsD10σis

JσE2=f1ζs1ζsD11f1ζs0f1ζs1ζsD12f1ζsD11f21ζsD1σ1+1ζs1+1ζsD100f21ζsD1σ1+1ζs1+1ζsD11.E10

The characteristic equation of (10) is given as

λ3+c1λ2+c2λ+c3=0,E11

where

c1=1D1D1ζs+1f1ζs1+D12ζs+ζs2D1ζs+f1ζsζsD11+D1f21ζsD1σζs+D1ζs2D11+3D11+D1ζs,
c2=1D1D1ζs+1f1ζs2ζs2+2+2D14ζs2D1ζsf21ζsD1σ1+2D12ζs+ζs23D1ζs+D1ζs2f1ζs2+2D12ζs+f21ζsD1σζs2D11+D1ζsf21ζsD1σ2D1+4D122D1ζs2D12ζs+3D1ζs3D123D1,
c3=1D1D1ζs+1f1ζs12ζs+D1+ζs2D1ζsf21ζsD1σ12ζs+ζs2D1+D1ζs2f1ζsf21ζsD1σζs12D1+D1ζs+1ζsf21ζsD1σD1D1ζs+D1ζsD1D12.

We applied Routh Hurwitz criterion (see [21, 22]) onto the characteristic equation (11) and obtain the matrices

M1=c1;M2=c11c3c2;M3=c110c3c2c100c3.

Thus, the characteristic equation (11) has all negative real parts of λif and only if

c1>0c3>0c1c2c3>0.E12

When the assumptions on functional responses f1s=f1ζsand f2x=f21ζsD1σas in (2) for the system (3), together with the hypothesesH1until H3;

H1:D1ζs+1>0,
H2:f1ζs1+D12ζs+ζs2D1ζs+f1ζsζsD11+D1f21ζsD1σζs+D1ζs2D11+3D11+D1ζs>0,
H3:f1ζs12ζs+D1+ζs2D1ζsf21ζsD1σ12ζs+ζs2D1+D1ζs2f1ζsf21ζsD1σζs12D1+D1ζs+1ζsf21ζsD1σD1D1ζs+D1ζsD1D12>0,

hold, we will have c1>0and c3>0.We shall obtain two pure imaginary roots for the characteristic equation (11) if and only if c1c2=c3for some values of σ,say, σ1.

Since at σ=σ1,there exists an interval σ1εσ1+εcontaining σ1for some ε>0such that σσ1εσ1+ε. Thus, for σσ1εσ1+ε,the characteristic equation (11) cannot have positive real roots. For σ=σ1,we acquire (see [3, 10])

λ2+c2λ+c1=0,E13

that consist of three roots; λ1=c2i,λ2=c2i,and λ3=c1.As for σσ1εσ1+ε,all the roots are in general of the form;

λ1σ=ασ+βσi,
λ2σ=ασβσi,
λ3σ=c1σ.

In order to apply the Hopf bifurcation theorem as stated in [9, 23] towards the system (8), we must verify the transversality condition

Redλjσ=σ10,j=1,2,3.E14

By substituting λ1σ=ασ+βσiand λ2σ=ασβσiinto characteristic equation (11) and calculating the implicit derivative, we obtain the following equations

KσασLσβσ+Mσ=0,Lσασ+Kσβσ+Nσ=0,E15

where

Kσ=3α2σ+2c1σασ+c2σ3β2σ;
Mσ=α2σc1σ+c2σασc1σβ2σ+c3σ;
Lσ=6ασβσ+2c1σβσ;
Nσ=2ασβσc1σ+c2σβσ.

Since Kσ1Mσ1+Lσ1Nσ10,we have

Redλjσ=σ1=Kσ1Mσ1+Lσ1Nσ1K2σ1+L2σ10,j=1,2,3,

and λ3σ1=c1σ10.We conclude the details above in the following theorem.

Suppose that the equilibrium point E2ζs1ζsD10σexists and those assumptions similar as in (2) for the system (3) together with hypothesis H1until H3hold. Then the system (8) undergoes Hopf bifurcation in the first octant, which leads to a family of periodic solutions bifurcating from E2ζs1ζsD10σfor some values of σin the neighbourhood of σ1.

Next, we determine the Hopf bifurcation at the equilibrium point E3sxyσ.The Jacobian matrix of the system (8) about the equilibrium point E3is given by

JσE3=f1sx1f1sf1sxD1f1sD1f2xσyD11+x+D1byD1x1+x+D1by2f2xσyD1x1+x+D1by+110f2xσyD11+x+D1byD1x1+x+D1by2+f2xσyD1x1+x+D1by+10f2xσD12bxy1+x+D1by2D1x1+x+D1by1f2xσD1x1+x+D1byD12bxy1+x+D1by2+11.

Hence, the characteristic equation for the Jacobian matrix JσE3is

λ3+c4λ2+c5λ+c6=0,E16

where

c4=f1sx+Q+R+VZ+3f1sD1,
c5=f1sxQ+R+VZ+2+f1sZR2D1+2Q+2R+2V2Z+3,
c6=f1sxQ+R+VZ+1+f1sZR1D1+Q+R+VZ+1,

and

Q=f2xσyD11+x+D1byD1x1+x+D1by2,
R=f2xσD12bxy1+x+D1by2,
V=f2xσyD1x1+x+D1by+1,
Z=f2xσD1x1+x+D1by+1.

By applying the Routh-Hurwitz criterion (see [21, 22]) towards the characteristic equation (16), we obtain the following Hurwitz matrices

M4=c4;M5=c41c6c5;M6=c410c6c5c400c6.

Thus, the characteristic equation (16) has all negative real parts of λif and only if

c4>0c6>0c4c5c6>0.E17

Suppose the assumptions of functional response f1sand f2xσsimilar as in (2) for the system (3), together with the hypotheses H4until H6;

H4:f1sx+Q+R+V+3>f1sD1+Z,
H5:Q+R+V+1>Z,
H6:Z>R+1,

hold, then clearly c4>0and c6>0. In particular, we shall have two pure imaginary roots for the characteristic equation (16) if and only if c4c5=c6for some values of σ,say, σ2.Since at σ=σ2,there exists an interval σ2εσ2+εcontaining σ2for some ε>0. Then, for σσ2εσ2+ε, the characteristic Eq. (16) cannot have positive real roots. For σ=σ2,we get (see [3, 10])

λ2+c5λ+c4=0,E18

Consisting of three roots; λ1=c5i,λ2=c5i,λ3=c4.As for σσ2εσ2+ε,all roots are in general of the form;

λ1σ=ασ+βσi,
λ2σ=ασβσi,
λ3σ=c4σ,

To establish Hopf bifurcation towards system (8), we must show that

Redλjσ=σ20,j=1,2,3.E19

By substituting λ1σ=ασ+βσiand λ2σ=ασβσiinto characteristic equation (18) and calculating the implicit derivative, we get the following equations;

A1σασA2σβσ+B1σ=0,A2σασ+A1σβσ+B2σ=0,E20

where

A1σ=3α2σ+2c1σασ+c2σ3β2σ;
A2σ=6ασβσ+2c1σβσ;
B1σ=α2σc1σ+c2σασc1σβ2σ+c3σ;
B2σ=2ασβσc1σ+c2σβσ.

Since A1σ2B1σ2+A2σ2B2σ20,we have

Redλjσ=σ2=A1σ2B1σ2+A2σ2B2σ2A1σ22+A2σ220,j=1,2,3,

and λ3σ2=c4σ20.We summarise the discussion above in the following theorem.

Suppose that the equilibrium point

E3sxyσ=E3s1sf1s2x+1f2xx1bD11f2xσ

exists and those assumptions similar as in (2) for the system (3) together with hypothesesH4until H6hold. Then the system (8) undergoes Hopf bifurcation in the first octant, which leads to a family of periodic solutions bifurcating from E3sxyσfor some values of σin the neighbourhood of σ2.

2.6. Discussion

We have proposed and analysed a simple nutrient-predator-prey model in a chemostat with intratrophic predation. This system consisted of the nutrients, prey organisms xand predator organisms y. We conclude that intratrophic predation denoted as bdoes not affected the local stability and existence of E1and E2when no predator organisms involved. Next, we have shown that the washout equilibrium point E1100(no prey and predator organisms present) is locally asymptotically stable if λ3=f11D11<0holds. For stability of E2and E3of the system (3), some sufficient criteria or conditions are derived and satisfied. The points E2and E3are said to be asymptotically stable if all of their eigenvalues are less than zero.

In particular, we investigated the global stability analysis for E2.A suitable Lyapunov function Lis defined and E2is globally asymptotically stable if and only if the conditions in Theorem 8 holds. Global stability of E2indicates that predator organisms ymight be washout in the chemostat despites the initial prey and predator organisms’ density levels. Next, in the study of the existence of positive interior equilibrium point E3,we presented some results regarding uniform persistence analysis. It has shown that the system (3) is uniformly persistence and thus, the positive interior equilibrium point E3exists.

In the analysis of occurrence of Hopf bifurcation, Hopf bifurcation theorems in Hassard et al. [9] are applied. We have shown that the system (8) undergoes Hopf bifurcation in the first octant, which leads to a family of periodic solutions bifurcating from E2ζs1ζsD10σand E3sxyσfor some values of σin the neighbourhoods of σ1and σ2, respectively. The method used to obtain these results is similar to the method used in [3, 10].

© 2017 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Zabidin Salleh and Liyana Abd Rahim (December 20th 2017). Stability and Hopf Bifurcation Analysis of a Simple Nutrient- Prey-Predator Model with Intratrophic Predation in Chemostat, Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals, Ricardo López-Ruiz, IntechOpen, DOI: 10.5772/intechopen.71624. Available from:

chapter statistics

272total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Sensitivity Analysis: A Useful Tool for Bifurcation Analysis

By Raheem Gul and Stefan Bernhard

Related Book

First chapter

Some Commonly Used Speech Feature Extraction Algorithms

By Sabur Ajibola Alim and Nahrul Khair Alang Rashid

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us