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

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.

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 $\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)$ are added to differential equations $x^{\prime }$ and $y^{\prime }$ 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.

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

$$\{\begin{array}{l}{s}^{\prime }\left(t\right)=\left(s^0-s\right)D_0-\frac{1}{{\beta }_1}{f}_1\left(s\right)x,\hfill \\ x^{\prime }\left(t\right)=\left( f_1\left(s\right)-D_1\right)x-\frac{1}{{\beta }_2}{f}_2\left(x\right)\left(y+\frac{D_{1}xy}{D_2+x+by}\right),\hfill \\ y^{\prime }\left(t\right)=f_2\left(x\right)\left(y+\frac{D_{1}xy}{D_2+x+by}\right)-y,\hfill \end{array}$$E1

where $$s\left(0\right)=s_0>0,$$ $$x\left(0\right)=x_0>0,$$ $$y\left(0\right)=y_0>0,$$ $$D_i=1;$$ $$i=0,1,2,$$ $${\beta }_j\le 1;$$ $$j=1,2$$ and $$0\le b\le 1.$$ This leads to the following assumptions on functional response $$f_n,n=1,2$$ :

$$\begin{array}{l}\left(i\right)f_n:{\mathbb{R}}^{+}\to {\mathbb{R}}^{+}\text{and}{f}_n,n=1,2\mathrm{are}\text{continuously differential equations},\\ \left(\mathit{ii}\right)f_n\left(0\right)=0,\\ \left(\mathit{iii}\right)f_1^{\prime }\left(s\right)>0,f_2^{\prime }\left(x\right)>0\text{for}\mathrm{all}s,x>0.\end{array}$$ 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 $${\beta }_1$$ and $${\beta }_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 $${\theta }_i$$ , i.e., $$D_i=D_0+{\theta }_i,i=1,2.$$ $$f_1\left(s\right)$$ and $$f_2\left(x\right)$$ 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.

$$\overline{s}=\frac{s}{s^0},$$	$$\overline{x}=\frac{x}{D_0{\beta }_1{s}^0},$$	$$\overline{y}=\frac{y}{D_0{D}_1{\beta }_1{\beta }_2{s}^0},$$
$$f_1\left(s\right)=\frac{\widetilde{f_1}\left(s^0\overline{s}\right)}{D_0},$$	$$f_2\left(x\right)=\frac{\widetilde{f_2}\left(D_0{\beta }_1{s}^0\overline{x}\right)}{D_0},$$	$$b=\frac{\overline{b}}{{\beta }_2},$$
$$\overline{t}=D_{i}t,i=0,1,$$	$$D_2=D_0{\beta }_1{s}^0.$$

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

$$\left\{\begin{array}{l}{s}^{\prime }=1-s-f_1\left(s\right)x,\\ x^{\prime }=\left(\frac{f_1\left(s\right)}{D_1}-1\right)x-f_2\left(x\right)y\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right),\\ y^{\prime }=f_2\left(x\right)y\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)-y,\end{array}\right.$$ E3

where

$$s\left(0\right)>0,x\left(0\right)>0,y\left(0\right)>0\text{and}{D}_1>0.$$ E4

Now, the variables are non-dimensional and the discussion is in $${\mathbb{R}}_{+}^3=\left\{\left(s,x,y\right):s>0,x>0,y>0\right\}$$ .

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 $${\mathbb{R}}_{+}^3$$ . Firstly, we consider dissipativity of the system (3). Definition 1

[13]. A system with differential equations $$x^{\prime }=f\left(x\right)$$ is defined to be dissipative if there exists a bounded subset $$\mathrm{\Gamma }$$ of $${\mathbb{R}}^3$$ , such that there is a time $$t_0$$ for any $$x^0\in {\mathbb{R}}^3$$ which depends on $$x^0$$ and $$\mathrm{\Gamma }$$ so that the solution of the system $$\varphi \left(t,x^0\right)\in \mathrm{\Gamma }$$ for $$t\ge t_0.$$

Theorem 1.

Let $$H$$ be the region

$$H=\left\{\left(s,x,y\right)\in {\mathbb{R}}_{+}^3:\frac{1}{P_{\mathit{max}}}-q\le s+x+y\le \frac{1}{P_{\mathit{min}}}+q\right\},$$

where $$P_{\mathit{max}}=\text{max}\left\{1,f_1\left(s\right)+\frac{f_1\left(s\right)}{D_1}-1\right\},P_{\mathit{min}}=\text{min}\left\{1,f_1\left(s\right)+\frac{f_1\left(s\right)}{D_1}-1\right\}$$ and $$q$$ is a positive constant.

Then,

1. $${\rm H}$$ is positively invariant.

2. All nonnegative solutions of (3) with initial values in $${\mathbb{R}}_{+}^3$$ are uniformly bounded and they eventually attracted into region $${\rm H}.$$

3. 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 $$\frac{1}{P_{\mathit{max}}}-q\le s+x+y\le \frac{1}{P_{\mathit{min}}}+q.$$ Let us define

$$P_{\mathit{max}}=\text{max}\left\{1,f_1\left(s\right)+\frac{f_1\left(s\right)}{D_1}-1\right\} \mathrm{and} P_{\mathit{min}}=\text{min}\left\{1,f_1\left(s\right)+\frac{f_1\left(s\right)}{D_1}-1\right\}.$$

By adding those differential equations $$s^{\prime },x^{\prime }$$ and $$y^{\prime }$$ in (3), we shall get

$$s^{\prime }+x^{\prime }+y^{\prime }=1-\left[s+\left(f_1\left(s\right)+\left(\frac{f_1\left(s\right)}{D_1}-1\right)\right)x+y\right].$$

Hence, this leads to

$$1-P_{\mathit{max}}\left(s+x+y\right)\le 1-\left[s+\left(f_1\left(s\right)+\left(\frac{f_1\left(s\right)}{D_1}-1\right)\right)x+y\right]\le 1-P_{\mathit{min}}\left(s+x+y\right).$$

Thus, solving the above inequality gives

$$\frac{1}{P_{\mathit{max}}}-q\le s+x+y\le \frac{1}{P_{\mathit{min}}}+q,$$

as $$q$$ is the positive constant

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

$$s^{\prime }=1-s-f_1\left(s\right)x<1-s.$$

Then $$s\left(t\right)<1-\frac{e^{-t}}{1-s_0}$$ for $$0<s_0<1,$$ and hence $$\underset{t\to \infty }{\text{lim}}\text{sup}s\left(t\right)<1$$ for $$0<s_0<1.$$

Now, consider the $$x^{\prime }$$ equation;

$$x^{\prime }=\left(\frac{f_1\left(s\right)}{D_1}-1\right)x-f_2\left(x\right)y\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)<\left(f_1\left(s\right)-1\right)x<{\alpha }_{1}x,$$

where $${\alpha }_1=\underset{s\in {\rm H}}{\mathit{max}}\left( f_1\left(s\right)-1\right).$$ We assume that $${\alpha }_1=\left(\underset{s\in {\rm H}}{\mathit{max}}\left( f_1\left(s\right)-1\right)\right)<0.$$ Hence, $$x\left(t\right)<x_0{e}^{{\alpha }_{1}t},{\alpha }_1<0,$$ and thus, $$\underset{t\to \infty }{\mathit{lim}}\mathit{sup}x\left(t\right)<x_0,$$ for $$x_0>0.$$

Similarly, consider the $$y^{\prime }$$ equation,

$$y^{\prime }=f_2\left(x\right)y\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)-y<{\alpha }_{2}y,$$

where $${\alpha }_2=\underset{x\in {\rm H}}{\mathit{max}}\left(f_2\left(x\right)+\frac{f_2\left(x\right)D_{1}x}{1+x+D_1\mathit{by}}\right).$$ We shall assume that $$\underset{x\in {\rm H}}{\mathit{max}}\left(f_2\left(x\right)+\frac{f_2\left(x\right)D_{1}x}{1+x+D_1\mathit{by}}\right)<0.$$ Then, $$y\left(t\right)<y_0{e}^{{\alpha }_{2}t}$$ where $${\alpha }_2<0.$$ Hence, $$\underset{t\to \infty }{\mathit{lim}}\mathit{sup} y\left(t\right)<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^{\prime },x^{\prime }$$ and $$y^{\prime }$$ to zero and solve the resulting equations simultaneously. The possible equilibrium points are as follows;

1. $$E_1\left(1,0,0\right)$$ 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\left(s\right)x=0$$ . Then, from equation $$x^{\prime }=0,$$ we get $$x=0$$ when $$y=0.$$

