Open access peer-reviewed chapter

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

Written By

Zabidin Salleh and Liyana Abd Rahim

Submitted: 08 April 2017 Reviewed: 11 October 2017 Published: 20 December 2017

DOI: 10.5772/intechopen.71624

From the Edited Volume

Complexity in Biological and Physical Systems - Bifurcations, Solitons and Fractals

Edited by Ricardo López-Ruiz

Chapter metrics overview

1,227 Chapter Downloads

View Full Metrics

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 b and term 1 + D 1 x 1 + x + D 1 by are added to differential equations x and y in 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.

Advertisement

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 s 0 = s 0 > 0 , x 0 = x 0 > 0 , y 0 = y 0 > 0 , D i = 1 ; i = 0 , 1 , 2 , β j 1 ; j = 1 , 2 and 0 b 1 . This leads to the following assumptions on functional response f n , n = 1 , 2 :

i f n : R + R + and f n , n = 1 , 2 are continuously differential equations , ii f n 0 = 0 , iii f 1 s > 0 , f 2 x > 0 for all s , 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, t with the density of prey organisms.

By referring to system (1), s represents the concentration of nutrient at time t while s 0 is the input nutrient concentration. The variables x and y represent the concentration of prey and predator at time t , respectively. Parameters β 1 and β 2 denote the yield constants, D 0 is the washout rate of the chemostat while D 1 and D 2 are the removal rates of x and y , respectively. Removal rate is the sum of washout rate and death rate θ i , i.e., D i = D 0 + θ i , i = 1 , 2 . f 1 s and f 2 x denote the specific growth rate of prey x and predator y , respectively, while b is 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 ¯ = s s 0 , x ¯ = x D 0 β 1 s 0 , y ¯ = y D 0 D 1 β 1 β 2 s 0 ,
f 1 s = f 1 ~ s 0 s ¯ D 0 , f 2 x = f 2 ~ D 0 β 1 s 0 x ¯ D 0 , b = b ¯ β 2 ,
t ¯ = D i t , i = 0 , 1 , D 2 = D 0 β 1 s 0 .

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

s = 1 s f 1 s x , x = f 1 s D 1 1 x f 2 x y 1 + D 1 x 1 + x + D 1 by , y = f 2 x y 1 + D 1 x 1 + x + D 1 by y , E3

where

s 0 > 0 , x 0 > 0 , y 0 > 0 and D 1 > 0 . E4

Now, the variables are non-dimensional and the discussion is in R + 3 = s x y : s > 0 x > 0 y > 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 in R + 3 . Firstly, we consider dissipativity of the system (3).

[13]. A system with differential equations x = f x is defined to be dissipative if there exists a bounded subset Γ of R 3 , such that there is a time t 0 for any x 0 R 3 which depends on x 0 and Γ so that the solution of the system ϕ t x 0 Γ for t t 0 .

Let H be the region

H = s x y R + 3 : 1 P max q s + x + y 1 P min + q ,

where P max = max 1 f 1 s + f 1 s D 1 1 , P min = min 1 f 1 s + f 1 s D 1 1 and q is a positive constant.

Then,

  1. Η is positively invariant.

  2. All nonnegative solutions of (3) with initial values in R + 3 are 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 1 P max q s + x + y 1 P min + q . Let us define

P max = max 1 f 1 s + f 1 s D 1 1   and   P min = min 1 f 1 s + f 1 s D 1 1 .

By adding those differential equations s , x and y in (3), we shall get

s + x + y = 1 s + f 1 s + f 1 s D 1 1 x + y .

Hence, this leads to

1 P max s + x + y 1 s + f 1 s + f 1 s D 1 1 x + y 1 P min s + x + y .

Thus, solving the above inequality gives

1 P max q s + x + y 1 P min + q ,

as q is the positive constant

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

s = 1 s f 1 s x < 1 s .

Then s t < 1 e t 1 s 0 for 0 < s 0 < 1 , and hence lim t sup s t < 1 for 0 < s 0 < 1 .

Now, consider the x equation;

x = f 1 s D 1 1 x f 2 x y 1 + D 1 x 1 + x + D 1 by < f 1 s 1 x < α 1 x ,

where α 1 = max s Η   f 1 s 1 . We assume that α 1 = max s Η   f 1 s 1 < 0 . Hence, x t < x 0 e α 1 t , α 1 < 0 , and thus, lim t sup x t < x 0 , for x 0 > 0 .

Similarly, consider the y equation,

y = f 2 x y 1 + D 1 x 1 + x + D 1 by y < α 2 y ,

