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

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.


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 þ D1x 1þxþD1by are added to differential equations x 0 and y 0 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.

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 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, 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., ð Þ 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.
To make our works more convenient, we just drop the bars and tildes. So after rescaling, the system (1) becomes where s 0 ð Þ > 0, x 0 ð Þ > 0, y 0 ð Þ > 0 and D 1 > 0: Now, the variables are non-dimensional and the discussion is in

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).
Definition 1 [13]. A system with differential equations 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 : ii.
All nonnegative solutions of (3) with initial values in R 3 þ are uniformly bounded and they eventually attracted into region Η: iii. The system is dissipative.
Remark 1. 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 Pmax À q ≤ s þ x þ y ≤ 1 Pmin þ q: Let us define By adding those differential equations s 0 , x 0 and y 0 in (3), we shall get Hence, this leads to Thus, solving the above inequality gives as q is the positive constant (ii) and (iii): Let 0 < s 0 < 1 and consider.

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 0 , x 0 and y 0 to zero and solve the resulting equations simultaneously. The possible equilibrium points are as follows; i. 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 Thus, E 2 is the one possible equilibrium point that consists only prey organisms and not predator organisms.
iii. 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 We solve for y, and get y * ¼ To show the existence of E 2 , we let the system (3) be restricted to R þ sx as where s 0 ð Þ > 0 and x 0 ð Þ > 0. Thus, Lemma 1 below shows the existence of non-trivial equilibrium point E 2 : Proof. We will get two surfaces by equating system (5) to zero; 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.

Definition 2. [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 E > 0, there is δ > 0, such that if 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.
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 .
. 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., DF 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.

Theorem 2. [14]
Let _ x ¼ F x ð Þ be a differential equation in n variables, with a hyperbolic fixed point x * .
Suppose thal F, ∂Fi ∂xj x ð Þ and ∂ 2 Fi ∂xj∂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.
a. If the real parts of all the eigenvalues of DF 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.
b. If there is at least one eigenvalue of DF 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.) c. If one of the eigenvalues of DF 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 AEβ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 Then, As Jacobian matrix J E1 above is an upper triangular, therefore the diagonal is the value of all of its eigenvalues; λ 1 ¼ À1, λ 2 ¼ À1 and λ 3 ¼ If all the eigenvalues of J E1 are negatives, then this leads to the following result.
(ii) Stability analysis of E 2 : Next, we analyse the local stability of the system (3)  ; 0 is given as We can see from both Jacobian matrices J E1 and J E2 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 Then, by using MATLAB R2015b, we get the eigenvalues of J E2 which are λ 4 ¼ À1 , where D 1 À ζ s þ 1 6 ¼ 0: We summarise the above discussion in the following theorems. ; 0 is locally asymptotically stable.
; 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 Proof. The results follow from inspections of the eigenvalues of the matrix J E2 and Theorem 2 (see [16,18]). ∎

Definition 4. [17]
The flow F will be called uniformly persistent if there exists ε 0 > 0 such that The following theorem shows the existence of the equilibrium point E 2 using persistence analysis. Theorem 6. Assume that i. Lemma 1 being holds, ii. E 2 is a unique hyperbolic saddle point in R þ sxy and repels locally in y-direction (as in Theorem 5), iii. no existence of periodic orbits in the planes of R þ sxy . Then where k > 0: 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 : The eigenvalues of J E3 are resulted from the characteristic equation below where and From (6), we obtain the eigenvalues; λ 7 ¼ À1, and λ 8 , λ 9 are Complexity in Biological and Physical Systems -Bifurcations, Solitons and Fractals These results lead to the following theorem.

Global stability and uniform persistence analysis
Now we will analyse the global stability of the equilibrium point E 2 ζ s ; 1Àζ s ð Þ D1 ; 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 D1 ; 0 Consider system (3) restricted to R þ sx as in (5). Let N be a neighbourhood of equilibrium point To analyse the global stability of the equilibrium point E 2 , a suitable , 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 where b x ¼ 1Àb s D1 . From (5), we have Hence (7) can be written as where n 11 ¼ Àn 1 < 0, Clearly that L 0 can be written as L 0 ¼ 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 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.

Lemma 2. [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.
; 0 is a hyperbolic saddle point and is repels locally in y direction (as in Theorem 5), ii.
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 ; 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 Þ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 ð Þ 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 ; 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]).
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 ð Þ D1 ; 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].