2. $$E_2\left({\zeta }_s,\frac{1-{\zeta }_s}{D_1},0\right)$$ . $$s={\zeta }_s$$ is the unique solution of $$\frac{f_1\left(s\right)}{D_1}-1=0.$$ Let $$y=0,$$ then Eq. (3) give $$1-s-f_1\left(s\right)x=0$$ and $$\left(\frac{f_1\left(s\right)}{D_1}-1\right)x=0.$$ When x $$\ne 0,$$ $$\frac{f_1\left(s\right)}{D_1}-1=0$$ . Therefore, $$f_1\left(s\right)=D_1.$$ Next, we find $$x.$$ From equation $$1-s-f_1\left(s\right)x=0$$ , $$x=\frac{1-s}{f_1\left(s\right)}.$$ By substituting $$f_1\left(s\right)=D_1$$ and $$s={\zeta }_s$$ , we get $$x=\frac{1-{\zeta }_s}{D_1}.$$ Thus, $$E_2$$ is the one possible equilibrium point that consists only prey organisms and not predator organisms.

3. $$E_3\left(s^{\ast },x^{\ast },y^{\ast }\right)$$ denotes the positive interior equilibrium point with $$s^{\ast },x^{\ast },y^{\ast }>0$$ . $$s^{\ast }$$ is a unique solution of $$1-s-f_1\left(s\right)x=0.$$ From this, we have $$x^{\ast }=\frac{1-s^{\ast }}{f_1\left(s^{\ast }\right)}$$ . Next, let the equation $$y^{\prime }=0$$ . Since $$y\ne 0,$$ $$f_2\left(x\right)\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)=1$$ . We solve for $$y,$$ and get $$y^{\ast }=\frac{\left(2x^{\ast }+1\right)f_2\left(x^{\ast }\right)-x^{\ast }-1}{b{D}_1\left(1-f_2\left(x^{\ast }\right)\right)}$$ . Thus, the equilibrium point $$E_3\left(s^{\ast },x^{\ast },y^{\ast }\right)$$ is $$\left(s^{\ast },\frac{1-s^{\ast }}{f_1\left(s^{\ast }\right)},\frac{\left(2x^{\ast }+1\right)f_2\left(x^{\ast }\right)-x^{\ast }-1}{b{D}_1\left(1-f_2\left(x^{\ast }\right)\right)}\right)$$ , where $$x^{\ast }=\frac{1-s^{\ast }}{f_1\left(s^{\ast }\right)}$$ .

   To show the existence of $$E_2,$$ we let the system (3) be restricted to $${\mathbb{R}}_{\mathit{sx}}^{+}$$ as

$$\left\{\begin{array}{l}{s}^{\prime }=1-s-f_1\left(s\right)x,\\ x^{\prime }=\left(\frac{f_1\left(s\right)}{D_1}-1\right)x,\end{array}\right.$$ E5

where $$s\left(0\right)>0$$ and $$x\left(0\right)>0$$ . Thus, Lemma 1 below shows the existence of non-trivial equilibrium point $$E_2.$$ Lemma 1.

Suppose that a point $$\left({\zeta }_s,\widehat{x}\right)$$ exists in $${\mathbb{R}}_{\mathit{sx}}^{+}$$ such that $${\zeta }_s+D_1\widehat{x}-1=0$$ as time, $$t$$ approaches $$\infty .$$ Then, the non-trivial equilibrium point $$E_2\left({\zeta }_s,\frac{1-{\zeta }_s}{D_1},0\right)$$ exists.

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

$$1-s-f_1\left(s\right)x=0,$$

$$\left(\frac{f_1\left(s\right)}{D_1}-1\right)x=0.$$

Then,

$$\frac{1-s}{x}=f_1\left(s\right)\text{and}{D}_1=f_1\left(s\right).$$

Thus, $$1-s=D_{1}x$$ , i.e., $$x=\frac{1-s}{D_1}$$ . By taking $$s={\zeta }_s$$ , we will have $$x=\widehat{x}=\frac{1-{\zeta }_s}{D_1}$$ and satisfying $${\zeta }_s+D_1\widehat{x}-1=0$$ . Hence, the equilibrium point $$E_2\left({\zeta }_s,\frac{1-{\zeta }_s}{D_1},0\right)$$ exists. $$∎$$

2.3. Stability analysis of equilibrium points $${\mathit{E}}_1,{\mathit{E}}_2$$ and $${\mathit{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^{\ast }$$ is Lyapunov stable or L-stable, provided that any solution $$\varphi \left(t;x_0\right)$$ stays near $$x^{\ast }$$ for all future time $$t\ge 0$$ if the initial condition $$x_0$$ starts near enough to $$x^{\ast }.$$ Specifically, a fixed point $$x^{\ast }$$ is L-stable, provided that for any $$ϵ>0,$$ there is $$\delta >0$$ , such that if $$‖x_0-x^{\ast }‖<\delta ,$$ then $$‖\varphi \left(t;x_0\right)-x^{\ast }‖<ϵ$$ for all $$t\ge 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 $$\omega$$ -limit set as follows. Definition 3.

[14] A point $$k$$ is an $$\omega$$ -limit point of the trajectory of $$x_0$$ , if $$\varphi \left(t;x_0\right)$$ keeps coming near $$k$$ as $$t\to \infty$$ (i.e., there is a sequence of times $$t_j$$ , with $$t_j\to \infty$$ as $$j\to \infty$$ such that $$\varphi \left(t_j;x_0\right)$$ converges to $$k$$ ). Certainly, if $$‖\varphi \left(t;x_0\right)-x^{\ast }‖\to 0$$ as $$t\to \infty$$ , then $$x^{\ast }$$ is the only $$\omega$$ -limit point of $$x_0$$ . There can be more than one point that is an $$\omega$$ -limit point of $$x_0$$ . The set of all $$\omega$$ -limit points of $$x_0$$ is denoted by $$\omega \left(x_0\right)$$ and is called the omega limit set of $$x_0$$ .

Definition 4.

[14] A fixed point $$x^{\ast }$$ is weakly asymptotically stable, if there exists $${\delta }_1>0$$ such that $$\omega \left(x_0\right)=\left\{x^{\ast }\right\}$$ for all $$‖x_0-x^{\ast }‖<{\delta }_1$$ (i.e., $$‖\varphi \left(t;x_0\right)-x^{\ast }‖\to 0$$ as time $$t\to \infty$$ for all $$‖x_0-x^{\ast }‖<{\delta }_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^{\ast }$$ 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^{\ast }$$ is hyperbolic (the real parts of all eigenvalues of the Jacobian matrix at $$x^{\ast }$$ , i.e., $$D{F}_{\left(x^{\ast }\right)}$$ 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 $$\dot{x}=F\left(x\right)$$ be a differential equation in $$n$$ variables, with a hyperbolic fixed point $$x^{\ast }$$ . Suppose thal $$F,\frac{\partial F_i}{\partial x_j}\left(x\right)$$ and $$\frac{{\partial }^2{F}_i}{\partial x_j\partial x_k}\left(x\right)$$ 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}_{\left(x^{\ast }\right)}$$ 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^{\ast }$$ is asymptotically stable for the nonlinear equation). In this case, the basin of attraction $$W^S\left(x^{\ast }\right)$$ is an open set that contains some solid ball about the fixed point.

2. If there is at least one eigenvalue of $$D{F}_{\left(x^{\ast }\right)}$$ has positive real part, then the fixed point $$x^{\ast }$$ 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}_{\left(x^{\ast }\right)}$$ has zero real part, then the situation is more complicated. In particular, for $$n=2,$$ if the fixed point is an elliptic center (eigenvalues $$\pm \mathit{\beta 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\left(1,0,0\right)$$ . The Jacobian matrix of the system (3) is

$$J=\left[\begin{array}{cc}-f_1^{\prime }\left(s\right)x-1& -f_1\left(s\right)\\ \frac{f_1^{\prime }\left(s\right)x}{D_1}& \frac{f_1\left(s\right)}{D_1}-f_2\left(x\right)y\left[\left(\frac{D_1}{1+x+D_1\mathit{by}}\right)-\left(\frac{D_{1}x}{{\left(1+x+D_1\mathit{by}\right)}^2}\right)\right]-{f_2}^{\prime }\left(x\right)y\left(\left(\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)+1\right)-1\\ 0& f_2\left(x\right)y\left[\left(\frac{D_1}{1+x+D_1\mathit{by}}\right)-\left(\frac{D_{1}x}{{\left(1+x+D_1\mathit{by}\right)}^2}\right)\right]+f_2^{\prime }\left(x\right)y\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)\end{array}\begin{array}{c}0\\ \left(\frac{f_2\left(x\right)D_1^2\mathit{bxy}}{{\left(1+x+D_1\mathit{by}\right)}^2}\right)-f_2\left(x\right)\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)\\ f_2\left(x\right)\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)-\left(\frac{f_2\left(x\right)D_1^2\mathit{bxy}}{{\left(1+x+D_1\mathit{by}\right)}^2}\right)-1\end{array}\right].$$

Then,

$$J_{E_1}=\left(\begin{array}{ccc}-1& -f_1\left(1\right)& \hfill 0\\ \hfill 0& \frac{f_1\left(1\right)}{D_1}-1& \hfill 0\\ \hfill 0& 0& -1\end{array}\right).$$

As Jacobian matrix $$J_{E_1}$$ above is an upper triangular, therefore the diagonal is the value of all of its eigenvalues; $${\lambda }_1=-1,{\lambda }_2=-1$$ and $${\lambda }_3=\frac{f_1\left(1\right)}{D_1}-1.$$ If all the eigenvalues of $${J_E}_1$$ are negatives $$,$$ then this leads to the following result. Theorem 3.

If the eigenvalue $${\lambda }_3$$ such that

$${\lambda }_3=\frac{f_1\left(1\right)}{D_1}-1<0,$$

then the trivial equilibrium point $$E_1\left(1,0,0\right)$$ 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\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ . The Jacobian matrix at $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ is given as

$$J_{E_2}=\left(\begin{array}{ccc}-\frac{f_1^{\prime }\left({\zeta }_s\right)\left(1-{\zeta }_s\right)}{D_1}-1& -f_1\left({\zeta }_s\right)& 0\\ \frac{f_1^{\prime }\left({\zeta }_s\right)\left(1-{\zeta }_s\right)}{D_1^2}& \frac{f_1\left({\zeta }_s\right)}{D_1}-1& -f_2\left(\frac{1-{\zeta }_s}{D_1}\right)\left(1+\left(\frac{1-{\zeta }_s}{1+\left(\frac{1-{\zeta }_S}{D_1}\right)}\right)\right)\\ 0& 0& f_2\left(\frac{1-{\zeta }_s}{D_1}\right)\left(1+\left(\frac{1-{\zeta }_s}{1+\left(\frac{1-{\zeta }_S}{D_1}\right)}\right)\right)-1\end{array}\right).$$

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 $${\lambda }^3+c_1{\lambda }^2+c_2\lambda +c_3=0$$ where

$$c_1=\frac{1}{D_1\left(D_1-{\zeta }_s+1\right)}\left(f_1^{\prime }\left({\zeta }_s\right)+3D_1-f_1\left({\zeta }_s\right)+f_1^{\prime }\left({\zeta }_s\right)D_1-D_1{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)-2{\zeta }_s{f}_1^{\prime }\left({\zeta }_s\right)-3D_1{\zeta }_s-2D_1^2{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)+{\zeta }_s^2{f}_1^{\prime }\left({\zeta }_s\right)-f_1\left({\zeta }_s\right)D_1+{\zeta }_s{f}_1\left({\zeta }_s\right)+3D_1^2-D_1{\zeta }_s{f}_1^{\prime }\left({\zeta }_s\right)+D_1{\zeta }_s{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)+D_1^2{\zeta }_s{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)\right.),$$

$$c_2=\frac{1}{D_1\left(D_1-{\zeta }_S+1\right)}\left(2f_1\left({\zeta }_s\right)-3D_1-2f_1^{\prime }\left({\zeta }_s\right)-2D_1{f}_1^{\prime }\left({\zeta }_s\right)+f_2\left(\frac{1-{\zeta }_s}{D_1}\right)f_1^{\prime }\left({\zeta }_s\right)+2D_1{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)+4{\zeta }_s{f}_1^{\prime }\left({\zeta }_s\right)+3D_1{\zeta }_s+4D_1^2{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)-2{\zeta }_s^2{f}_1^{\prime }\left({\zeta }_s\right)+2D_1{f}_1\left({\zeta }_s\right)-f_1\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)-2{\zeta }_s{f}_1\left({\zeta }_s\right)-3D_1^2-2D_1{f}_1\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+{\zeta }_s{f}_1\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+2D_1{f}_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+2D_1{\zeta }_s{f}_1^{\prime }\left({\zeta }_s\right)-2{\zeta }_s{f}_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)-2D_1{\zeta }_s{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)+{\zeta }_s^2{f}_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)-2D_1^2{\zeta }_s{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)-3D_1{\zeta }_s{f}_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+D_1{\zeta }_s^2{f}_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+D_1{\zeta }_s{f}_1\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)\right),$$