where α 2 = max x Η f 2 x + f 2 x D 1 x 1 + x + D 1 by . We shall assume that max x Η f 2 x + f 2 x D 1 x 1 + x + D 1 by < 0 . Then, y t < y 0 e α 2 t where α 2 < 0 . Hence, lim t sup   y t < y 0 for y 0 > 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); E 1 , E 2 and E 3 , we equate the differential equations s , x and y to zero and solve the resulting equations simultaneously. The possible equilibrium points are as follows;

  1. E 1 1 0 0 where no predator organism y and prey x exist. By equating system (3) to zero, we will get s = 1 from the equation 1 s f 1 s x = 0 . Then, from equation x = 0 , we get x = 0 when y = 0 .

  2. E 2 ζ s 1 ζ s D 1 0 . s = ζ s is the unique solution of f 1 s D 1 1 = 0 . Let y = 0 , then Eq. (3) give 1 s f 1 s x = 0 and f 1 s D 1 1 x = 0 . When x 0 , f 1 s D 1 1 = 0 . Therefore, f 1 s = D 1 . Next, we find x . From equation 1 s f 1 s x = 0 , x = 1 s f 1 s . By substituting f 1 s = D 1 and s = ζ s , we get x = 1 ζ s D 1 . Thus, E 2 is the one possible equilibrium point that consists only prey organisms and not predator organisms.

  3. E 3 s x y denotes the positive interior equilibrium point with s , x , y > 0 . s is a unique solution of 1 s f 1 s x = 0 . From this, we have x = 1 s f 1 s . Next, let the equation y = 0 . Since y 0 , f 2 x 1 + D 1 x 1 + x + D 1 by = 1 . We solve for y , and get y = 2 x + 1 f 2 x x 1 b D 1 1 f 2 x . Thus, the equilibrium point E 3 s x y is s 1 s f 1 s 2 x + 1 f 2 x x 1 b D 1 1 f 2 x , where x = 1 s f 1 s .

    To show the existence of E 2 , we let the system (3) be restricted to R sx + as

s = 1 s f 1 s x , x = f 1 s D 1 1 x , E5

where s 0 > 0 and x 0 > 0 . Thus, Lemma 1 below shows the existence of non-trivial equilibrium point E 2 .

Suppose that a point ζ s x ̂ exists in R sx + such that ζ s + D 1 x ̂ 1 = 0 as time, t approaches . Then, the non-trivial equilibrium point E 2 ζ s 1 ζ s D 1 0 exists.

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

1 s f 1 s x = 0 ,
f 1 s D 1 1 x = 0 .

Then,

1 s x = f 1 s and D 1 = f 1 s .

Thus, 1 s = D 1 x , i.e., x = 1 s D 1 . By taking s = ζ s , we will have x = x ̂ = 1 ζ s D 1 and satisfying ζ s + D 1 x ̂ 1 = 0 . Hence, the equilibrium point E 2 ζ s 1 ζ s D 1 0 exists.

2.3. Stability analysis of equilibrium points E 1 , E 2 and E 3

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 x is Lyapunov stable or L-stable, provided that any solution ϕ t x 0 stays near x for all future time t 0 if the initial condition x 0 starts near enough to x . Specifically, a fixed point x is L-stable, provided that for any ϵ > 0 , there is δ > 0 , such that if x 0 x < δ , then ϕ t x 0 x < ϵ for all t 0 .

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 k is an ω -limit point of the trajectory of x 0 , if ϕ t x 0 keeps coming near k as t (i.e., there is a sequence of times t j , with t j as j such that ϕ t j x 0 converges to k ). Certainly, if ϕ t x 0 x 0 as t , then x is the only ω -limit point of x 0 . There can be more than one point that is an ω -limit point of x 0 . The set of all ω -limit points of x 0 is denoted by ω x 0 and is called the omega limit set of x 0 .