$$c_3=\frac{1}{D_1\left(D_1-{\zeta }_s+1\right)}\left(f_1\left({\zeta }_s\right)-D_1-f_1^{\prime }\left({\zeta }_s\right)-D_1{f}_1^{\prime }\left({\zeta }_s\right)+f_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+D_1{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)+2{\zeta }_s{f}_1^{\prime }\left({\zeta }_s\right)+D_1{\zeta }_s-{\zeta }_s^2{f}_1^{\prime }\left({\zeta }_s\right)+f_1\left({\zeta }_s\right)D_1-f_1\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)-{\zeta }_s{f}_1\left({\zeta }_s\right)-D_1^2-2D_1{f}_1\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+{\zeta }_s{f}_1\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+2D_1{f}_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+D_1{\zeta }_s{f}_1^{\prime }\left({\zeta }_s\right)-2{\zeta }_s{f}_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)-D_1{\zeta }_s{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)+{\zeta }_s^2{f}_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)-D_1^2{\zeta }_s{f}_2\left(\frac{1-{\zeta }_s}{D_1}\right)-3D_1{\zeta }_s{f}_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+D_1{\zeta }_s^2{f}_1^{\prime }\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)+D_1{\zeta }_s{f}_1\left({\zeta }_s\right)f_2\left(\frac{1-{\zeta }_s}{D_1}\right)\right).$$

Then, by using MATLAB R2015b, we get the eigenvalues of $$J_{E_2}$$ which are $${\lambda }_4=-1$$ , $${\lambda }_5=\frac{f_1^{\prime }\left({\zeta }_s\right)\left({\zeta }_s-1\right)+f_1\left({\zeta }_s\right)-D_1}{D_1}$$ and $${\lambda }_6=\frac{f_2\left(\frac{1-{\zeta }_s}{D_1}\right)\left(1+2D_1-{\zeta }_s-D_1{\zeta }_s\right)-1}{D_1-{\zeta }_s+1}$$ where $$D_1-{\zeta }_s+1\ne 0.$$ We summarise the above discussion in the following theorems. Theorem 4.

Suppose the assumptions such as

$$\frac{f_1^{\prime }\left({\zeta }_s\right)\left({\zeta }_s-1\right)+f_1\left({\zeta }_s\right)-D_1}{D_1}<0$$ and $$\frac{f_2\left(\frac{1-{\zeta }_s}{D_1}\right)\left(1+2D_1-{\zeta }_s-D_1{\zeta }_s\right)-1}{D_1-{\zeta }_s+1}<0$$

are satisfied, then the equilibrium point, $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ is locally asymptotically stable.

Theorem 5.

If

$$\frac{f_1^{\prime }\left({\zeta }_s\right)\left({\zeta }_s-1\right)+f_1\left({\zeta }_s\right)-D_1}{D_1}<0$$ and $$\frac{f_2\left(\frac{1-{\zeta }_s}{D_1}\right)\left(1+2D_1-{\zeta }_s-D_1{\zeta }_s\right)-1}{D_1-{\zeta }_s+1}<0$$ ,

then the equilibrium point $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ is a hyperbolic saddle and is repels locally in $$y$$ -direction. Particularly, the dimension of the stable manifold $$W^s\left(E_2\right)$$ and the unstable manifold $$W^u\left(E_2\right)$$ are given by

$$\text{dim}{W}^s\left(E_2\right)=2$$ and $$\text{dim}{W}^u\left(E_2\right)=1$$ .

Proof. The results follow from inspections of the eigenvalues of the matrix $$J_{E_2}$$ and Theorem 2 (see [16, 18]). $$∎$$ Definition 4.

[17] The flow $$\mathcal{F}$$ will be called uniformly persistent if there exists $${\epsilon }_0>0$$ such that for all $$x\in E^0$$ , $$\underset{\mathrm{t}\to \infty }{\text{lim}}d\left(\pi \left(x,t\right),\partial E\right)\ge {\epsilon }_0.$$

The following theorem shows the existence of the equilibrium point $$E_2$$ using persistence analysis. Theorem 6.

Assume that

1. Lemma 1 being holds,

2. $$E_2$$ is a unique hyperbolic saddle point in $${\mathbb{R}}_{\mathit{sxy}}^{+}$$ and repels locally in $$y$$ -direction (as in Theorem 5),

3. no existence of periodic orbits in the planes of $${\mathbb{R}}_{\mathit{sxy}}^{+}$$ .

Then

$$\underset{t\to +\infty }{\mathit{lim}}\mathit{inf}s\left(t\right)>k,\underset{t\to +\infty }{\mathit{lim}}\mathit{inf}x\left(t\right)>k,\underset{t\to +\infty }{\mathit{lim}}\mathit{inf}y\left(t\right)>k,$$

where $$k>0.$$

Particularly, the system (3) exhibits uniform persistence and the equilibrium point $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ 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\left(s^{\ast },x^{\ast },y^{\ast }\right).$$

The Jacobian matrix at $$E_3$$ is

$$J_{E_3}=\left[\begin{array}{cc}-f_1^{\prime }\left(s^{\ast }\right)x^{\ast }-1& -f_1\left(s^{\ast }\right)\\ \frac{f_1^{\prime }\left(s^{\ast }\right)x^{\ast }}{D_1}& \frac{f_1\left(s^{\ast }\right)}{D_1}-f_2\left(x^{\ast }\right)y^{\ast }\left[\left(\frac{D_1}{1+x^{\ast }+D_{1}b{y}^{\ast }}\right)-\left(\frac{D_1{x}^{\ast }}{{\left(1+x^{\ast }+D_{1}b{y}^{\ast }\right)}^2}\right)\right]-{f_2}^{\prime }\left(x^{\ast }\right)y^{\ast }\left(\frac{D_1{x}^{\ast }}{1+x^{\ast }+D_{1}b{y}^{\ast }}+1\right)-1\\ 0& f_2\left(x^{\ast }\right)y^{\ast }\left[\left(\frac{D_1}{1+x^{\ast }+D_{1}b{y}^{\ast }}\right)-\left(\frac{D_1{x}^{\ast }}{{\left(1+x^{\ast }+D_{1}b{y}^{\ast }\right)}^2}\right)\right]+f_2^{\prime }\left(x^{\ast }\right)y^{\ast }\left(\frac{D_1{x}^{\ast }}{1+x^{\ast }+D_{1}b{y}^{\ast }}+1\right)\end{array}\begin{array}{c}0\\ \left(\frac{f_2\left(x^{\ast }\right)D_1^{2}b{x}^{\ast }{y}^{\ast }}{{\left(1+x^{\ast }+D_{1}b{y}^{\ast }\right)}^2}\right)-f_2\left(x^{\ast }\right)\left(\frac{D_1{x}^{\ast }}{1+x^{\ast }+D_{1}b{y}^{\ast }}+1\right)\\ f_2\left(x^{\ast }\right)\left(\frac{D_1{x}^{\ast }}{1+x^{\ast }+D_{1}b{y}^{\ast }}+1\right)-\left(\frac{f_2\left(x^{\ast }\right)D_1^{2}b{x}^{\ast }{y}^{\ast }}{{\left(1+x^{\ast }+D_{1}b{y}^{\ast }\right)}^2}\right)-1\end{array}\right].$$

The eigenvalues of $$J_{E_3}$$ are resulted from the characteristic equation below

$${\lambda }^3+c_4{\lambda }^2+c_5\lambda +c_6=0,$$ E6

where

$$c_4=f_1^{\prime }\left(s^{\ast }\right)x^{\ast }+Q+R+V-Z+3-\frac{f_1\left(s^{\ast }\right)}{D_1},$$

$$c_5=f_1^{\prime }\left(s^{\ast }\right)x^{\ast }\left(Q+R+V-Z+2\right)+\frac{f_1\left(s^{\ast }\right)\left(Z-R-2\right)}{D_1}+2Q+2R+2V-2Z+3,$$

$$c_6=f_1^{\prime }\left(s^{\ast }\right)x^{\ast }\left(Q+R+V-Z+1\right)+\frac{f_1\left(s^{\ast }\right)\left(Z-R-1\right)}{D_1}+Q+R+V-Z+1,$$

and

$$Q=f_2\left(x^{\ast }\right)y^{\ast }\left[\left(\frac{D_1}{1+x^{\ast }+D_{1}b{y}^{\ast }}\right)-\left(\frac{D_1{x}^{\ast }}{{\left(1+x^{\ast }+D_{1}b{y}^{\ast }\right)}^2}\right)\right],$$

$$R=\left(\frac{f_2\left(x^{\ast }\right)D_1^{2}b{x}^{\ast }{y}^{\ast }}{{\left(1+x^{\ast }+D_{1}b{y}^{\ast }\right)}^2}\right),$$

$$V=f_2^{\prime }\left(x^{\ast }\right)y^{\ast }\left(\frac{D_1{x}^{\ast }}{1+x^{\ast }+D_{1}b{y}^{\ast }}+1\right),$$

$$Z=f_2\left(x^{\ast }\right)\left(\frac{D_1{x}^{\ast }}{1+x^{\ast }+D_{1}b{y}^{\ast }}+1\right).$$

From (6), we obtain the eigenvalues; $${\lambda }_7=-1,$$ and $${\lambda }_8$$ , $${\lambda }_9$$ are

$${\lambda }_8=-\frac{1}{2D_1}\left(D_1\left(f_1^{\prime }\left(s^{\ast }\right)x^{\ast }+Q+R+V-Z+2\right)-f_1\left(s^{\ast }\right)+\left(\sqrt{\begin{array}{c}{D}_1\left(\begin{array}{c}{Q}^2+2\mathit{QR}+2\mathit{QV}-2\mathit{QZ}-2Q{f}_1\left(s^{\ast }\right)-f_1^{\prime }\left(s^{\ast }\right)x^{\ast }\left(2Q+2R+2V-2Z-f_1^{\prime }\left(s^{\ast }\right)x^{\ast }\right)+R^2+\\ 2\mathit{RV}-2\mathit{RZ}+V^2-2\mathit{VZ}+Z^2\end{array}\right)+\\ f_1\left(s^{\ast }\right)\left(2D_{1}R-2D_{1}V-2D_{1}Z-f_1^{\prime }\left(s^{\ast }\right)x^{\ast }\right)+f_1{\left(s^{\ast }\right)}^2\end{array}}\right)+f_1^{\prime }\left(s^{\ast }\right)D_1{x}^{\ast }\right),$$

$${\lambda }_9=-\frac{1}{2D_1}\left(D_1\left(f_1^{\prime }\left(s^{\ast }\right)x^{\ast }+Q+R+V-Z+2\right)-f_1\left(s^{\ast }\right)-\left(\sqrt{\begin{array}{c}{D}_1\left(\begin{array}{c}{Q}^2+2\mathit{QR}+2\mathit{QV}-2\mathit{QZ}-2Q{f}_1\left(s^{\ast }\right)-f_1^{\prime }\left(s^{\ast }\right)x^{\ast }\left(2Q+2R+2V-2Z-f_1^{\prime }\left(s^{\ast }\right)x^{\ast }\right)+R^2+\\ 2\mathit{RV}-2\mathit{RZ}+V^2-2\mathit{VZ}+Z^2\end{array}\right)+\\ f_1\left(s^{\ast }\right)\left(2D_{1}R-2D_{1}V-2D_{1}Z-f_1^{\prime }\left(s^{\ast }\right)x^{\ast }\right)+f_1{\left(s^{\ast }\right)}^2\end{array}}\right)+f_1^{\prime }\left(s^{\ast }\right)D_1{x}^{\ast }\right).$$

These results lead to the following theorem. Theorem 7.

Suppose that $${\lambda }_8<0$$ and $${\lambda }_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\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ and the existence of positive interior equilibrium point $$E_3\left(s^{\ast },x^{\ast },y^{\ast }\right)$$ . For global stability analysis, we use the similar technique as in [3, 18].

(i) Global stability analysis of $$E_2\left({\zeta }_s,\frac{1-{\zeta }_s}{D_1},0\right)$$