[14] A fixed point x is weakly asymptotically stable, if there exists δ 1 > 0 such that ω x 0 = x for all x 0 x < δ 1 (i.e., ϕ t x 0 x 0 as time t for all x 0 x < δ 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 x is 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 x is hyperbolic (the real parts of all eigenvalues of the Jacobian matrix at x , i.e., D F x are 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 x ̇ = F x be a differential equation in n variables, with a hyperbolic fixed point x . Suppose thal F , F i x j x and 2 F i x j x k x are 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 D F x are 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 x is asymptotically stable for the nonlinear equation). In this case, the basin of attraction W S x is an open set that contains some solid ball about the fixed point.

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

  3. If one of the eigenvalues of D F x has 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 E 1 , E 2 and E 3 .

(i) Stability analysis of E 1 .

Now, we will discuss the stability type of the equilibrium point E 1 1 0 0 . The Jacobian matrix of the system (3) is

J = f 1 s x 1 f 1 s f 1 s x D 1 f 1 s D 1 f 2 x y D 1 1 + x + D 1 by D 1 x 1 + x + D 1 by 2 f 2 x y D 1 x 1 + x + D 1 by + 1 1 0 f 2 x y D 1 1 + x + D 1 by D 1 x 1 + x + D 1 by 2 + f 2 x y 1 + D 1 x 1 + x + D 1 by 0 f 2 x D 1 2 bxy 1 + x + D 1 by 2 f 2 x 1 + D 1 x 1 + x + D 1 by f 2 x 1 + D 1 x 1 + x + D 1 by f 2 x D 1 2 bxy 1 + x + D 1 by 2 1 .

Then,

J E 1 = 1 f 1 1 0 0 f 1 1 D 1 1 0 0 0 1 .

As Jacobian matrix J E 1 above is an upper triangular, therefore the diagonal is the value of all of its eigenvalues; λ 1 = 1 , λ 2 = 1 and λ 3 = f 1 1 D 1 1 . If all the eigenvalues of J E 1 are negatives , then this leads to the following result.

If the eigenvalue λ 3 such that

λ 3 = f 1 1 D 1 1 < 0 ,

then the trivial equilibrium point E 1 1 0 0 is locally asymptotically stable.

(ii) Stability analysis of E 2 .

Next, we analyse the local stability of the system (3) that restricted to the neighbourhood of the equilibrium point E 2 ζ s 1 ζ s D 1 0 . The Jacobian matrix at E 2 ζ s 1 ζ s D 1 0 is given as

J E 2 = f 1 ζ s 1 ζ s D 1 1 f 1 ζ s 0 f 1 ζ s 1 ζ s D 1 2 f 1 ζ s D 1 1 f 2 1 ζ s D 1 1 + 1 ζ s 1 + 1 ζ S D 1 0 0 f 2 1 ζ s D 1 1 + 1 ζ s 1 + 1 ζ S D 1 1 .

We can see from both Jacobian matrices J E 1 and J E 2 that there are no parameter b involved when the predator organisms y is zero. This means that intratrophic predation does not affect the local stability and existence of E 1 and E 2 when no predator organisms involved. The characteristic equation is λ 3 + c 1 λ 2 + c 2 λ + c 3 = 0 where

c 1 = 1 D 1 D 1 ζ s + 1 f 1 ζ s + 3 D 1 f 1 ζ s + f 1 ζ s D 1 D 1 f 2 1 ζ s D 1 2 ζ s f 1 ζ s 3 D 1 ζ s 2 D 1 2 f 2 1 ζ s D 1 + ζ s 2 f 1 ζ s f 1 ζ s D 1 + ζ s f 1 ζ s + 3 D 1 2 D 1 ζ s f 1 ζ s + D 1 ζ s f 2 1 ζ s D 1 + D 1 2 ζ s f 2 1 ζ s D 1 ) ,
c 2 = 1 D 1 D 1 ζ S + 1 2 f 1 ζ s 3 D 1 2 f 1 ζ s 2 D 1 f 1 ζ s + f 2 1 ζ s D 1 f 1 ζ s + 2 D 1 f 2 1 ζ s D 1 + 4 ζ s f 1 ζ s + 3 D 1 ζ s + 4 D 1 2 f 2 1 ζ s D 1 2 ζ s 2 f 1 ζ s + 2 D 1 f 1 ζ s f 1 ζ s f 2 1 ζ s D 1 2 ζ s f 1 ζ s 3 D 1 2 2 D 1 f 1 ζ s f 2 1 ζ s D 1 + ζ s f 1 ζ s f 2 1 ζ s D 1 + 2 D 1 f 1 ζ s f 2 1 ζ s D 1 + 2 D 1 ζ s f 1 ζ s 2 ζ s f 1 ζ s f 2 1 ζ s D 1 2 D 1 ζ s f 2 1 ζ s D 1 + ζ s 2 f 1 ζ s f 2 1 ζ s D 1 2 D 1 2 ζ s f 2 1 ζ s D 1 3 D 1 ζ s f 1 ζ s f 2 1 ζ s D 1 + D 1 ζ s 2 f 1 ζ s f 2 1 ζ s D 1 + D 1 ζ s f 1 ζ s f 2 1 ζ s D 1 ,
c 3 = 1 D 1 D 1 ζ s + 1 f 1 ζ s D 1 f 1 ζ s D 1 f 1 ζ s + f 1 ζ s f 2 1 ζ s D 1 + D 1 f 2 1 ζ s D 1 + 2 ζ s f 1 ζ s + D 1 ζ s ζ s 2 f 1 ζ s + f 1 ζ s D 1 f 1 ζ s f 2 1 ζ s D 1 ζ s f 1 ζ s D 1 2 2 D 1 f 1 ζ s f 2 1 ζ s D 1 + ζ s f 1 ζ s f 2 1 ζ s D 1 + 2 D 1 f 1 ζ s f 2 1 ζ s D 1 + D 1 ζ s f 1 ζ s 2 ζ s f 1 ζ s f 2 1 ζ s D 1 D 1 ζ s f 2 1 ζ s D 1 + ζ s 2 f 1 ζ s f 2 1 ζ s D 1 D 1 2 ζ s f 2 1 ζ s D 1 3 D 1 ζ s f 1 ζ s f 2 1 ζ s D 1 + D 1 ζ s 2 f 1 ζ s f 2 1 ζ s D 1 + D 1 ζ s f 1 ζ s f 2 1 ζ s D 1 .

Then, by using MATLAB R2015b, we get the eigenvalues of J E 2 which are λ 4 = 1 , λ 5 = f 1 ζ s ζ s 1 + f 1 ζ s D 1 D 1 and λ 6 = f 2 1 ζ s D 1 1 + 2 D 1 ζ s D 1 ζ s 1 D 1 ζ s + 1 where D 1 ζ s + 1 0 . We summarise the above discussion in the following theorems.

Suppose the assumptions such as

f 1 ζ s ζ s 1 + f 1 ζ s D 1 D 1 < 0 and f 2 1 ζ s D 1 1 + 2 D 1 ζ s D 1 ζ s 1 D 1 ζ s + 1 < 0