Consider system (3) restricted to $${\mathbb{R}}_{\mathit{sx}}^{+}$$ as in (5). Let $$\mathcal{N}$$ be a neighbourhood of equilibrium point $$E_2\left({\zeta }_s,\frac{1-{\zeta }_s}{D_1},0\right)$$ in $${\mathbb{R}}_{\mathit{sx}}^{+}$$ . To analyse the global stability of the equilibrium point $$E_2$$ , a suitable Lyapunov function $$L=\frac{1}{2}\left(n_1{\left(s-\widehat{s}\right)}^2+n_2{\left(x-\widehat{x}\right)}^2\right)$$ is used, where $$\widehat{s}$$ and $$\widehat{x}$$ denote the components of $$E_2\left(\widehat{s},\widehat{x}\right)$$ , i.e., $$\widehat{s}={\zeta }_s$$ and $$\widehat{x}=\frac{1-{\zeta }_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 $${\mathbb{R}}_{\mathit{sx}}^{+}$$ and a Lyapunov function for system (5) in $$\mathcal{N}$$ . By differentiating $$L$$ with respect to time $$t,$$ we get

$$\frac{\mathit{dL}}{\mathit{dt}}=n_1\left(s-\widehat{s}\right)s^{\prime }+n_2\left(x-\widehat{x}\right)x^{\prime },$$ E7

where $${\displaystyle \widehat{x}}=\frac{1-{\displaystyle \widehat{s}}}{D_1}$$ . From (5), we have

$$1=\widehat{s}+\widehat{x}{f}_1\left(\widehat{s}\right)\text{and}{D}_1=f_1\left(\widehat{s}\right).$$

Hence (7) can be written as

$$\begin{array}{l}{L}^{\prime }=n_1\left(s-\widehat{s}\right)\left(\widehat{s}+\widehat{x}{f}_1\left(\widehat{s}\right)-s-f_1\left(s\right)x\right)+n_2\left(x-\widehat{x}\right)\left(\frac{f_1\left(s\right)}{f_1\left(\widehat{s}\right)}-1\right)x\\ =-n_1{\left(s-\widehat{s}\right)}^2+n_1\left(s-\widehat{s}\right)\left(f_1\left(\widehat{s}\right)\widehat{x}-f_1\left(s\right)x\right)+n_2\left(x-\widehat{x}\right)x\left(\frac{f_1\left(s\right)}{f_1\left(\widehat{s}\right)}-1\right)\\ =n_{11}{\left(s-\widehat{s}\right)}^2+\frac{1}{2}{n}_{12}\left(s-\widehat{s}\right)\left(x-\widehat{x}\right)+\frac{1}{2}{n}_{21}\left(s-\widehat{s}\right)\left(x-\widehat{x}\right)+n_{22}{\left(x-\widehat{x}\right)}^2,\end{array}$$

where

$$n_{11}=-n_1<0,$$

$$n_{12}=n_{21}=n_1\frac{f_1\left(\widehat{s}\right)\widehat{x}-f_1\left(s\right)x}{x-\widehat{x}},$$

$$n_{22}=n_2\frac{x\left(f_1\left(s\right)-f_1\left(\widehat{s}\right)\right)}{\left(x-\widehat{x}\right)f_1\left(\widehat{s}\right)}.$$

Clearly that $$L^{\prime }$$ can be written as $$L^{\prime }=X^T\mathit{NX},$$ which $$T$$ denotes the transpose and the matrix $$N$$ is particularly a real, symmetric $$2\times 2$$ matrix, where $$X$$ and $$N$$ can be represented by

$$X=\left(\begin{array}{c}{v}_1\\ v_2\end{array}\right)=\left(\begin{array}{c}s-{\displaystyle \widehat{s}}\\ x-{\displaystyle \widehat{x}}\end{array}\right) \text{and} N=\left(\begin{array}{cc}{n}_{11}& \frac{1}{2}{n}_{12}\\ \frac{1}{2}{n}_{21}& n_{22}\end{array}\right).$$

Thus, it leads to the following theorem. Theorem 8.

The equilibrium point $$E_2$$ is global asymptotically stable with respect to solution trajectories are initiated from $$\text{int}{\mathbb{R}}_{\mathit{sx}}^{+}$$ if the assumptions $$n_{22}<0$$ and $$\text{det}N>0$$ are satisfied.

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

$$n_{22}<0, \text{and} \mathit{det}N=\mathit{det}\left(\begin{array}{cc}{n}_{11}& \frac{1}{2}{n}_{12}\\ \frac{1}{2}{n}_{21}& n_{22}\end{array}\right)>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. Lemma 2.

[19] Let $$G$$ be an isolated hyperbolic equilibrium point in the omega limit set, $$\omega \left(X\right)$$ of the orbit $$\mathcal{O}\left(X\right).$$ Then either $$\omega \left(X\right)=G$$ or there exist points $$J^{+}$$ and $$J^{-}$$ in $$\omega \left(X\right)$$ with $$J^{+}\in W^s\left(G\right)$$ and $$J^{-}\in W^u\left(G\right)$$ , where $$W^s\left(G\right)$$ and $$W^u\left(G\right)$$ denote stable and unstable manifolds of $$G$$ .

Theorem 9.

Assume that

1. $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ 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 $${\mathbb{R}}_{\mathit{sxy}}^{+}$$ are uniformly bounded and attracted into region $${\rm H}$$ (as in Theorem 1), and

3. equilibrium point $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ is globally asymptotically stable with respect to $${\mathbb{R}}_{\mathit{sx}}^{+}$$ .

Then, the system (3) is uniformly persistence.

Proof. This proof strictly depends on Lemma 2. Suppose $${\rm H}$$ is the region as stated in Theorem 1. It showed that region $${\rm H}$$ is positive invariant set and any solutions of system (3) emanating at a point in $${\mathbb{R}}_{+}^3$$ is uniformly bounded. Despites that, the only compact invariant set on $$\partial {\mathbb{R}}_{+}^3$$ is $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ . Let $$P=E_3\left(s^{\ast },x^{\ast },y^{\ast }\right)$$ belongs to the interior of $${\mathbb{R}}_{+}^3$$ , i.e., $$P\in \text{int}{\mathbb{R}}_{+}^3$$ . We shall show that there is no points $$J_i\in \partial {\mathbb{R}}_{+}^3$$ where $$\partial {\mathbb{R}}_{+}^3$$ belongs to $$\mathrm{\omega }\left(P\right)$$ , the omega limit set of $$P.$$ Now, we prove that the equilibrium point $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)\notin \mathrm{\omega }\left(P\right).$$ Assume that $$E_2\in \mathrm{\omega }\left(P\right)$$ is true. Then, there is a point $$J_1^{+}\in W^s\left(E_2\right)\backslash \left\{E_2\right\}$$ such that $$J_1^{+}\in \mathrm{\omega }\left(P\right)$$ by Lemma 2. But $$W^s\left(E_2\right)\cap \left({\mathbb{R}}_{+}^3\backslash \left\{E_2\right\}\right)$$ is empty which is a contradiction for the positive invariance property of $${\rm H}\subset {\mathbb{R}}_{+}^3.$$ Thus, the equilibrium point $$E_2$$ is not in the omega limit set of $$P$$ ; $$E_2\notin \mathrm{\omega }\left(P\right).$$ Next, we shall show $$\partial {\mathbb{R}}_{+}^3\cap \mathrm{\omega }\left(P\right)=\varnothing .$$ Assume that $$\partial {\mathbb{R}}_{+}^3\cap \mathrm{\omega }\left(P\right)\ne \varnothing$$ , and let $$J\in \partial {\mathbb{R}}_{+}^3$$ and $$J\in \mathrm{\omega }\left(P\right).$$ Then the closure of the orbit of the point $$J$$ , i.e., $$\mathrm{cl}\left(O\left(J\right)\right)$$ must either contains $$E_2$$ or unbounded. This is a contradiction, and hence it is proved that $$\partial {\mathbb{R}}_{+}^3\cap \mathrm{\omega }\left(P\right)=\varnothing .$$ We deduce that if $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ is unstable, then, for stable manifold $$W^s\left(E_2\right);W^s\left(E_2\right)\cap \left(\text{int}\left({\mathbb{R}}_{+}^3\right)\right)=\varnothing$$ , and for unstable manifold $$W^u\left(E_2\right);$$ $$W^u\left(E_2\right)\cap \left(\text{int}\left({\mathbb{R}}_{+}^3\right)\right)\ne \varnothing$$ . Therefore, the result of uniform persistence follows since the omega limit set of $$P$$ , $$\mathrm{\omega }\left(P\right)$$ must be in $$\text{int}\left({\mathbb{R}}_{+}^3\right).$$ This completes the proof. $$∎$$ Remark 2.

Global stability of equilibrium point $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0\right)$$ with respect to $${\mathbb{R}}_{\mathit{sx}}^{+}$$ indicates that the boundary flow is isolated and a cyclic with respect to region $${\rm H}$$ . Thus, the system (3) undergoes uniform persistence and implies that a positive interior equilibrium point $$E_3\left(s^{\ast },x^{\ast },y^{\ast }\right)$$ exists (see [20]).

2.5. Hopf bifurcation

In this section, we investigate Hopf bifurcation on the system (3) with a bifurcation real parameter, $$\sigma .$$ Particularly, $$\sigma$$ is selected in such a way that the growth rate function $$f_2$$ is a function of $$x$$ and $$\sigma .$$ Therefore, system (3) takes of the form

$$\left\{\begin{array}{l}{s}^{\prime }=1-s-f_1\left(s\right)x,\\ x^{\prime }=\left(\frac{f_1\left(s\right)}{D_1}-1\right)x-f_2\left(x;\sigma \right)y\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right),\\ y^{\prime }=f_2\left(x;\sigma \right)y\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)-y,\end{array}\right.$$ E8

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

$$\left\{\begin{array}{c}S=s-h_1;\\ X=x-h_2;\\ Y=y-h_3;\end{array}\right.\Rightarrow \left\{\begin{array}{c}s=S+h_1;\\ x=X+h_2;\\ y=Y+h_3;\end{array}\right.$$

where $$\left(h_1,h_2,h_3\right)$$ is the non-trivial equilibrium point. Then, we obtain the following differential equations

$$\left\{\begin{array}{l}{S}^{\prime }=1-S-h_1-f_1\left(S+h_1\right)\left(X+h_2\right),\\ X^{\prime }=\left(\frac{f_1\left(S+h_1\right)}{D_1}-1\right)\left(X+h_2\right)-f_2\left(X+h_2;\sigma \right)\left(Y+h_3\right)\left(1+\frac{D_1\left(X+h_2\right)}{1+\left(X+h_2\right)+D_{1}b\left(Y+h_3\right)}\right),\\ Y^{\prime }=f_2\left(X+h_2;\sigma \right)\left(Y+h_3\right)\left(1+\frac{D_1\left(X+h_2\right)}{1+\left(X+h_2\right)+D_{1}b\left(Y+h_3\right)}\right)-\left(Y+h_3\right).\end{array}\right.$$ E9

The Jacobian matrix of system (8) is given by

$$J_{\sigma }=\left[\begin{array}{cc}-f_1^{\prime }\left(s\right)x-1& -f_1\left(s\right)\\ \frac{f_1^{\prime }\left(s\right)x}{D_1}& \frac{f_1\left(s\right)}{D_1}-f_2\left(x;\sigma \right)y\left[\left(\frac{D_1}{1+x+D_1\mathit{by}}\right)-\left(\frac{D_{1}x}{{\left(1+x+D_1\mathit{by}\right)}^2}\right)\right]-{f_2}^{\prime }\left(x;\sigma \right)y\left(\left(\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)+1\right)-1\\ 0& f_2\left(x;\sigma \right)y\left[\left(\frac{D_1}{1+x+D_1\mathit{by}}\right)-\left(\frac{D_{1}x}{{\left(1+x+D_1\mathit{by}\right)}^2}\right)\right]+f_2^{\prime }\left(x;\sigma \right)y\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)\end{array}\begin{array}{c}0\\ f_2\left(x;\sigma \right)\left[\left(\frac{D_1^2\mathit{bxy}}{{\left(1+x+D_1\mathit{by}\right)}^2}\right)-\left(\frac{D_{1}x}{1+x+D_1\mathit{by}}+1\right)\right]\\ f_2\left(x;\sigma \right)\left[\left(1+\frac{D_{1}x}{1+x+D_1\mathit{by}}\right)-\left(\frac{D_1^2\mathit{bxy}}{{\left(1+x+D_1\mathit{by}\right)}^2}\right)\right]-1\end{array}\right].$$

Thus, the Jacobian matrix of system (8) about $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0;\sigma \right)$$ is

$$J_{\sigma \left(E_2\right)}=\left(\begin{array}{ccc}-\frac{f_1^{\prime }\left({\zeta }_s\right)\left(1-{\zeta }_s\right)}{D_1}-1& -f_1\left({\zeta }_s\right)& 0\\ \frac{f_1^{\prime }\left({\zeta }_s\right)\left(1-{\zeta }_s\right)}{D_1^2}& \frac{f_1\left({\zeta }_s\right)}{D_1}-1& -f_2\left(\frac{1-{\zeta }_s}{D_1};\sigma \right)\left(1+\left(\frac{1-{\zeta }_s}{1+\left(\frac{1-{\zeta }_s}{D_1}\right)}\right)\right)\\ 0& 0& f_2\left(\frac{1-{\zeta }_s}{D_1};\sigma \right)\left(1+\left(\frac{1-{\zeta }_s}{1+\left(\frac{1-{\zeta }_s}{D_1}\right)}\right)\right)-1\end{array}\right).$$ E10

The characteristic equation of (10) is given as

$${\lambda }^3+c_1{\lambda }^2+c_2\lambda +c_3=0,$$ E11

where

$$c_1=\frac{1}{{\mathit{D}}_1\left({\mathit{D}}_1-{\mathrm{\zeta }}_{\mathrm{s}}+1\right)}\left(f_1^{\prime }\left({\mathrm{\zeta }}_{\mathrm{s}}\right)\left(1+{\mathit{D}}_1-2{\mathrm{\zeta }}_{\mathrm{s}}+{\mathrm{\zeta }}_{\mathrm{s}}^2-{\mathit{D}}_1{\mathrm{\zeta }}_{\mathrm{s}}\right)+{\mathrm{f}}_1\left({\mathrm{\zeta }}_{\mathrm{s}}\right)\left({\mathrm{\zeta }}_{\mathrm{s}}-{\mathit{D}}_1-1\right)+{\mathit{D}}_1{\mathrm{f}}_2\left(\frac{1-{\mathrm{\zeta }}_{\mathrm{s}}}{{\mathit{D}}_1};\mathrm{\sigma }\right)\left({\mathrm{\zeta }}_{\mathrm{s}}+{\mathit{D}}_1{\mathrm{\zeta }}_{\mathrm{s}}-2{\mathit{D}}_1-1\right)+3{\mathit{D}}_1\left(1+{\mathit{D}}_1-{\mathrm{\zeta }}_{\mathrm{s}}\right)\right),$$

$$c_2=\frac{1}{D_1\left(D_1-{\zeta }_s+1\right)}\left(f_1^{\prime }\left({\zeta }_s\right)\left(2{\zeta }_s^2+2+2D_1-4{\zeta }_s-2D_1{\zeta }_s-f_2\left(\frac{1-{\zeta }_s}{D_1};\sigma \right)\left(1+2D_1-2{\zeta }_s+{\zeta }_s^2-3D_1{\zeta }_s+D_1{\zeta }_s^2\right)\right)-f_1\left({\zeta }_s\right)\left(2+2D_1-2{\zeta }_s+f_2\left(\frac{1-{\zeta }_s}{D_1};\sigma \right)\left({\zeta }_s-2D_1-1+D_1{\zeta }_s\right)\right)-f_2\left(\frac{1-{\zeta }_s}{D_1};\sigma \right)\left(2D_1+4D_1^2-2D_1{\zeta }_s-2D_1^2{\zeta }_s\right)+3D_1{\zeta }_s-3D_1^2-3D_1\right),$$

$$c_3=\frac{1}{{\mathit{D}}_1\left({\mathit{D}}_1-{\mathrm{\zeta }}_{\mathrm{s}}+1\right)}\left({\mathit{f}}_1^{\prime }\left({\mathrm{\zeta }}_{\mathrm{s}}\right)\left(1-2{\mathrm{\zeta }}_{\mathrm{s}}+{\mathit{D}}_1+{\mathrm{\zeta }}_{\mathrm{s}}^2-{\mathit{D}}_1{\mathrm{\zeta }}_{\mathrm{s}}-{\mathit{f}}_2\left(\frac{1-{\mathrm{\zeta }}_{\mathrm{s}}}{{\mathit{D}}_1};\mathrm{\sigma }\right)\left(1-2{\mathrm{\zeta }}_{\mathrm{s}}+{\mathrm{\zeta }}_{\mathrm{s}}^2-{\mathit{D}}_1+{\mathit{D}}_1{\mathrm{\zeta }}_{\mathrm{s}}^2\right)\right)-{\mathit{f}}_1\left({\mathrm{\zeta }}_{\mathrm{s}}\right)\left({\mathit{f}}_2\left(\frac{1-{\mathrm{\zeta }}_{\mathrm{s}}}{{\mathit{D}}_1};\mathrm{\sigma }\right)\left({\mathrm{\zeta }}_{\mathrm{s}}-1-2{\mathit{D}}_1+{\mathit{D}}_1{\mathrm{\zeta }}_{\mathrm{s}}\right)+1-{\mathrm{\zeta }}_{\mathrm{s}}\right)-{\mathit{f}}_2\left(\frac{1-{\mathrm{\zeta }}_{\mathrm{s}}}{{\mathit{D}}_1};\mathrm{\sigma }\right)\left({\mathit{D}}_1-{\mathit{D}}_1{\mathrm{\zeta }}_{\mathrm{s}}\right)+{\mathit{D}}_1{\mathrm{\zeta }}_{\mathrm{s}}-{\mathit{D}}_1-{\mathit{D}}_1^2\right).$$

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

$$M_1=\left[c_1\right];M_2=\left[\begin{array}{cc}{c}_1& 1\\ c_3& c_2\end{array}\right];M_3=\left[\begin{array}{ccc}{c}_1& 1& 0\\ c_3& c_2& c_1\\ 0& 0& c_3\end{array}\right].$$

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

$$\left.\begin{array}{c}{c}_1>0\\ c_3>0\\ c_1{c}_2-c_3>0.\end{array}\right\}$$ E12