are satisfied, then the equilibrium point, E 2 ζ s 1 ζ s D 1 0 is locally asymptotically stable.

If

f 1 ζ s ζ s 1 + f 1 ζ s D 1 D 1 < 0 and f 2 1 ζ s D 1 1 + 2 D 1 ζ s D 1 ζ s 1 D 1 ζ s + 1 < 0 ,

then the equilibrium point E 2 ζ s 1 ζ s D 1 0 is a hyperbolic saddle and is repels locally in y -direction. Particularly, the dimension of the stable manifold W s E 2 and the unstable manifold W u E 2 are given by

dim W s E 2 = 2 and dim W u E 2 = 1 .

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

[17] The flow F will be called uniformly persistent if there exists ε 0 > 0 such that for all x E 0 , lim t d π x t E ε 0 .

The following theorem shows the existence of the equilibrium point E 2 using persistence analysis.

Assume that

  1. Lemma 1 being holds,

  2. E 2 is a unique hyperbolic saddle point in R sxy + and repels locally in y -direction (as in Theorem 5),

  3. no existence of periodic orbits in the planes of R sxy + .

Then

lim t + inf s t > k , lim t + inf x t > k , lim t + inf y t > k ,

where k > 0 .

Particularly, the system (3) exhibits uniform persistence and the equilibrium point E 2 ζ s 1 ζ s D 1 0 exists.

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 E 3 s x y .

The Jacobian matrix at E 3 is

J E 3 = f 1 s x 1 f 1 s f 1 s x D 1 f 1 s D 1 f 2 x y D 1 1 + x + D 1 b y D 1 x 1 + x + D 1 b y 2 f 2 x y D 1 x 1 + x + D 1 b y + 1 1 0 f 2 x y D 1 1 + x + D 1 b y D 1 x 1 + x + D 1 b y 2 + f 2 x y D 1 x 1 + x + D 1 b y + 1 0 f 2 x D 1 2 b x y 1 + x + D 1 b y 2 f 2 x D 1 x 1 + x + D 1 b y + 1 f 2 x D 1 x 1 + x + D 1 b y + 1 f 2 x D 1 2 b x y 1 + x + D 1 b y 2 1 .

The eigenvalues of J E 3 are resulted from the characteristic equation below

λ 3 + c 4 λ 2 + c 5 λ + c 6 = 0 , E6

where

c 4 = f 1 s x + Q + R + V Z + 3 f 1 s D 1 ,
c 5 = f 1 s x Q + R + V Z + 2 + f 1 s Z R 2 D 1 + 2 Q + 2 R + 2 V 2 Z + 3 ,
c 6 = f 1 s x Q + R + V Z + 1 + f 1 s Z R 1 D 1 + Q + R + V Z + 1 ,

and

Q = f 2 x y D 1 1 + x + D 1 b y D 1 x 1 + x + D 1 b y 2 ,
R = f 2 x D 1 2 b x y 1 + x + D 1 b y 2 ,
V = f 2 x y D 1 x 1 + x + D 1 b y + 1 ,
Z = f 2 x D 1 x 1 + x + D 1 b y + 1 .

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

λ 8 = 1 2 D 1 D 1 f 1 s x + Q + R + V Z + 2 f 1 s + D 1 Q 2 + 2 QR + 2 QV 2 QZ 2 Q f 1 s f 1 s x 2 Q + 2 R + 2 V 2 Z f 1 s x + R 2 + 2 RV 2 RZ + V 2 2 VZ + Z 2 + f 1 s 2 D 1 R 2 D 1 V 2 D 1 Z f 1 s x + f 1 s 2 + f 1 s D 1 x ,
λ 9 = 1 2 D 1 D 1 f 1 s x + Q + R + V Z + 2 f 1 s D 1 Q 2 + 2 QR + 2 QV 2 QZ 2 Q f 1 s f 1 s x 2 Q + 2 R + 2 V 2 Z f 1 s x + R 2 + 2 RV 2 RZ + V 2 2 VZ + Z 2 + f 1 s 2 D 1 R 2 D 1 V 2 D 1 Z f 1 s x + f 1 s 2 + f 1 s D 1 x .

These results lead to the following theorem.

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

2.4. Global stability and uniform persistence analysis

Now we will analyse the global stability of the equilibrium point E 2 ζ s 1 ζ s D 1 0 and the existence of positive interior equilibrium point E 3 s x y . For global stability analysis, we use the similar technique as in [3, 18].

(i) Global stability analysis of E 2 ζ s 1 ζ s D 1 0

Consider system (3) restricted to R sx + as in (5). Let N be a neighbourhood of equilibrium point E 2 ζ s 1 ζ s D 1 0 in R sx + . To analyse the global stability of the equilibrium point E 2 , a suitable Lyapunov function L = 1 2 n 1 s s ̂ 2 + n 2 x x ̂ 2 is used, where s ̂ and x ̂ denote the components of E 2 s ̂ x ̂ , i.e., s ̂ = ζ s and x ̂ = 1 ζ s D 1 , while n 1 and n 2 are positive constants. Note that L is a positive definite function with respect to E 2 in R sx + and a Lyapunov function for system (5) in N . By differentiating L with respect to time t , we get