When the assumptions on functional responses $$f_1\left(s\right)=f_1\left({\zeta }_s\right)$$ and $$f_2\left(x\right)=f_2\left(\frac{1-{\zeta }_s}{D_1};\sigma \right)$$ as in (2) for the system (3), together with the hypotheses $$H_1$$ until $$H_3$$ ;

$$H_1:D_1-{\zeta }_s+1>0,$$

$$H_2:f_1^{\prime }\left({\zeta }_s\right)\left(1+D_1-2{\zeta }_s+{\zeta }_s^2-D_1{\zeta }_s\right)+f_1\left({\zeta }_s\right)\left({\zeta }_s-D_1-1\right)+D_1{f}_2\left(\frac{1-{\zeta }_s}{D_1};\sigma \right)\left({\zeta }_s+D_1{\zeta }_s-2D_1-1\right)+3D_1\left(1+D_1-{\zeta }_s\right)>0,$$

$$H_3:f_1^{\prime }\left({\zeta }_s\right)\left(1-2{\zeta }_s+D_1+{\zeta }_s^2-D_1{\zeta }_s-f_2\left(\frac{1-{\zeta }_s}{D_1};\sigma \right)\left(1-2{\zeta }_s+{\zeta }_s^2-D_1+D_1{\zeta }_s^2\right)\right)-f_1\left({\zeta }_s\right)\left(f_2\left(\frac{1-{\zeta }_s}{D_1};\sigma \right)\left({\zeta }_s-1-2D_1+D_1{\zeta }_s\right)+1-{\zeta }_s\right)-f_2\left(\frac{1-{\zeta }_s}{D_1};\sigma \right)\left(D_1-D_1{\zeta }_s\right)+D_1{\zeta }_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 $$\sigma ,$$ say, $${\sigma }_1^{\ast }.$$

Since at $$\sigma ={\sigma }_1^{\ast },$$ there exists an interval $$\left({\sigma }_1^{\ast }-\epsilon ,{\sigma }_1^{\ast }+\epsilon \right)$$ containing $${\sigma }_1^{\ast }$$ for some $$\epsilon >0$$ such that $$\sigma \in \left({\sigma }_1^{\ast }-\epsilon ,{\sigma }_1^{\ast }+\epsilon \right)$$ . Thus, for $$\sigma \in \left({\sigma }_1^{\ast }-\epsilon ,{\sigma }_1^{\ast }+\epsilon \right),$$ the characteristic equation (11) cannot have positive real roots. For $$\sigma ={\sigma }_1^{\ast },$$ we acquire (see [3, 10])

$$\left({\lambda }^2+c_2\right)\left(\lambda +c_1\right)=0,$$ E13

that consist of three roots; $${\lambda }_1=\sqrt{c_2}i,{\lambda }_2=-\sqrt{c_2}i,$$ and $${\lambda }_3=-c_1.$$ As for $$\sigma \in \left({\sigma }_1^{\ast }-\epsilon ,{\sigma }_1^{\ast }+\epsilon \right),$$ all the roots are in general of the form;

$${\lambda }_1\left(\sigma \right)=\alpha \left(\sigma \right)+\beta \left(\sigma \right)i,$$

$${\lambda }_2\left(\sigma \right)=\alpha \left(\sigma \right)-\beta \left(\sigma \right)i,$$

$${\mathrm{\lambda }}_3\left(\mathrm{\sigma }\right)=-{\mathrm{c}}_1\left(\mathrm{\sigma }\right).$$

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

$$\text{Re}{\left[\frac{d{\lambda }_j}{\mathit{d\sigma }}\right]}_{\sigma ={\sigma }_1^{\ast }}\ne 0,j=1,2,3.$$ E14

By substituting $${\lambda }_1\left(\sigma \right)=\alpha \left(\sigma \right)+\beta \left(\sigma \right)i$$ and $${\lambda }_2\left(\sigma \right)=\alpha \left(\sigma \right)-\beta \left(\sigma \right)i$$ into characteristic equation (11) and calculating the implicit derivative, we obtain the following equations

$$\left.\begin{array}{c}K\left(\sigma \right){\alpha }^{\prime }\left(\sigma \right)-L\left(\sigma \right){\beta }^{\prime }\left(\sigma \right)+M\left(\sigma \right)=0,\\ L\left(\sigma \right){\alpha }^{\prime }\left(\sigma \right)+K\left(\sigma \right){\beta }^{\prime }\left(\sigma \right)+N\left(\sigma \right)=0,\end{array}\right\}$$ E15

where

$$K\left(\sigma \right)=3{\alpha }^2\left(\sigma \right)+2c_1\left(\sigma \right)\alpha \left(\sigma \right)+c_2\left(\sigma \right)-3{\beta }^2\left(\sigma \right);$$

$$M\left(\sigma \right)={\alpha }^2\left(\sigma \right)c_1^{\prime }\left(\sigma \right)+c_2^{\prime }\left(\sigma \right)\alpha \left(\sigma \right)-c_1^{\prime }\left(\sigma \right){\beta }^2\left(\sigma \right)+c_3^{\prime }\left(\sigma \right);$$

$$L\left(\sigma \right)=6\alpha \left(\sigma \right)\beta \left(\sigma \right)+2c_1\left(\sigma \right)\beta \left(\sigma \right);$$

$$N\left(\sigma \right)=2\alpha \left(\sigma \right)\beta \left(\sigma \right)c_1^{\prime }\left(\sigma \right)+c_2^{\prime }\left(\sigma \right)\beta \left(\sigma \right).$$

Since $$K\left({\sigma }_1^{\ast }\right)M\left({\sigma }_1^{\ast }\right)+L\left({\sigma }_1^{\ast }\right)N\left({\sigma }_1^{\ast }\right)\ne 0,$$ we have

$$\text{Re}{\left[\frac{d{\lambda }_j}{\mathit{d\sigma }}\right]}_{\sigma ={\sigma }_1^{\ast }}=\frac{K\left({\sigma }_1^{\ast }\right)M\left({\sigma }_1^{\ast }\right)+L\left({\sigma }_1^{\ast }\right)N\left({\sigma }_1^{\ast }\right)}{K^2\left({\sigma }_1^{\ast }\right)+L^2\left({\sigma }_1^{\ast }\right)}\ne 0,j=1,2,3,$$

and $${\lambda }_3\left({\sigma }_1^{\ast }\right)=-c_1\left({\sigma }_1^{\ast }\right)\ne 0.$$ We conclude the details above in the following theorem. Theorem 10.

Suppose that the equilibrium point $$E_2\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0;\sigma \right)$$ 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\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0;\sigma \right)$$ for some values of $$\sigma$$ in the neighbourhood of $${\sigma }_1^{\ast }.$$

Next, we determine the Hopf bifurcation at the equilibrium point $$E_3\left(s^{\ast },x^{\ast },y^{\ast };\sigma \right).$$ The Jacobian matrix of the system (8) about the equilibrium point $$E_3$$ is given by

$$J_{\sigma \left(E_3\right)}=\left[\begin{array}{cc}-f_1^{\prime }\left(s^{\ast }\right)x^{\ast }-1& -f_1\left(s^{\ast }\right)\\ \frac{f_1^{\prime }\left(s^{\ast }\right)x^{\ast }}{{\mathit{D}}_1}& \frac{f_1\left(s^{\ast }\right)}{{\mathit{D}}_1}-f_2\left(x^{\ast };\sigma \right)y^{\ast }\left[\left(\frac{{\mathit{D}}_1}{1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }}\right)-\left(\frac{{\mathit{D}}_1{x}^{\ast }}{{\left(1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }\right)}^2}\right)\right]-{f_2}^{\prime }\left(x^{\ast };\sigma \right)y^{\ast }\left(\frac{{\mathit{D}}_1{x}^{\ast }}{1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }}+1\right)-1\\ 0& f_2\left(x^{\ast };\sigma \right)y^{\ast }\left[\left(\frac{{\mathit{D}}_1}{1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }}\right)-\left(\frac{{\mathit{D}}_1{x}^{\ast }}{{\left(1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }\right)}^2}\right)\right]+f_2^{\prime }\left(x^{\ast };\sigma \right)y^{\ast }\left(\frac{{\mathit{D}}_1{x}^{\ast }}{1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }}+1\right)\end{array},,\begin{array}{c}0\\ f_2\left(x^{\ast };\sigma \right)\left(\frac{{\mathit{D}}_1^{2}b{x}^{\ast }{y}^{\ast }}{{\left(1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }\right)}^2}-\frac{{\mathit{D}}_1{x}^{\ast }}{1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }}-1\right)\\ f_2\left(x^{\ast };\sigma \right)\left(\frac{{\mathit{D}}_1{x}^{\ast }}{1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }}-\frac{{\mathit{D}}_1^{2}b{x}^{\ast }{y}^{\ast }}{{\left(1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }\right)}^2}+1\right)-1\end{array}\right].$$

Hence, the characteristic equation for the Jacobian matrix $$J_{\sigma \left(E_3\right)}$$ is