dL dt = n 1 s s ̂ s + n 2 x x ̂ x , E7

where x ̂ = 1 s ̂ D 1 . From (5), we have

1 = s ̂ + x ̂ f 1 s ̂ and D 1 = f 1 s ̂ .

Hence (7) can be written as

L = n 1 s s ̂ s ̂ + x ̂ f 1 s ̂ s f 1 s x + n 2 x x ̂ f 1 s f 1 s ̂ 1 x = n 1 s s ̂ 2 + n 1 s s ̂ f 1 s ̂ x ̂ f 1 s x + n 2 x x ̂ x f 1 s f 1 s ̂ 1 = n 11 s s ̂ 2 + 1 2 n 12 s s ̂ x x ̂ + 1 2 n 21 s s ̂ x x ̂ + n 22 x x ̂ 2 ,

where

n 11 = n 1 < 0 ,
n 12 = n 21 = n 1 f 1 s ̂ x ̂ f 1 s x x x ̂ ,
n 22 = n 2 x f 1 s f 1 s ̂ x x ̂ f 1 s ̂ .

Clearly that L can be written as L = X T NX , which T denotes the transpose and the matrix N is particularly a real, symmetric 2 × 2 matrix, where X and N can be represented by

X = v 1 v 2 = s s ̂ x x ̂   and   N = n 11 1 2 n 12 1 2 n 21 n 22 .

Thus, it leads to the following theorem.

The equilibrium point E 2 is global asymptotically stable with respect to solution trajectories are initiated from int R sx + if the assumptions n 22 < 0 and det N > 0 are satisfied.

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

n 22 < 0 ,   and   det N = det n 11 1 2 n 12 1 2 n 21 n 22 > 0 .

This completes the proof of the theorem.

(ii) Existence of positive interior equilibrium point E 3

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 E 3 . The following lemma from Ref. [19] is applied to obtain the persistence results.

[19] Let G be an isolated hyperbolic equilibrium point in the omega limit set, ω X of the orbit O X . Then either ω X = G or there exist points J + and J in ω X with J + W s G and J W u G , where W s G and W u G denote stable and unstable manifolds of G .

Assume that

  1. E 2 ζ s 1 ζ s D 1 0 is a hyperbolic saddle point and is repels locally in y direction (as in Theorem 5),

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

  3. equilibrium point E 2 ζ s 1 ζ s D 1 0 is globally asymptotically stable with respect to R sx + .

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 + 3 is uniformly bounded. Despites that, the only compact invariant set on R + 3 is E 2 ζ s 1 ζ s D 1 0 . Let P = E 3 s x y belongs to the interior of R + 3 , i.e., P int R + 3 . We shall show that there is no points J i R + 3 where R + 3 belongs to ω P , the omega limit set of P . Now, we prove that the equilibrium point E 2 ζ s 1 ζ s D 1 0 ω P . Assume that E 2 ω P is true. Then, there is a point J 1 + W s E 2 \ E 2 such that J 1 + ω P by Lemma 2. But W s E 2 R + 3 \ E 2 is empty which is a contradiction for the positive invariance property of Η R + 3 . Thus, the equilibrium point E 2 is not in the omega limit set of P ; E 2 ω P . Next, we shall show R + 3 ω P = . Assume that R + 3 ω P , and let J R + 3 and J ω P . Then the closure of the orbit of the point J , i.e., cl O J must either contains E 2 or unbounded. This is a contradiction, and hence it is proved that R + 3 ω P = . We deduce that if E 2 ζ s 1 ζ s D 1 0 is unstable, then, for stable manifold W s E 2 ; W s E 2 int R + 3 = , and for unstable manifold W u E 2 ; W u E 2 int R + 3 . Therefore, the result of uniform persistence follows since the omega limit set of P , ω P must be in int R + 3 . This completes the proof.

Global stability of equilibrium point E 2 ζ s 1 ζ s D 1 0 with respect to R sx + 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 E 3 s x y exists (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 f 2 is a function of x and σ . Therefore, system (3) takes of the form

s = 1 s f 1 s x , x = f 1 s D 1 1 x f 2 x σ y 1 + D 1 x 1 + x + D 1 by , y = f 2 x σ y 1 + D 1 x 1 + x + D 1 by y , E8

where s 0 > 0 , x 0 > 0 , y 0 > 0 . Next, we do linearization on the system (8). First, let

S = s h 1 ; X = x h 2 ; Y = y h 3 ; s = S + h 1 ; x = X + h 2 ; y = Y + h 3 ;

where h 1 h 2 h 3 is the non-trivial equilibrium point. Then, we obtain the following differential equations

S = 1 S h 1 f 1 S + h 1 X + h 2 , X = f 1 S + h 1 D 1 1 X + h 2 f 2 X + h 2 σ Y + h 3 1 + D 1 X + h 2 1 + X + h 2 + D 1 b Y + h 3 , Y = f 2 X + h 2 σ Y + h 3 1 + D 1 X + h 2 1 + X + h 2 + D 1 b Y + h 3 Y + h 3 . E9

The Jacobian matrix of system (8) is given by

J σ = f 1 s x 1 f 1 s f 1 s x D 1 f 1 s D 1 f 2 x σ y D 1 1 + x + D 1 by D 1 x 1 + x + D 1 by 2 f 2 x σ y D 1 x 1 + x + D 1 by + 1 1 0 f 2 x σ y D 1 1 + x + D 1 by D 1 x 1 + x + D 1 by 2 + f 2 x σ y 1 + D 1 x 1 + x + D 1 by 0 f 2 x σ D 1 2 bxy 1 + x + D 1 by 2 D 1 x 1 + x + D 1 by + 1 f 2 x σ 1 + D 1 x 1 + x + D 1 by D 1 2 bxy 1 + x + D 1 by 2 1 .

Thus, the Jacobian matrix of system (8) about E 2 ζ s 1 ζ s D 1 0 σ is

J σ E 2 = f 1 ζ s 1 ζ s D 1 1 f 1 ζ s 0 f 1 ζ s 1 ζ s D 1 2 f 1 ζ s D 1 1 f 2 1 ζ s D 1 σ 1 + 1 ζ s 1 + 1 ζ s D 1 0 0 f 2 1 ζ s D 1 σ 1 + 1 ζ s 1 + 1 ζ s D 1 1 . E10

The characteristic equation of (10) is given as

λ 3 + c 1 λ 2 + c 2 λ + c 3 = 0 , E11

where

c 1 = 1 D 1 D 1 ζ s + 1 f 1 ζ s 1 + D 1 2 ζ s + ζ s 2 D 1 ζ s + f 1 ζ s ζ s D 1 1 + D 1 f 2 1 ζ s D 1 σ ζ s + D 1 ζ s 2 D 1 1 + 3 D 1 1 + D 1 ζ s ,
c 2 = 1 D 1 D 1 ζ s + 1 f 1 ζ s 2 ζ s 2 + 2 + 2 D 1 4 ζ s 2 D 1 ζ s f 2 1 ζ s D 1 σ 1 + 2 D 1 2 ζ s + ζ s 2 3 D 1 ζ s + D 1 ζ s 2 f 1 ζ s 2 + 2 D 1 2 ζ s + f 2 1 ζ s D 1 σ ζ s 2 D 1 1 + D 1 ζ s f 2 1 ζ s D 1 σ 2 D 1 + 4 D 1 2 2 D 1 ζ s 2 D 1 2 ζ s + 3 D 1 ζ s 3 D 1 2 3 D 1 ,
c 3 = 1 D 1 D 1 ζ s + 1 f 1 ζ s 1 2 ζ s + D 1 + ζ s 2 D 1 ζ s f 2 1 ζ s D 1 σ 1 2 ζ s + ζ s 2 D 1 + D 1 ζ s 2 f 1 ζ s f 2 1 ζ s D 1 σ ζ s 1 2 D 1 + D 1 ζ s + 1 ζ s f 2 1 ζ s D 1 σ D 1 D 1 ζ s + D 1 ζ s D 1 D 1 2 .

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

M 1 = c 1 ; M 2 = c 1 1 c 3 c 2 ; M 3 = c 1 1 0 c 3 c 2 c 1 0 0 c 3 .

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

c 1 > 0 c 3 > 0 c 1 c 2 c 3 > 0 . E12

When the assumptions on functional responses f 1 s = f 1 ζ s and f 2 x = f 2 1 ζ s D 1 σ as in (2) for the system (3), together with the hypotheses H 1 until H 3 ;

H 1 : D 1 ζ s + 1 > 0 ,
H 2 : f 1 ζ s 1 + D 1 2 ζ s + ζ s 2 D 1 ζ s + f 1 ζ s ζ s D 1 1 + D 1 f 2 1 ζ s D 1 σ ζ s + D 1 ζ s 2 D 1 1 + 3 D 1 1 + D 1 ζ s > 0 ,
H 3 : f 1 ζ s 1 2 ζ s + D 1 + ζ s 2 D 1 ζ s f 2 1 ζ s D 1 σ 1 2 ζ s + ζ s 2 D 1 + D 1 ζ s 2 f 1 ζ s f 2 1 ζ s D 1 σ ζ s 1 2 D 1 + D 1 ζ s + 1 ζ s f 2 1 ζ s D 1 σ D 1 D 1 ζ s + D 1 ζ s D 1 D 1 2 > 0 ,

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

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

λ 2 + c 2 λ + c 1 = 0 , E13

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

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

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

Re d λ j σ = σ 1 0 , j = 1 , 2 , 3 . E14

By substituting λ 1 σ = α σ + β σ i and λ 2 σ = α σ β σ i into 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 σ + 2 c 1 σ α σ + c 2 σ 3 β 2 σ ;
M σ = α 2 σ c 1 σ + c 2 σ α σ c 1 σ β 2 σ + c 3 σ ;
L σ = 6 α σ β σ + 2 c 1 σ β σ ;
N σ = 2 α σ β σ c 1 σ + c 2 σ β σ .

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

Re d λ j σ = σ 1 = K σ 1 M σ 1 + L σ 1 N σ 1 K 2 σ 1 + L 2 σ 1 0 , j = 1 , 2 , 3 ,

and λ 3 σ 1 = c 1 σ 1 0 . We conclude the details above in the following theorem.

Suppose that the equilibrium point E 2 ζ s 1 ζ s D 1 0 σ exists and those assumptions similar as in (2) for the system (3) together with hypothesis H 1 until H 3 hold. Then the system (8) undergoes Hopf bifurcation in the first octant, which leads to a family of periodic solutions bifurcating from E 2 ζ s 1 ζ s D 1 0 σ for some values of σ in the neighbourhood of σ 1 .

Next, we determine the Hopf bifurcation at the equilibrium point E 3 s x y σ . The Jacobian matrix of the system (8) about the equilibrium point E 3 is given by

J σ E 3 = f 1 s x 1 f 1 s f 1 s x D 1 f 1 s D 1 f 2 x σ y D 1 1 + x + D 1 b y D 1 x 1 + x + D 1 b y 2 f 2 x σ y D 1 x 1 + x + D 1 b y + 1 1 0 f 2 x σ y D 1 1 + x + D 1 b y D 1 x 1 + x + D 1 b y 2 + f 2 x σ y D 1 x 1 + x + D 1 b y + 1 0 f 2 x σ D 1 2 b x y 1 + x + D 1 b y 2 D 1 x 1 + x + D 1 b y 1 f 2 x σ D 1 x 1 + x + D 1 b y D 1 2 b x y 1 + x + D 1 b y 2 + 1 1 .

Hence, the characteristic equation for the Jacobian matrix J σ E 3 is

λ 3 + c 4 λ 2 + c 5 λ + c 6 = 0 , E16

where

c 4 = f 1 s x + Q + R + V Z + 3 f 1 s D 1 ,
c 5 = f 1 s x Q + R + V Z + 2 + f 1 s Z R 2 D 1 + 2 Q + 2 R + 2 V 2 Z + 3 ,
c 6 = f 1 s x Q + R + V Z + 1 + f 1 s Z R 1 D 1 + Q + R + V Z + 1 ,

and

Q = f 2 x σ y D 1 1 + x + D 1 b y D 1 x 1 + x + D 1 b y 2 ,
R = f 2 x σ D 1 2 b x y 1 + x + D 1 b y 2 ,
V = f 2 x σ y D 1 x 1 + x + D 1 b y + 1 ,
Z = f 2 x σ D 1 x 1 + x + D 1 b y + 1 .

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

M 4 = c 4 ; M 5 = c 4 1 c 6 c 5 ; M 6 = c 4 1 0 c 6 c 5 c 4 0 0 c 6 .

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

c 4 > 0 c 6 > 0 c 4 c 5 c 6 > 0 . E17

Suppose the assumptions of functional response f 1 s and f 2 x σ similar as in (2) for the system (3), together with the hypotheses H 4 until H 6 ;

H 4 : f 1 s x + Q + R + V + 3 > f 1 s D 1 + Z ,
H 5 : Q + R + V + 1 > Z ,
H 6 : Z > R + 1 ,

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

λ 2 + c 5 λ + c 4 = 0 , E18

Consisting of three roots; λ 1 = c 5 i , λ 2 = c 5 i , λ 3 = c 4 . As for σ σ 2 ε σ 2 + ε , all roots are in general of the form;

λ 1 σ = α σ + β σ i ,
λ 2 σ = α σ β σ i ,
λ 3 σ = c 4 σ ,

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

Re d λ j σ = σ 2 0 , j = 1 , 2 , 3 . E19

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

A 1 σ α σ A 2 σ β σ + B 1 σ = 0 , A 2 σ α σ + A 1 σ β σ + B 2 σ = 0 , E20

where

A 1 σ = 3 α 2 σ + 2 c 1 σ α σ + c 2 σ 3 β 2 σ ;
A 2 σ = 6 α σ β σ + 2 c 1 σ β σ ;
B 1 σ = α 2 σ c 1 σ + c 2 σ α σ c 1 σ β 2 σ + c 3 σ ;
B 2 σ = 2 α σ β σ c 1 σ + c 2 σ β σ .

Since A 1 σ 2 B 1 σ 2 + A 2 σ 2 B 2 σ 2 0 , we have

Re d λ j σ = σ 2 = A 1 σ 2 B 1 σ 2 + A 2 σ 2 B 2 σ 2 A 1 σ 2 2 + A 2 σ 2 2 0 , j = 1 , 2 , 3 ,

and λ 3 σ 2 = c 4 σ 2 0 . We summarise the discussion above in the following theorem.

Suppose that the equilibrium point

E 3 s x y σ = E 3 s 1 s f 1 s 2 x + 1 f 2 x x 1 b D 1 1 f 2 x σ

exists and those assumptions similar as in (2) for the system (3) together with hypotheses H 4 until H 6 hold. Then the system (8) undergoes Hopf bifurcation in the first octant, which leads to a family of periodic solutions bifurcating from E 3 s x y σ 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 nutrient s , prey organisms x and predator organisms y . We conclude that intratrophic predation denoted as b does not affected the local stability and existence of E 1 and E 2 when no predator organisms involved. Next, we have shown that the washout equilibrium point E 1 1 0 0 (no prey and predator organisms present) is locally asymptotically stable if λ 3 = f 1 1 D 1 1 < 0 holds. For stability of E 2 and E 3 of the system (3), some sufficient criteria or conditions are derived and satisfied. The points E 2 and E 3 are said to be asymptotically stable if all of their eigenvalues are less than zero.

In particular, we investigated the global stability analysis for E 2 . A suitable Lyapunov function L is defined and E 2 is globally asymptotically stable if and only if the conditions in Theorem 8 holds. Global stability of E 2 indicates that predator organisms y might 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 E 3 , 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 E 3 exists.

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 E 2 ζ s 1 ζ s D 1 0 σ and E 3 s x y σ for some values of σ in the neighbourhoods of σ 1 and σ 2 , respectively. The method used to obtain these results is similar to the method used in [3, 10].

References

  1. 1. Smith HL, Waltman P. The Theory of the Chemostat: Dynamics of Microbial Competition. New York: Cambridge University Press; 1995. pp. 1-7
  2. 2. Li B, Kuang Y. Simple food chain in a chemostat with distinct removal rates. Journal of Mathematical Analysis and Applications. 2000;242:75-92. DOI: 10.1006 jmaa.1999.6655
  3. 3. El-Sheikh MMA, Mahrouf SAA. Stability and bifurcation of a simple food chain in a chemostat with removal rates. Chaos, Solitons & Fractals. 2005;23:1475-1489. DOI: 10.1016/j.chaos.2004.06.079
  4. 4. Hamra MA, Yadi K. Asymptotic behavior of a chemostat model with constant recycle sludge concentration. Acta Biotheoretica. 2017;65:233-252. DOI: 10.1007/s10441-017-9309-4
  5. 5. Yang J, Tang G, Tang S. Modelling the regulatory system of a chemostat model with a threshold window. Mathematics and Computers in Simulation. 2017;132:220-235. DOI: 10.1016/j.matcom.2016.08.005
  6. 6. Pitchford J, Brindley J. Intratrophic predation in simple predator-prey models. Bulletin of Mathematical Biology. 1998;60:937-953; bu980053
  7. 7. Mada Sanjaya WS, Salleh Z, Mamat M. Mathematical model of three species food chain with Holling Type-III functional response. International Journal of Pure and Applied Mathematics. 2013;89(5):647-657. DOI: 10.12737/ijpam.v89i5.1
  8. 8. Tee LS, Salleh Z. Hopf bifurcation of a nonlinear system derived from Lorenz system using normal form theory. International Journal of Applied Mathematics and Statistics. 2016;55(3):122-132
  9. 9. Hassard BD, Kazarinoff ND, Wan YH. Theory and Applications of Hopf Bifurcation. USA: Cambridge University Press; 1981
  10. 10. Freedman HI, Ruan S. Hopf bifurcation in three-species food chain models with group defense. Mathematical Biosciences. 1992;111:73-87
  11. 11. Rana SMS, Nasrin F, Hossain MI. Dynamics of a predator-prey interaction in chemostat with variable yield. Journal of Sustainability Science and Management. 2015;10(2):16-23
  12. 12. Solomon ME. The natural control of animal populations. The Journal of Animal Ecology. 1949;18:1-32
  13. 13. Hale JK, Kocak H. Dynamics and Bifurcations. New York: Springer; 1991
  14. 14. Robinson CR. An Introduction to Dynamical Systems: Continuous and Discrete. Pearson Prentice Hall, New Jersey: Pearson Education, Inc; 2004
  15. 15. Hirsch MW, Smale S, Devaney RL. Differential Equations, Dynamical Systems & an Introduction to Chaos. 2nd ed. Elsevier Academic Press, San Diego California: Elsevier; 2004
  16. 16. Freedman HI, Mathsen RM. Persistence in predator-prey systems with ratio-dependent predator influence. Bulletin of Mathematical Biology. 1993;55(4):817-827
  17. 17. Butler GJ, Hsu SB, Waltman P. Coexistence of competing predators in a chemostat. Journal of Mathematical Biology. 1983;17:133-151
  18. 18. Nani F, Freedman HI. A mathematical model of cancer treatment by immunotherapy. Mathematical Biosciences. 2000;163:159-199
  19. 19. Freedman HI, Waltman P. Persistence in models of three interacting predator-prey populations. Mathematical Biosciences. 1984;68:213-231
  20. 20. Freedman HI, Ruan S. Uniform persistence in functional differential equations. Mathematical Biosciences. 1995;115:173-192
  21. 21. Gantmacher FR. The Theory of Matrices. USA: Chelsea Publishing Company; 1964
  22. 22. Tee LS, Salleh Z. Hopf bifurcation analysis of Zhou system. International Journal of Pure and Applied Mathematics. 2015;104(1):1-18
  23. 23. Marsden JE, McKracken M. The Hopf Bifurcation and its Applications. New York: Springer-Verlag; 1976

Written By

Zabidin Salleh and Liyana Abd Rahim

Submitted: 08 April 2017 Reviewed: 11 October 2017 Published: 20 December 2017