$${\lambda }^3+c_4{\lambda }^2+c_5\lambda +c_6=0,$$ E16

where

$$c_4=f_1^{\prime }\left(s^{\ast }\right)x^{\ast }+Q+R+V-Z+3-\frac{f_1\left(s^{\ast }\right)}{D_1},$$

$$c_5=f_1^{\prime }\left(s^{\ast }\right)x^{\ast }\left(Q+R+V-Z+2\right)+\frac{f_1\left(s^{\ast }\right)\left(Z-R-2\right)}{D_1}+2Q+2R+2V-2Z+3,$$

$$c_6=f_1^{\prime }\left(s^{\ast }\right)x^{\ast }\left(Q+R+V-Z+1\right)+\frac{f_1\left(s^{\ast }\right)\left(Z-R-1\right)}{D_1}+Q+R+V-Z+1,$$

and

$$Q=f_2\left(x^{\ast };\sigma \right)y^{\ast }\left[\left(\frac{D_1}{1+x^{\ast }+D_{1}b{y}^{\ast }}\right)-\left(\frac{D_1{x}^{\ast }}{{\left(1+x^{\ast }+D_{1}b{y}^{\ast }\right)}^2}\right)\right],$$

$$R=\left(\frac{f_2\left(x^{\ast };\sigma \right)D_1^{2}b{x}^{\ast }{y}^{\ast }}{{\left(1+x^{\ast }+D_{1}b{y}^{\ast }\right)}^2}\right),$$

$$V=f_2^{\prime }\left(x^{\ast };\sigma \right)y^{\ast }\left(\frac{D_1{x}^{\ast }}{1+x^{\ast }+D_{1}b{y}^{\ast }}+1\right),$$

$$Z=f_2\left(x^{\ast };\sigma \right)\left(\frac{{\mathit{D}}_1{x}^{\ast }}{1+x^{\ast }+{\mathit{D}}_{1}b{y}^{\ast }}+1\right).$$

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

$$M_4=\left[c_4\right];$$ $$M_5=\left[\begin{array}{cc}{c}_4& 1\\ c_6& c_5\end{array}\right];$$ $$M_6=\left[\begin{array}{ccc}{c}_4& 1& 0\\ c_6& c_5& c_4\\ 0& 0& c_6\end{array}\right].$$

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

$$\left.\begin{array}{c}{c}_4>0\\ c_6>0\\ c_4{c}_5-c_6>0.\end{array}\right\}$$ E17

Suppose the assumptions of functional response $$f_1\left(s^{\ast }\right)$$ and $$f_2\left(x^{\ast };\sigma \right)$$ similar as in (2) for the system (3), together with the hypotheses $$H_4$$ until $$H_6$$ ;

$$H_4:f_1^{\prime }\left(s^{\ast }\right)x^{\ast }+Q+R+V+3>\frac{f_1\left(s^{\ast }\right)}{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 $$\sigma ,$$ say, $${\sigma }_2^{\ast }.$$ Since at $$\sigma ={\sigma }_2^{\ast },$$ there exists an interval $$\left({\sigma }_2^{\ast }-\epsilon ,{\sigma }_2^{\ast }+\epsilon \right)$$ containing $${\sigma }_2^{\ast }$$ for some $$\epsilon >0$$ . Then, for $$\sigma \in \left({\sigma }_2^{\ast }-\epsilon ,{\sigma }_2^{\ast }+\epsilon \right)$$ , the characteristic Eq. (16) cannot have positive real roots. For $$\sigma ={\sigma }_2^{\ast },$$ we get (see [3, 10])

$$\left({\lambda }^2+c_5\right)\left(\lambda +c_4\right)=0,$$ E18

Consisting of three roots; $${\lambda }_1=\sqrt{c_5}i,{\lambda }_2=-\sqrt{c_5}i,{\lambda }_3=-c_4.$$ As for $$\sigma \in \left({\sigma }_2^{\ast }-\epsilon ,{\sigma }_2^{\ast }+\epsilon \right),$$ all roots are in general of the form;

$${\lambda }_1\left(\sigma \right)=\alpha \left(\sigma \right)+\beta \left(\sigma \right)i,$$

$${\lambda }_2\left(\sigma \right)=\alpha \left(\sigma \right)-\beta \left(\sigma \right)i,$$

$${\lambda }_3\left(\sigma \right)=-c_4\left(\sigma \right),$$

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

$$\mathit{Re}{\left[\frac{d{\lambda }_j}{\mathit{d\sigma }}\right]}_{\sigma ={\sigma }_2^{\ast }}\ne 0,j=1,2,3.$$ E19

By substituting $${\lambda }_1\left(\sigma \right)=\alpha \left(\sigma \right)+\beta \left(\sigma \right)i$$ and $${\lambda }_2\left(\sigma \right)=\alpha \left(\sigma \right)-\beta \left(\sigma \right)i$$ into characteristic equation (18) and calculating the implicit derivative, we get the following equations;

$$\left.\begin{array}{c}{A}_1\left(\sigma \right){\alpha }^{\prime }\left(\sigma \right)-A_2\left(\sigma \right){\beta }^{\prime }\left(\sigma \right)+B_1\left(\sigma \right)=0,\\ A_2\left(\sigma \right){\alpha }^{\prime }\left(\sigma \right)+A_1\left(\sigma \right){\beta }^{\prime }\left(\sigma \right)+B_2\left(\sigma \right)=0,\end{array}\right\}$$ E20

where

$$A_1\left(\sigma \right)=3{\alpha }^2\left(\sigma \right)+2c_1\left(\sigma \right)\alpha \left(\sigma \right)+c_2\left(\sigma \right)-3{\beta }^2\left(\sigma \right);$$

$$A_2\left(\sigma \right)=6\alpha \left(\sigma \right)\beta \left(\sigma \right)+2c_1\left(\sigma \right)\beta \left(\sigma \right);$$

$$B_1\left(\sigma \right)={\alpha }^2\left(\sigma \right)c_1^{\prime }\left(\sigma \right)+c_2^{\prime }\left(\sigma \right)\alpha \left(\sigma \right)-c_1^{\prime }\left(\sigma \right){\beta }^2\left(\sigma \right)+c_3^{\prime }\left(\sigma \right);$$

$$B_2\left(\sigma \right)=2\alpha \left(\sigma \right)\beta \left(\sigma \right)c_1^{\prime }\left(\sigma \right)+c_2^{\prime }\left(\sigma \right)\beta \left(\sigma \right).$$

Since $$A_1\left({\sigma }_2^{\ast }\right)B_1\left({\sigma }_2^{\ast }\right)+A_2\left({\sigma }_2^{\ast }\right)B_2\left({\sigma }_2^{\ast }\right)\ne 0,$$ we have

$$\text{Re}{\left[\frac{d{\lambda }_j}{\mathit{d\sigma }}\right]}_{\sigma ={\sigma }_2^{\ast }}=\frac{A_1\left({\sigma }_2^{\ast }\right)B_1\left({\sigma }_2^{\ast }\right)+A_2\left({\sigma }_2^{\ast }\right)B_2\left({\sigma }_2^{\ast }\right)}{{\left(A_1\left({\sigma }_2^{\ast }\right)\right)}^2+{\left(A_2\left({\sigma }_2^{\ast }\right)\right)}^2}\ne 0,j=1,2,3,$$

and $${\lambda }_3\left({\sigma }_2^{\ast }\right)=-c_4\left({\sigma }_2^{\ast }\right)\ne 0.$$ We summarise the discussion above in the following theorem. Theorem 11.

Suppose that the equilibrium point

$$E_3\left(s^{\ast },x^{\ast },y^{\ast };\sigma \right)=E_3\left(s^{\ast },\frac{1-s^{\ast }}{f_1\left(s^{\ast }\right)},\frac{\left(2x^{\ast }+1\right)f_2\left(x^{\ast }\right)-x^{\ast }-1}{b{D}_1\left(1-f_2\left(x^{\ast }\right)\right)};\sigma \right)$$

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\left(s^{\ast },x^{\ast },y^{\ast };\sigma \right)$$ for some values of $$\sigma$$ in the neighbourhood of $${\sigma }_2^{\ast }.$$

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\left(1,0,0\right)$$ (no prey and predator organisms present) is locally asymptotically stable if $${\lambda }_3=\frac{f_1\left(1\right)}{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\left({\zeta }_s,\frac{\left(1-{\zeta }_s\right)}{D_1},0;\sigma \right)$$ and $$E_3\left(s^{\ast },x^{\ast },y^{\ast };\sigma \right)$$ for some values of $$\sigma$$ in the neighbourhoods of $${\sigma }_1^{\ast }$$ and $${\sigma }_2^{\ast }$$ , respectively. The method used to obtain these results is similar to the method used in [3, 10].