Mathematical Analysis of a Prey–Predator System: An Adaptive Back-Stepping Control and Stochastic Approach

Abstract

In this paper, stochastic analysis of a diseased prey–predator system involving adaptive back-stepping control is studied. The system was investigated for its dynamical behaviours, such as boundedness and local stability analysis. The global stability of the system was derived using the Lyapunov function. The uniform persistence condition for the system is obtained. The proposed system was studied with adaptive back-stepping control, and it is proved that the system stabilizes to its steady state in nonlinear feedback control. The value of the system is described mostly by the environmental stochasticity in the form of Gaussian white noise. We also established some conditions for oscillations of all positive solutions of the delayed system. Numerical simulations are illustrated, and sustained our analytical findings. We concluded that controlled harvesting on the susceptible and infected prey is able to control prey infection.

1. Introduction

Theoretical research and field observations have established the prevalence of various infectious diseases amongst the majority of the ecosystem population. In the ecological system, the impact of such infectious diseases is an important area of research for ecologists and mathematicians. The processes of merging ecology and epidemiology in the past few decades have been challenging and interesting. By nature, species are always dependent on other species for its food and living space. It is responsible for spreading infectious diseases and also competes against and is predated by other species. The dynamical behavior of such systems is analyzed using mathematical models that are described by differential equations. Mathematical epidemic models have gained much attention from researchers after the pioneering work of Kermack and McKendrick [1] on the SIRS (Susceptible-Infective-Removal-Susceptible) system, in which the evolution of a disease which gets transmitted upon contact is described. The influence of epidemics on predation was first studied by Anderson and May [2,3]. Hadler and Freedman [4] considered the prey–predator model in which predation is more likely on infected prey. In their model, they considered that predators only became infected from infected prey by predation. Haque and Venturino [5,6] discussed the models of diseases spreading in symbolic communities. Mukhopadhyay [7] studied the role of harvesting and switching on the dynamics of disease propagation and/or eradication. The role of prey infection on the stability aspects of a prey–predator models with different functional responses was studied by Bairagi et al. [8]. Han et al. [9] analyzed four epidemiological models for SIS (Susceptible-Infectious-Susceptible) and SIR (Susceptible-Infectious-Recovered) diseases with standard and mass action incidents. Das [10] showed that parasite infection in predator populations stabilized prey–predator oscillations. Pal and Samanta [11] studied the dynamical behavior of a ratio-dependent prey–predator system with infection in the prey population. They proved that prey refuge had a stabilizing effect on the prey–predator interaction. Numerous examples of the prey–predator relationship with infection in the prey population have been found in various studies [12,13,14,15,16,17].

Adaptation is a fundamental characteristic of living organisms, such as in the prey–predator system and other ecological models, since they attempt to maintain physiological equilibrium in the midst of changing environmental conditions. Adaptive control is an active field in the design of control systems that helps deal with uncertainties. Back-stepping is a technique for designing stability controls for a nonlinear dynamical system, and this approach is a recursive method for stabilizing the origin of the system. The control process terminates when the final external control reaches. El-Gohary and Al-Ruzaiza [18,19] discussed the chaos and adaptive control of a three-species, continuous-time prey–predator model. Recently, Madhusudanan et al. [20] studied back-stepping control in a diseased prey–predator system. They proved that the system was globally asymptotically stable at the origin with the help of nonlinear feedback control. Numerous examples of control techniques in the prey–predator system have been found in various studies [21,22].

The rest of the paper is structured as follows. In Section 2, we formulate a mathematical model with an assumption, and the positivity and boundedness of the deterministic model is also discussed. Section 3 deals with the existence of equilibrium points with a feasible condition. In Section 4, local stability analysis of equilibrium points is discussed. Section 5 deals with global stability analysis of the interior equilibrium point $E_3(x^{*},y^{*},z^{*})$. We discuss the condition for permanence of the system in Section 6. In Section 7, we introduce adaptive back-stepping control in the prey–predator system. In Section 8, we compute the population intensity of fluctuation due to the incorporation of noise, which leads to chaos in reality. In Section 9, we propose and analyze a delayed prey–predator system. Numerical simulation of the proposed model is presented in Section 10. Finally, the discussion is presented in Section 11 and conclusions are presented in the final section.

2. Mathematical Model

In this section, a continuous-time prey–predator system with susceptible, infected prey and a predator is considered. It is assumed that the susceptible prey population was developed on the basis of logistic law, and that only infected prey are predated. The disease is inherited only from the prey population, and they remain infected and do not recover.

The model becomes:

$${\displaystyle \frac{dx}{dt}}=x\left(1-{\displaystyle \frac{x}{K}}\right)-{\displaystyle \frac{axy}{1+x}}-h_{1}x,$$

$${\displaystyle \frac{dy}{dt}}={\displaystyle \frac{axy}{1+x}}-byz-h_{2}y,$$

(1)

$${\displaystyle \frac{dz}{dt}}=fyz-dz.$$

Here, the parameters $x(t)$, $y(t)$, and $z(t)$ denote the susceptible prey, infected prey, and predator populations, respectively. The parameters a, d, b, f, $h_1$, and $h_2$ denote the rate of transmission from the susceptible to infected prey population, death rate of predators, searching efficiency of the predator, conversion-efficiency rate of the predator, and constant harvesting rate of susceptible prey and infected prey, respectively. Now, we will analyze system (1) with the following initial conditions:

$$x(0)\ge 0,y(0)\ge 0,z(0)\ge 0.$$

(2)

Positiveness and Boundedness of the System

In theoretical eco-epidemiology, boundedness of the system implies that the system is biologically well-behaved. The following theorems ensure the positivity and boundedness of the system (1):

Theorem 1.

All solutions of $(x(t),y(t),z(t))$ of system (1) with the initial condition (2) are positive for all $t\ge 0.$

Proof. 

From (1), it is observed that

$${\displaystyle \frac{dx}{x}}=\left[\left(1-{\displaystyle \frac{x}{K}}\right)-{\displaystyle \frac{ay}{1+x}}-h_1\right]dt={\varphi }_1(x,y)dt,$$

where ${\varphi }_1(x,y)=\left(1-{\displaystyle \frac{x}{K}}\right)-{\displaystyle \frac{ay}{1+x}}-h_1$.

Integrating in the region $[0,t]$, we get $x(t)=x(0)exp\left(\int {\varphi }_1(x,y)dt\right)>0$ for all t. From (1), it is observed that

$${\displaystyle \frac{dy}{y}}=\left[{\displaystyle \frac{ax}{1+x}}-bz-h_2\right]dt={\varphi }_2(x,z)dt,$$

where ${\varphi }_2(x,z)={\displaystyle \frac{ax}{1+x}}-bz-h_2.$

Integrating in region $[0,t],$ we get $y(t)=y(0)exp\left(\int {\varphi }_2(x,z)dt\right)>0$ for all t. From (1), it is observed that

$${\displaystyle \frac{dz}{z}}=[fy-d]dt={\varphi }_3(y)dt,$$

where ${\varphi }_3(y)=fy-d$.

Integrating in the region $[0,t]$, we get $z(t)=z(0)exp\left(\int {\varphi }_3(y)dt\right)>0$ for all t. Hence, all solutions starting from interior of the first octant $(In{\Re }_{+}^3)$ remain positive for the future. □

Theorem 2.

All the non-negative solutions of the model system (1) that initiate in ${\Re }_{+}^3$ are uniformly bounded.

Proof. 

Let $x(t),y(t),z(t)$ be any solution of system (1). Since from (1) ${\displaystyle \frac{dx}{dt}}\le x\left(1-{\displaystyle \frac{x}{K}}\right),$ we have $\underset{t\to \infty }{lim}supx(t)\le K.$ Let $\xi =x+y+{\displaystyle \frac{b}{f}}z;$ therefore,

$${\displaystyle \frac{d\xi }{dt}}={\displaystyle \frac{dx}{dt}}+{\displaystyle \frac{dy}{dt}}+{\displaystyle \frac{b}{f}}{\displaystyle \frac{dz}{dt}}.$$

(3)

Substituting Equation (1) in Equation (3), we get

$${\displaystyle \frac{d\xi }{dt}}+m\xi =x\left((1+m-h_1)-{\displaystyle \frac{x}{K}}\right)+(m-h_2)y+{\displaystyle \frac{bz}{f}}(m-d)\le x\left((1+m-h_1)-{\displaystyle \frac{x}{K}}\right),$$

$${\displaystyle \frac{d\xi }{dt}}+m\xi \le \mu \mathrm{since}K(1+m-h_1)=\mu ,$$

where m and $\mu$ are positive constants. Applying Lemma on differential inequalities [23], we obtain

$$0\le \xi (x,y,z)\le (\mu /m)(1-e^{-mt})+(\xi (x(0),y(0),z(0))/e^{mt})$$

and, for $t\to \infty$, we have $0\le \xi (x,y,z)\le (\mu /m).$ Thus, all solutions of system (1) enter into the region

$$\Gamma =\left\{(x,y,z)\in {\Re }_{+}^3:0\le x\le K,0\le \xi \le (\mu /m)+ϵ,\forall ϵ>0\right\}.$$

(4)

This completes the proof. □

3. Equilibrium Points and Criteria for Existence

The possible steady states for system (1) and their existence conditions for each of them are as follows:

• The vanishing equilibrium point, $E_0=(0,0,0)$, always exists.
• The disease and predator free equilibrium point, $E_1=(\tilde{x},0,0)$, where $\tilde{x}=K(1-h_1)$, exists provided that $h_1<1.$
• The predator free equilibrium point, $E_2=(\widehat{x},\widehat{y},0)$, where $\widehat{x}=(h_2/(a-h_2))$ and $\widehat{y}=[(K(a-h_2)(1-h_1)-h_2)/({(a-h_2)}^{2}K)]$, exists provided that $a>h_2,h_1<1,K(a-h_2)(1-h_1)>h_2.$
• The steady state, $E_3=(x^{*},y^{*},z^{*})$, where $y^{*}=d/f$ and $z^{*}=[((a-h_2)x^{*}-h_2)/(b(1+x^{*}))]$, exists.

However, $x^{*}$ is a positive root of (5)

$$A{x^{*}}^2+B{x}^{*}+C=0,$$

(5)

where $A=f,B=Kf{h}_1+f-Kf,C=Kf{h}_1+adK-Kf$.

By Discard’s rule of sign, Equation (5) has a unique positive root, if $K{h}_1+1<K,f{h}_1+ad>f.$

4. Stability Analysis

In this section, we analyzed the local stability of system (1) that is examined by constructing the Jacobian matrix relating to every equilibrium point. The Jacobian matrix of the system at any point $(x,y,z)$ is given by

$$J(x,y,z)=\left(\begin{array}{ccc}1-{\displaystyle \frac{2x}{K}}-h_1-{\displaystyle \frac{ay}{{(1+x)}^2}}& {\displaystyle \frac{-ax}{(1+x)}}& 0\\ {\displaystyle \frac{ay}{{(1+x)}^2}}& {\displaystyle \frac{ax}{1+x}}-bz-h_2& -by\\ 0& fz& fy-d\end{array}\right).$$

• The variational matrix for the equilibrium point at $E_0(0,0,0)$ is

  $$J(E_0)=\left(\begin{array}{ccc}1-h_1& 0& 0\\ 0& -h_2& 0\\ 0& 0& -d\end{array}\right).$$
  The eigenvalues of $J(E_0)$ are ${\lambda }_1=1-h_1,{\lambda }_2=-h_2,{\lambda }_3=-d.$ All the eigenvalues are negative if $h_1>1$. Hence, $E_0$ is locally asymptotically stable in the x–y–z direction.
• The variational matrix for the equilibrium point at $E_1(K(1-h_1),0,0)$ is

  $$J(E_1)=\left(\begin{array}{ccc}{h}_1-1& {\displaystyle \frac{-aK(1-h_1)}{1+K(1-h_1)}}& 0\\ 0& {\displaystyle \frac{aK(1-h_1)}{1+K(1-h_1)}}-h_2& 0\\ 0& 0& -d\end{array}\right).$$
  The eigenvalues of $J(E_1)$ are ${\lambda }_1=h_1-1,{\lambda }_2={\displaystyle \frac{aK(1-h_1)}{1+K(1-h_1)}}-h_2,$ ${\lambda }_3=-d.$ All the eigenvalues are negative if $h_1<1\mathrm{and}aK(1-h_1)<h_2(1+K(1-h_1))$. Hence, $E_1$ is asymptotically stable in the x–y–z direction.
• The variational matrix for the equilibrium point at $E_2(\widehat{x},\widehat{y},0)$ is

  $$J(E_2)=\left(\begin{array}{ccc}1-{\displaystyle \frac{2\widehat{x}}{K}}-{\displaystyle \frac{a\widehat{y}}{{(1+\widehat{x})}^2}}-h_1& -h_2& 0\\ {\displaystyle \frac{K(a-h_2)(1-h_1)-h_2}{Ka}}& 0& -b\widehat{y}\\ 0& 0& f\widehat{y}-d\end{array}\right).$$
  The eigenvalues of $J(E_3)$ are negative if it satisfies $f\widehat{y}<d$ and $K(1-h_1){(1+\widehat{x})}^2<(2\widehat{x}{(1+\widehat{x})}^2+aK\widehat{y})$. Hence, the equilibrium point $E_2$ is locally asymptotically stable in the x–y–z direction.

Theorem 3.

The co-existent equilibrium point $E_3(x^{*},y^{*},z^{*})$ of system (1) exists. Then, $E_3$ is locally asymptotically stable if the following conditions satisfy

1. 

$a{x}^{*}<(b{z}^{*}+h_2)(1+x^{*})$ and

2. 

$K(1-h_1){(1+x^{*})}^2<(2x^{*}{(1+x^{*})}^2+aK{y}^{*}).$

Proof. 

The variational matrix at the interior point $E_3(x^{*},y^{*},z^{*})$ is

$$J(E_3)=\left(\begin{array}{ccc}{a}_{11}& a_{12}& 0\\ a_{21}& a_{22}& a_{23}\\ 0& a_{32}& 0\end{array}\right),$$

where

$$\begin{array}{c}{a}_{11}=1-{\displaystyle \frac{2x^{*}}{K}}-h_1-{\displaystyle \frac{a{y}^{*}}{{(1+x^{*})}^2}}{a}_{12}={\displaystyle \frac{-a{x}^{*}}{1+x^{*}}}{a}_{21}={\displaystyle \frac{a{y}^{*}}{{(1+x^{*})}^2}},\hfill \\ a_{22}={\displaystyle \frac{a{x}^{*}}{1+x^{*}}}-b{z}^{*}-h_2{a}_{23}=-b{y}^{*}{a}_{32}=f{z}^{*}.\hfill \end{array}$$

The characteristic equations of Jacobian matrix $J(E_3)$ is given by ${\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3=0,$ where

$$A_1=-a_{11}-a_{22},A_2=a_{11}{a}_{22}-a_{12}{a}_{21}-a_{23}{a}_{32},A_3=a_{11}{a}_{23}{a}_{32}$$

and

$$A_1{A}_2-A_3=a_{11}{a}_{12}{a}_{21}+a_{12}{a}_{21}{a}_{22}+a_{22}{a}_{23}{a}_{32}-a_{11}^2{a}_{22}-a_{11}{a}_{22}^2.$$

The sufficient conditions for $A_1{A}_2-A_3>0$ are $a_{11}\le 0,$ $a_{22}\le 0$, which implies

$$a{x}^{*}<(b{z}^{*}+hz)(1+x^{*}),K(1-h_1){(1+x^{*})}^2<[2x^{*}{(1+x^{*})}^2+aK{y}^{*}].$$

□

5. Global Stability Analysis

In this section, we investigated the global stability behavior of the system (1) at the interior equilibrium $E_3(x^{*},y^{*},z^{*})$ by using the Lyapunov stability theorem.

Theorem 4.

If $(y{x}^{*}/y^{*})<x<x^{*}$ or $x^{*}<x<(y{x}^{*}/y^{*}),$ then $E_3(x^{*},y^{*},z^{*})$ is globally asymptotically stable.

Proof. 

Let us define

$$\begin{array}{c}\hfill V(x,y,z)=P\left(x-x^{*}-x^{*}ln\left({\displaystyle \frac{x}{x^{*}}}\right)\right)+\left(y-y^{*}-y^{*}ln\left({\displaystyle \frac{y}{y^{*}}}\right)\right)\\ \hfill +Q\left(z-z^{*}-z^{*}ln\left({\displaystyle \frac{z}{z^{*}}}\right)\right),\end{array}$$

(6)

where $P,Q$ are positive constants to be chosen later.

Differentiating (6) along the solution of the system (1) with respect to $t,$ we get

$$\begin{array}{cc}\hfill {\displaystyle \frac{dV}{dt}}& =P\left({\displaystyle \frac{x-x^{*}}{x}}\right){\displaystyle \frac{dx}{dt}}+\left({\displaystyle \frac{y-y^{*}}{y}}\right){\displaystyle \frac{dy}{dt}}+Q\left({\displaystyle \frac{z-z^{*}}{z}}\right){\displaystyle \frac{dz}{dt}}\hfill \\ & =P\left(1-{\displaystyle \frac{x}{K}}-{\displaystyle \frac{ay}{1+x}}-h_1\right)(x-x^{*})+(y-y^{*})\left({\displaystyle \frac{ax}{1+x}}-bz-h_2\right)+Q(z-z^{*})(fy-d)\hfill \\ & =P\left(-{\displaystyle \frac{1}{K}}(x-x^{*})-a\left({\displaystyle \frac{y}{1+x}}-{\displaystyle \frac{y^{*}}{1+x}}\right)\right)(x-x^{*})\hfill \\ & +(y-y^{*})\left({\displaystyle \frac{ax}{1+x}}-{\displaystyle \frac{a{x}^{*}}{1+x^{*}}}-b(z-z^{*})\right)+Q(f(z-z^{*})(y-y^{*})).\hfill \end{array}$$

Choosing $P=1,Q=b/f,$ and then simplified to

$${\displaystyle \frac{dV}{dt}}=-{\displaystyle \frac{1}{K}}{(x-x^{*})}^2-a{\displaystyle \frac{(y{x}^{*}-x{y}^{*})}{(1+x)(1+x^{*})}}(x-x^{*}).$$

Now,

$${\displaystyle \frac{dV}{dt}}<0\mathrm{if}{x}^{*}<x<{\displaystyle \frac{y{x}^{*}}{y^{*}}}\mathrm{or}{\displaystyle \frac{y{x}^{*}}{y^{*}}}<x<x^{*}.$$

Then, ${\displaystyle \frac{dV}{dt}}$ is negative definite. Consequently, V is a Lyapunov function with respect to all solutions in the interior of the positive octant. □

6. Permanence of the System

In this section, our main intuition is that the long time survival of species in an ecological system. Many notions and terms are identified in the literature to discuss and analyze the long-term survival of populations. Out of such, permanence and persistence are the ones to better analyze the system. From an ecological point of view, permanence of a system means that the long-term survival of all populations of the system.

Definition 1.

The system (1) is said to be permanent if ∃$M\ge m>0$, such that for any solution of $(x(t),y(t),z(t))$ of system (1), $(x(0),y(0),z(0))>0$,

$$m\le \underset{t\to \infty }{lim}(inf(x(t))\le \underset{t\to \infty }{lim}(sup(x(t))\le M,$$

$$m\le \underset{t\to \infty }{lim}(inf(y(t))\le \underset{t\to \infty }{lim}(sup(y(t))\le M,$$

$$m\le \underset{t\to \infty }{lim}(inf(z(t))\le \underset{t\to \infty }{lim}(sup(z(t))\le M.$$

Now, we show that system (1) is uniformly persistent. The survival of all populations of the system in the future time is nothing but persistence in the view of ecology.

In the mathematical point of view, persistence of a system means that a strictly positive solution does not have omega limit points on the boundary of the non-negative cone.

Definition 2.

A population is said to be uniformly persistent if there exists $\delta >0$, independent of $x(0)$ where $x(0)>0$, such that

$$\underset{t\to \infty }{lim}\left(inf(x(t))\right)>\delta .$$

Theorem 5.

The system (1) is uniformly persistent if

$$f\widehat{y}-d>0\mathit{holds}.$$

(7)

Proof. 

We will prove this theorem by the method of Lyapunov average function.

Let the average Lyapunov function for the system (1) be $\sigma (X)=x^p{y}^q{z}^r$, where p, q, r are positive constants. Clearly, $\sigma (X)$ is non-negative function defined in D of ${\Re }_{+}^3$, where

$${\Re }_{+}^3=\{(x,y,z),x>0,y>0,{(1+x)}^2-ayK>0\}.$$

Then, we have $\mathrm{\Psi }(X)={\displaystyle \frac{\dot{\sigma }(X)}{\sigma (X)}}=p{\displaystyle \frac{\dot{x}}{x}}+q{\displaystyle \frac{\dot{y}}{y}}+r{\displaystyle \frac{\dot{z}}{z}},$

$$\mathrm{\Psi }(X)=p\left(\left(1-{\displaystyle \frac{x}{K}}\right)-{\displaystyle \frac{ay}{1+x}}-h_1\right)+q\left({\displaystyle \frac{ax}{1+x}}-bz-h_2\right)+r(fy-d).$$

(8)

Furthermore, there are no periodic orbits in the interior of positive quadrant of x–y plane. Thus, to prove the uniform persistence of the system, it is enough to show that $\mathrm{\Psi }(X)>0$ in ${\Re }_{+}^3$ for a suitable choice of $p,q,r>0:$

$$\mathrm{\Psi }(E_0)=p(1-h_1)-q(h_2)-rd>0,$$

(9)

$$\mathrm{\Psi }(E_1)=q\left({\displaystyle \frac{aK(1-h_1)}{1+K(1-h_1)}}-h_2\right)-rd>0,$$

(10)

$$\mathrm{\Psi }(E_2)=r(f\widehat{y}-d)>0.$$

(11)

We noticed that, by increasing the value of p, while $h_1<1$, $p(1-h_1)>(q{h}_2+rd)$, $\mathrm{\Psi }(E_0)$ can be made positive. Thus, the inequality (9) holds. If $qaK(1-h_1)>(q{h}_2+rd)(1+K(1-h_1))$, then $\mathrm{\Psi }(E_1)$ is positive. Thus, the inequality (10) holds. If the inequality in Equation (7) holds, then (11) is satisfied. □

7. Introduction of Adaptive Back-Stepping Control in a Prey–Predator System

Adaptive back-stepping method is the back-stepping control and adaptive laws. The back-stepping design is initiated with the first state equation whose state variable has the highest integration order from the control input. The second state variable is considered as the virtual control and the stabilizing function is replaced by it. This stabilizing function can stabilize the first state variables and we set the error between the virtual control and stabilizing function as $\eta .$ Then, for the second state equation, we will design a new stabilizing law to replace the third state variable for the second order system, and step back to control the signal. From the steps above, we can see that the term back-stepping means that we use the latter state as a virtual control to stabilize the previous one. The Lyapunov direct method is utilized from the stabilization method for the error between virtual control and stabilizing function. The control Lyapunov function is to be used which will be a positive definite and includes the quadratic form of the errors. In this section, the system with susceptible prey, infected prey and a predator population controlled by back-stepping using a nonlinear feedback control approach is studied. We initiate the study by assuming that system (1) can be written in the suitable form by introducing nonlinear feedback control inputs $u_1,u_2,u_3$ into the system to better analyze the prey–predator interactions:

$${\displaystyle \frac{dx}{dt}}=x\left(1-{\displaystyle \frac{x}{K}}\right)-{\displaystyle \frac{axy}{1+x}}-h_{1}x+u_1,$$

(12)

$${\displaystyle \frac{dy}{dt}}={\displaystyle \frac{axy}{1+x}}-byz-h_{2}y+u_2,$$

(13)

$${\displaystyle \frac{dz}{dt}}=fyz-dz+u_3,$$

(14)

where $u_1,u_2,u_3$ are back-stepping nonlinear feedback controllers that are the functions of state variables and will be suitable choices to make the trajectories of the whole system (12)–(14). As long as this feedback stabilizes the system (1), $\underset{t\to \infty }{lim}||x(t)||=0$ converges.

Theorem 6.

A diseased prey–predator system (12)–(14) is globally asymptotically stable provided the following adaptive nonlinear controls

$$u_1={\displaystyle \frac{x^2}{K}}-2x+{\widehat{h}}_{1}x,$$

(15)

$$u_2={\displaystyle \frac{a{x}^2}{1+x}}-{\displaystyle \frac{ax{\eta }_1}{1+x}}+{\widehat{h}}_2{\eta }_1-{\eta }_1,$$

(16)

$$u_3=(b{\eta }_1^2-f{\eta }_1{\eta }_2+d{\eta }_2-{\eta }_2),$$

(17)

with the errors

$${\eta }_1=y-{\nu }_1(x,y),$$

(18)

$${\eta }_2=z-{\nu }_2(x,y,z).$$

(19)

Proof. 

Consider the parameter estimators

$$e_a=a-\widehat{a},e_b=b-\widehat{b},e_{h_1}=h_1-{\widehat{h}}_1,e_{h_2}=h_2-{\widehat{h}}_2,e_d=d-\widehat{d},e_f=f-\widehat{f}.$$

(20)

Considering Equation (12), the Lyapunov function of x is given by

$$V_1(x,e_a,e_{h_1})={\displaystyle \frac{1}{2}}{x}^2+{\displaystyle \frac{1}{2}}{e}_a^2+{\displaystyle \frac{1}{2}}{e}_{h_1}^2.$$

(21)

By using the derivative of (20),

$${\dot{e}}_a=-\dot{\widehat{a}},{\dot{e}}_b=-\dot{\widehat{b}},{\dot{e}}_{h_1}=-\dot{\widehat{h_1}},{\dot{e}}_{h_2}=-\dot{\widehat{h_2}},{\dot{e}}_d=-\dot{\widehat{d}},{\dot{e}}_f=-\dot{\widehat{f}}.$$

Differentiating (21) with respect to $t,$ we have

$${\dot{V}}_1=x\left[x\left(1-{\displaystyle \frac{x}{K}}\right)-{\displaystyle \frac{axy}{1+x}}-h_{1}x+u_1\right]+e_a(-\dot{\widehat{a}})+e_{h_1}(-\dot{\widehat{h_1}}).$$

Considering y as a virtual controller, then $y={\nu }_1(x)$

$$\dot{\widehat{V_1}}=x\left[x\left(1-{\displaystyle \frac{x}{K}}\right)-{\displaystyle \frac{ax{\gamma }_1}{1+x}}-h_{1}x+u_1\right]+e_a(-\dot{\widehat{a}})+e_{h_1}(\dot{\widehat{h_1}}).$$

(22)

Choosing ${\nu }_1=0$ and using the controller (15), Ref. (22) becomes

$$\dot{\widehat{V_1}}=x(-x-x{e}_{h_1})+e_a(\dot{\widehat{a}})+e_{h_1}(\dot{\widehat{h_1}}).$$

(23)

The updated law by the unknown parameter

$$\dot{\widehat{a}}=e_{a}and\dot{\widehat{h_1}}=-x^2+e_{h_1}.$$

(24)

Substituting (24) in (22), we get $\dot{V_1}=-x^2-e_a^2-e_{h_1}^2$ is the negative definite function, where

$${\nu }_1=0\Rightarrow {\eta }_1=y.$$

(25)

Differentiating (25), we get

$${\dot{\eta }}_1=\dot{y}={\displaystyle \frac{axy}{1+x}}-byz-h_{2}y+u_2.$$

(26)

Now, Equation (13) along with Equation (26), we get that

$$\dot{x}=-x-{\displaystyle \frac{ax{\eta }_1}{1+x}}-e_{h_1}x.$$

(27)

Considering the Lyapunov function of $(x,{\eta }_1),$

$$V_2=V_1+{\displaystyle \frac{1}{2}}{\eta }_1^2+{\displaystyle \frac{1}{2}}{e}_b^2+{\displaystyle \frac{1}{2}}{e}_{h_2}^2.$$

(28)

Differentiating (28) with respect to $t,$ we get

$$\begin{array}{cc}\hfill {\dot{V}}_2=& x\left(-x-{\displaystyle \frac{ax{\eta }_1}{1+x}}-e_{h_1}x\right)+{\eta }_1\left({\displaystyle \frac{ax{\eta }_1}{1+x}}-b{\eta }_{1}z-h_2{\eta }_1+u_2\right)\hfill \\ & e_a(-\dot{\widehat{a}})+e_{h_1}(-\dot{\widehat{h_1}})+e_b(-\dot{\widehat{b}})+e_{h_2}(-\dot{\widehat{h_2}}).\hfill \end{array}$$

(29)

Again, considering a new virtual controller $z={\nu }_2(x,y)$ where ${\nu }_2=0,$ and using this in (24), we have

$$\dot{V_2}=-x^2-e_{h_1}^2-e_a^2+e_{h_2}(-\dot{{\widehat{h}}_2})+{\eta }_1\left({\displaystyle \frac{-a{x}^2}{1+x}}+{\displaystyle \frac{ax{\eta }_1}{1+x}}-h_2{\eta }_1+u_2\right).$$

(30)

Now, choosing the controller (16) along with (30), we get

$${\dot{V}}_2=-x^2-e_a^2-e_{h_1}^2+e_{h_2}(\dot{\widehat{h_2}})+e_b(-\dot{\widehat{b}})-{\eta }_1^2.$$

(31)

The updated law for the unknown parameter $\dot{\widehat{b}}$ and $\dot{\widehat{h_2}}$ is

$$\dot{\widehat{b}}=e_{b}and\dot{\widehat{h_2}}=-{\eta }_1^2+e_{h_2}.$$

(32)

Substituting (32) in (31), we get

$${\dot{V}}_2=-x^2-e_b^2-e_a^2-e_{h_1}^2-e_{h_2}^2-{\eta }_1^2,$$

(33)

which is again a negative definite function where

$${\nu }_2=0\Rightarrow {\eta }_2=z.$$

(34)

Differentiating (34) with respect to $t,$ we have ${\dot{\eta }}_2=\dot{z}.$

Now,

$${\dot{\eta }}_1={\displaystyle \frac{ax{\eta }_1}{1+x}}-b{\eta }_{1}z-h_2{\eta }_1+u_2,$$

(35)

where the controller (17) along with (35) gives

$${\dot{\eta }}_1=-b{\eta }_1{\eta }_2-e_{h_2}{\eta }_1+{\displaystyle \frac{a{x}^2}{1+x}},$$

(36)

$${\dot{\eta }}_2=fyz-dz+u_3.$$

(37)

Now, considering the Lyapunov function of $(x,{\eta }_1,{\eta }_2),$

$$V_3=V_2+{\displaystyle \frac{1}{2}}{\eta }_2^2+{\displaystyle \frac{1}{2}}{e}_d^2+{\displaystyle \frac{1}{2}}{e}_f^2.$$

(38)

Differentiating (38) with respect to t gives

$$\begin{array}{cc}\hfill {\dot{V}}_3=& -x^2-e_{h_1}(x^2)-e_{h_2}{\eta }_1^2-{\eta }_1^2-{\eta }_1^2+{\eta }_2(-b{\eta }_1^2+f{\eta }_1{\eta }_2-d{\eta }_2+u_3)\hfill \\ & +e_a(-\dot{\widehat{a}})+e_b(-\dot{\widehat{b}})+e_f(-\dot{\widehat{f}})+e_{h_1}(-\dot{\widehat{h_1}})+e_{h_2}(\dot{\widehat{h_2}})+e_d(\dot{\widehat{d}}).\hfill \end{array}$$

(39)

The unknown parameters $\dot{\widehat{a}},\dot{\widehat{b}},\dot{\widehat{f}},\dot{\widehat{h_2}}$ are updated by

$$\dot{\widehat{a}}=e_a,\dot{\widehat{b}}=e_b,\dot{\widehat{f}}=e_f,\dot{\widehat{h_2}}=-{\eta }_1^2+e_{h_2},(\dot{\widehat{d}})=+d.$$

(40)

Substituting the updated parameters along with choosing the controller (17) and by updating law (40) in Equation (39), we get ${\dot{V}}_3=-x^2-{\eta }_1^2-{\eta }_2^2-e_a^2-e_b^2-e_f^2-e_{h_1}^2-e_{h_2}^2-e_d^2$ is the negative definite function. Thus, by the Lyapunov stability theorem, systems (10)–(12) is globally asymptotically stable with adaptive back-stepping controllers. □

8. Stochastic Analysis

All usual occurrences explicitly in the ecosystem are continuously under random fluctuations of the environment. The stochastic examination of any ecosystem gives an enhanced vision on the dynamic forces of the populace by means of population variances. In a stochastic model, the model parameters oscillate about their average values [24,25,26,27]. Therefore, the steady state which we anticipated as permanent will now oscillate around the mean state. The method to measure the mean-square fluctuations of population is proposed by [24] and it was applied by [28] nicely. Furthermore, many researchers like Samanta [29], Maiti, Jana and Samanta [30] have investigated critically the stochastic analysis to interpret local as well as global stability using mean-square fluctuations on population variances.

Now, this segment is meant for the extension of the deterministic model (1), which is formed by adding a noisy term. There are several ways in which environmental noise may be incorporated in the model system (1). External noise may arise from random fluctuations of a finite number of parameters around some known mean values of the population densities around some fixed values. Since the aquatic ecosystem always has unsystematic fluctuations of the environment, it is difficult to define the usual phenomenon as a deterministic ideal. The stochastic investigation enables us to get extra intuition about the continuous changing aspects of any ecological unit. The deterministic model (1) with the effect of random noise of the environment results in a stochastic system (41)–(43) given in the following discussion:

$${\displaystyle \frac{dx}{dt}}=x\left(1-{\displaystyle \frac{x}{K}}\right)-{\displaystyle \frac{axy}{1+x}}-h_{1}x+{\alpha }_1{\xi }_1(t),$$

(41)

$${\displaystyle \frac{dy}{dt}}={\displaystyle \frac{axy}{1+x}}-byz-h_{2}y+{\alpha }_2{\xi }_2(t),$$

(42)

$${\displaystyle \frac{dz}{dt}}=fyz-dz+{\alpha }_3{\xi }_3(t),$$

(43)

where ${\alpha }_1,{\alpha }_2,{\alpha }_3$ are the real constants and ${\xi }_i(t)=[{\xi }_1(t),{\xi }_2(t),{\xi }_3(t)]$ is a three-dimensional Gaussian white noise process. $E({\xi }_i(t))=0,$ where $i=1,2,3$; $E[{\xi }_i(t){\xi }_j(t)]={\delta }_{ij}\delta (t-t^{\prime }),$ where $i=j=1,2,3$; ${\delta }_{ij}$ is the Kronecker delta function; $\delta$ is the Dirac delta function.

Let

$$x(t)=u_1(t)+S^{*};y(t)=u_2(t)+P^{*};z(t)=u_3(t)+T^{*},$$

(44)

$${\displaystyle \frac{dx}{dt}}={\displaystyle \frac{d{u}_1(t)}{dt}};{\displaystyle \frac{dy}{dt}}={\displaystyle \frac{d{u}_2(t)}{dt}};{\displaystyle \frac{dz}{dt}}={\displaystyle \frac{d{u}_3(t)}{dt}}.$$

(45)

The linear parts of (41), (42) and (43) are (using (44) and (45))

$$u_1^{\prime }(t)=u_1(t)S^{*}-a{u}_2(t)S^{*}+{\alpha }_1{\xi }_1(t),$$

(46)

$$u_2^{\prime }(t)=a{u}_1(t)P^{*}-b{u}_3(t)P^{*}+{\alpha }_2{\xi }_2(t),$$

(47)

$$u_3^{\prime }(t)=f{u}_2(t)T^{*}+{\alpha }_3{\xi }_3(t).$$

(48)

Taking the Fourier transform on both sides of (46), (47) and (48), we get

$$i\omega {\tilde{u}}_1(\omega )=-S^{*}{\tilde{u}}_1(\omega )-a{S}^{*}{\tilde{u}}_2(\omega )+{\alpha }_1{\xi }_1(\omega ),$$

(49)

$$i\omega {\tilde{u}}_2(\omega )=a{P}^{*}{\tilde{u}}_1(\omega )-b{P}^{*}{\tilde{u}}_3(\omega )+{\alpha }_2{\xi }_2(\omega ),$$

(50)

$$i\omega {\tilde{u}}_3(\omega )=f{T}^{*}{\tilde{u}}_2(\omega )+{\alpha }_3{\xi }_3(\omega ).$$

(51)

The matrix form of (49)–(51) is

$$M(\omega )\tilde{u}(\omega )=\tilde{\xi }(\omega ),$$

(52)

where

$$M(\omega )=\left(\begin{array}{ccc}i\omega +S^{*}& a{S}^{*}& 0\\ -a{P}^{*}& i\omega & b{P}^{*}\\ 0& -f{T}^{*}& i\omega \end{array}\right);\tilde{u}(\omega )=\left[\begin{array}{c}{\tilde{u}}_1(\omega )\\ {\tilde{u}}_2(\omega )\\ {\tilde{u}}_3(\omega )\end{array}\right];\tilde{\xi }(\omega )=\left[\begin{array}{c}{\eta }_1{\tilde{\xi }}_1(\omega )\\ {\eta }_2{\tilde{\xi }}_2(\omega )\\ {\eta }_3{\tilde{\xi }}_3(\omega )\end{array}\right].$$

Equation (52) can also be written into

$$\tilde{u}(\omega )={[M(\omega )]}^{-1}\tilde{\xi }(\omega ),$$

(53)

where

$${[M(\omega )]}^{-1}={\displaystyle \frac{1}{R(\omega )+iI(\omega )}}\left(\begin{array}{ccc}{A}_1& D_1& G_1\\ B_1& E_1& H_1\\ C_1& F_1& I_1\end{array}\right)$$

(54)

and

$$A_1=-{\omega }^2+f{b}_3{T}^{*}{P}^{*};C_1=af{P}^{*}{T}^{*};D_1=i\omega {\alpha }_1{S}^{*};$$

$$E_1=-{\omega }^2+i\omega S^{*};F_1=i\omega f{T}^{*}+f{T}^{*}{S}^{*};G_1={\alpha }_1{b}_3{S}^{*}{P}^{*};$$

$$H_1=-i\omega b{P}^{*}-b{S}^{*}{P}^{*};I_1=-{\omega }^2+i\omega S^{*}+a{\alpha }_1{P}^{*}{S}^{*}.$$

Here,

$$|A_1{|}^2=X_1^2+Y_1^2;|B_1{|}^1=X_2^2+Y_2^2;{|C_1|}^2=X_3^2+Y_3^2;$$

$$|D_1{|}^2=X_4^2+Y_4^2;|E_1{|}^2=X_5^2+Y_5^2;{|F_1|}^2=X_6^2+Y_6^2;$$

$$|G_1{|}^2=X_7^2+Y_7^2;|H_1{|}^2=X_8^2+Y_8^2;{|I_1|}^2=X_9^2+Y_9^2,$$

where

$$\begin{array}{c}\hfill X_1=-{\omega }^2+fb{T}^{*}{P}^{*};Y_1=0;X_2=0;Y_2=(a\omega P^{*});X_3=af{P}^{*}{T}^{*};\\ \hfill Y_3=0;X_4=0;Y_4=\omega {\alpha }_1{S}^{*};X_5=\omega S^{*};X_6=\omega f{T}^{*};\\ \hfill Y_6=f{T}^{*}{S}^{*};X_7=0;Y_7={\alpha }_{1}b{S}^{*}{P}^{*};X_8=b{S}^{*}{P}^{*};Y_8=\omega b{P}^{*};\\ \hfill X_9=-{\omega }^2+a{\alpha }_1{S}^{*}{P}^{*};Y_9=\omega S^{*}.\end{array}$$

(55)

${|M(\omega )|}^2={[R(\omega )]}^2+{[I(\omega )]}^2$, where $R(\omega )=b_{3}f{T}^{*}{P}^{*}{S}^{*}-S^{*}{\omega }^2$ and $I(\omega )={\omega }^2+\omega b_{3}f{T}^{*}{P}^{*}+{\alpha }_{1}a\omega S^{*}{P}^{*}.$

If the function $Y(t)$ has a zero mean value, then the fluctuation intensity (variance) of its components in the frequency interval $[\omega ,\omega +d\omega ]$ is $S_Y(\omega )d\omega .$ $S_Y(\omega )$ is the spectral density of Y and is defined as

$$S_Y(\omega )=\underset{\tilde{T}\to \infty }{lim}{\displaystyle \frac{{\left|\tilde{Y}(\omega )\right|}^2}{\tilde{T}}}.$$

(56)

If Y has a zero mean value, the inverse transform of $S_Y(\omega )$ is the auto covariance function

$$C_Y(\tau )={\displaystyle \frac{1}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{S}_Y(\omega )e^{i\omega \tau }d\omega .$$

(57)

The corresponding variance of fluctuations in $Y(t)$ is given by

$${\sigma }_Y^2=C_Y(0)={\displaystyle \frac{1}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{S}_Y(\omega )d\omega$$

(58)

and the auto-correlation function is the normalized auto-covariance

$$P_Y(\tau )={\displaystyle \frac{C_Y(\tau )}{C_Y(0)}}.$$

(59)

For a Gaussian white noise process, it is

$$\begin{array}{cc}\hfill S_{{\xi }_i{\xi }_j}(\omega )=& \underset{\tilde{T}\to +\infty }{lim}{\displaystyle \frac{E\left[\widetilde{{\xi }_i}(\omega \widetilde{{\xi }_j}(\omega ))\right]}{\tilde{T}}}\hfill \\ \hfill =& \underset{\widehat{T}\to +\infty }{lim}{\displaystyle \frac{1}{\widehat{T}}}\underset{\frac{-\tilde{T}}{2}}{\overset{\frac{\tilde{T}}{2}}{\int }}\underset{\frac{-\tilde{T}}{2}}{\overset{\frac{\tilde{T}}{2}}{\int }}\left[\widetilde{{\xi }_i}(t)\widetilde{{\xi }_j}(t^{\prime })\right]e^{-i\omega (t-t^{\prime })}dtd{t}^{\prime }\hfill \\ \hfill =& {\delta }_{ij}.\hfill \end{array}$$

(60)

From (54), we have

$$\widetilde{u_i}(\omega )=\sum _{j=1}^3{K}_{ij}(\omega )\widetilde{{\xi }_j}(\omega );i=1,2,3.$$

(61)

From (59), we have

$$S_{ui}(\omega )=\sum _{j=1}^3{\eta }_j{|K_{ij}(\omega )|}^2;i=1,2,3,$$

(62)

where

$$K_{ij}(\omega )={[M(\omega )]}^{-1}.$$

Hence, by (61) and (62), the intensities of fluctuations in the variable $u_i$ $(i=1,2,3)$ are given by

$${\sigma }_{u_i}^2={\displaystyle \frac{1}{2\pi }}\sum _{j=1}^3\underset{-\infty }{\overset{\infty }{\int }}{\eta }_j{|K_{ij}(\omega )|}^{2}d\omega ;i=1,2,3$$

(63)

and by (54), we obtain

$${\sigma }_{{\mu }_1}^2={\displaystyle \frac{1}{2\pi }}\left\{\underset{-\infty }{\overset{\infty }{\int }}{\eta }_1{\left|{\displaystyle \frac{Adj(1)}{|M(\omega )|}}\right|}^{2}d\omega +\underset{-\infty }{\overset{\infty }{\int }}{\eta }_2{\left|{\displaystyle \frac{Adj(2)}{|M(\omega )|}}\right|}^{2}d\omega +\underset{-\infty }{\overset{\infty }{\int }}{\eta }_3{\left|{\displaystyle \frac{Adj(3)}{|M(\omega )|}}\right|}^{2}d\omega \right\},$$

(64)

$${\sigma }_{{\mu }_2}^2={\displaystyle \frac{1}{2\pi }}\left\{\underset{-\infty }{\overset{\infty }{\int }}{\eta }_1{\left|{\displaystyle \frac{Adj(4)}{|M(\omega )|}}\right|}^{2}d\omega +\underset{-\infty }{\overset{\infty }{\int }}{\eta }_2{\left|{\displaystyle \frac{Adj(5)}{|M(\omega )|}}\right|}^{2}d\omega +\underset{-\infty }{\overset{\infty }{\int }}{\eta }_3{\left|{\displaystyle \frac{Adj(6)}{|M(\omega )|}}\right|}^{2}d\omega \right\},$$

(65)

$${\sigma }_{{\mu }_1}^2={\displaystyle \frac{1}{2\pi }}\left\{\underset{-\infty }{\overset{\infty }{\int }}{\eta }_1{\left|{\displaystyle \frac{Adj(7)}{|M(\omega )|}}\right|}^{2}d\omega +\underset{-\infty }{\overset{\infty }{\int }}{\eta }_2{\left|{\displaystyle \frac{Adj(8)}{|M(\omega )|}}\right|}^{2}d\omega +\underset{-\infty }{\overset{\infty }{\int }}{\eta }_3{\left|{\displaystyle \frac{Adj(9)}{|M(\omega )|}}\right|}^{2}d\omega \right\}.$$

(66)

From (55), (64), (65) and (66),

$${\sigma }_{u_1}^2={\displaystyle \frac{1}{2\pi }}\left\{\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{1}{R^2(\omega )+I^2(\omega )}}\left[{\alpha }_1(X_1^2+Y_1^2)+{\alpha }_2(X_2^2+Y_2^2)+{\alpha }_3(X_3^2+Y_3^2)\right]d\omega \right\},$$

(67)

$${\sigma }_{u_2}^2={\displaystyle \frac{1}{2\pi }}\left\{\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{1}{R^2(\omega )+I^2(\omega )}}\left[{\alpha }_1(X_4^2+Y_4^2)+{\alpha }_2(X_5^2+Y_5^2)+{\alpha }_3(X_6^2+Y_6^2)\right]d\omega \right\},$$

(68)

$${\sigma }_{u_3}^2={\displaystyle \frac{1}{2\pi }}\left\{\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{1}{R^2(\omega )+I^2(\omega )}}\left[{\alpha }_1(X_7^2+Y_7^2)+{\alpha }_2(X_8^2+Y_8^2)+{\alpha }_3(X_9^2+Y_9^2)\right]d\omega \right\},$$

(69)

where

$$|M(\omega )|=R(\omega )+iI(\omega ).$$

If we are interested in the dynamics of stochastic process (41)–(69) with either ${\alpha }_1=0$ or ${\alpha }_2=0$ or ${\alpha }_3=0,$ the population variances are

if ${\alpha }_1=0,{\alpha }_2=0,$ then ${\sigma }_{u1}^2={\displaystyle \frac{{\alpha }_3}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{(X_3^2+Y_3^2)}{R^2(\omega )+I^2(\omega )}}d\omega ;{\sigma }_{u2}^2={\displaystyle \frac{{\alpha }_3}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{(X_6^2+Y_6^2)}{R^2(\omega )+I^2(\omega )}}d\omega ;$

${\sigma }_{u3}^2={\displaystyle \frac{{\alpha }_3}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{(X_9^2+Y_9^2)}{R^2(\omega )+I^2(\omega )}}d\omega .$

If ${\alpha }_2=0,{\alpha }_3=0,$ then ${\sigma }_{u1}^2={\displaystyle \frac{{\alpha }_1}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{(X_1^2+Y_1^2)}{R^2(\omega )+I^2(\omega )}}d\omega ;{\sigma }_{u2}^2={\displaystyle \frac{{\alpha }_1}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{(X_4^2+Y_4^2)}{R^2(\omega )+I^2(\omega )}}d\omega ;$

${\sigma }_{u3}^2={\displaystyle \frac{{\alpha }_1}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{(X_7^2+Y_7^2)}{R^2(\omega )+I^2(\omega )}}d\omega .$

If ${\alpha }_3=0,{\alpha }_1=0,$ then ${\sigma }_{u1}^2={\displaystyle \frac{{\alpha }_2}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{(X_2^2+Y_2^2)}{R^2(\omega )+I^2(\omega )}}d\omega ;{\sigma }_{u2}^2={\displaystyle \frac{{\alpha }_2}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{(X_5^2+Y_5^2)}{R^2(\omega )+I^2(\omega )}}d\omega ;$

${\sigma }_{u3}^2={\displaystyle \frac{{\alpha }_2}{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \frac{(X_8^2+Y_8^2)}{R^2(\omega )+I^2(\omega )}}d\omega .$

Equations (67)–(69) give three variations of the inhabitants. The integrations over the real line can be estimated, which gives the variations of the inhabitants.

9. Mathematical Model with Delay

In this section, we establish some conditions for oscillations of all positive solutions of the delay system

$${\displaystyle \frac{dx}{dt}}\left(t\right)=x\left(t\right)\left(1-h_1-{\displaystyle \frac{x\left(t-\tau \right)}{K}}-{\displaystyle \frac{ay\left(t-\tau \right)}{1+x\left(t-\tau \right)}}\right),$$

(70)

$${\displaystyle \frac{dy}{dt}}\left(t\right)=y\left(t\right)\left(-h_2+{\displaystyle \frac{ax\left(t-\tau \right)}{1+x\left(t-\tau \right)}}-bz\left(t-\tau \right)\right),$$

(71)

$${\displaystyle \frac{dz}{dt}}\left(t\right)=z\left(t\right)\left(-d+fy\left(t-\tau \right)\right).$$

(72)

Here, the parameter $\tau \in {\mathbb{R}}^{+}$ is the delay. This proposed system is concerned not only with the present number of predator and prey but also with the number of predator and prey in the past. If t is the present time, then (t-$\tau$) is the past.

According to Krisztin [31], a solution of (70)–(72), $(x(t),y(t),z(t))$ is called oscillatory if every component has arbitrary large zeros; otherwise, $(x(t),y(t),z(t))$ is said to be a non-oscillatory solution. Whenever all solutions of (70)–(72) are oscillatory, we will say that (70)–(72) is an oscillatory system.

In [32], Kubiaczyk and Saker studied the oscillatory behavior of the delay differential equation

$$x^{\prime }(t)+px(t)-\frac{qx(t)}{r+x^n(t-\tau )}=0,$$

where $p,q,r,\tau \in {\mathbb{R}}^{+}.$ Using similar methods to liberalize each equation of the delay system, we will establish conditions for oscillations of all positive solutions of the system.

Now, we will analyze the system of (70)–(72) with the following initial conditions:

$$x(\theta )\ge 0,y(\theta )\ge 0,z(\theta )\ge 0,\theta \in (-\tau ,0),(\tau \ne 0).$$

(73)

Using the same arguments that we got in Theorem 1, we can establish the following theorem:

Theorem 7.

All solutions of $(x(t),y(t),z(t))$ of systems (70)–(72) with the initial condition (73) are positive for all $t\ge 0.$

Easily, we can see that the equilibrium point remains the same when we have the delay system. However, it is important to know the oscillatory behavior of the solutions around the equilibrium points.

Theorem 8.

If there exist a $\lambda \in \mathbb{R}$ such that

$$\underset{\lambda \in \mathbb{R}}{min}\left({\lambda }^3-\alpha {\lambda }^2{e}^{-\lambda \tau }-\left(d\sigma +\beta \gamma \right)\lambda e^{-2\lambda \tau }+d\sigma \alpha e^{-3\lambda \tau }\right)>0,$$

(74)

where $\alpha =-\frac{x^{*}}{k}$, $\beta =-\frac{ad}{f\left(1+x^{*}\right)},$ $\gamma =\frac{a{x}^{*}}{1+x^{*}}$ and $\sigma =-d{z}^{*}$; then, all solutions of the system (70)–(72) oscillate around $E_3$.

Proof. 

Let us consider the system (70)–(72). Let

$$x(t)=x^{*}{e}^{{\varphi }_1\left(t\right)},$$

(75)

$$y(t)=y^{*}{e}^{{\varphi }_2\left(t\right)},$$

(76)

$$z(t)=z^{*}{e}^{{\varphi }_3\left(t\right)}.$$

(77)

Then, $\left(x(t),y(t),z(t)\right)$ oscillate around $\left(x^{*},y^{*},z^{*}\right)$ if $\left({\varphi }_1(t),{\varphi }_2(t),{\varphi }_3(t)\right)$ oscillate around $\left(0,0,0\right).$ From (70)–(72) and (75)–(77), we have

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{\varphi }_1}{dt}}\left(t\right)& =& 1-h_1-{\displaystyle \frac{x^{*}{e}^{{\varphi }_1\left(t-\tau \right)}}{K}}-{\displaystyle \frac{a{y}^{*}{e}^{{\varphi }_2\left(t-\tau \right)}}{1+x^{*}{e}^{{\varphi }_1\left(t-\tau \right)}}}\hfill \\ & =& 1-h_1-{\displaystyle \frac{x^{*}}{K}}-{\displaystyle \frac{a{y}^{*}}{1+x^{*}{e}^{{\varphi }_1\left(t-\tau \right)}}}-{\displaystyle \frac{x^{*}}{K}}\left(e^{{\varphi }_1\left(t-\tau \right)}-1\right)\hfill \\ & & -{\displaystyle \frac{a{y}^{*}}{1+x^{*}{e}^{{\varphi }_1\left(t-\tau \right)}}}\left(e^{{\varphi }_2\left(t-\tau \right)}-1\right),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{\varphi }_2}{dt}}\left(t\right)& =& -h_2+\frac{a{x}^{*}{e}^{{\varphi }_1\left(t-\tau \right)}}{1+x^{*}{e}^{{\varphi }_1\left(t-\tau \right)}}-b{z}^{*}{e}^{{\varphi }_3\left(t-\tau \right)}\hfill \\ & =& -h_2+\frac{a{x}^{*}}{1+x^{*}{e}^{{\varphi }_1\left(t-\tau \right)}}-b{z}^{*}+\frac{a{x}^{*}}{1+x^{*}{e}^{{\varphi }_1\left(t-\tau \right)}}\left(e^{{\varphi }_1\left(t-\tau \right)}-1\right)\hfill \\ & & -b{z}^{*}\left(e^{{\varphi }_3\left(t-\tau \right)}-1\right),\hfill \end{array}$$

$${\displaystyle \frac{d{\varphi }_3}{dt}}\left(t\right)=-d+f{y}^{*}+f{y}^{*}\left(e^{{\varphi }_2\left(t-\tau \right)}-1\right).$$

Let

$$f\left(u\right)=e^u-1.$$

Note that

$$uf\left(u\right)>0\mathrm{for}u>0,\mathrm{and}\underset{u\to 0}{lim}\frac{f(u)}{u}=1.$$

(78)

Moreover, since $\left(x^{*},y^{*},z^{*}\right)$ is the equilibrium point $E_3$, we have

$$1-h_1-{\displaystyle \frac{x^{*}}{K}}-{\displaystyle \frac{a{y}^{*}}{1+x^{*}}}=0,-h_2+\frac{a{x}^{*}}{1+x^{*}}-b{z}^{*}=0,\mathrm{and}-d+f{y}^{*}=0.$$

Thus, the linearized system associated with the system (75)–(77) is given by

$${\displaystyle \frac{d{m}_1}{dt}}\left(t\right)=-{\displaystyle \frac{x^{*}}{K}}{m}_1\left(t-\tau \right)-{\displaystyle \frac{a{y}^{*}}{1+x^{*}}}{m}_2\left(t-\tau \right),$$

(79)

$${\displaystyle \frac{d{m}_2}{dt}}\left(t\right)=\frac{a{x}^{*}}{1+x^{*}}{m}_1\left(t-\tau \right)-b{z}^{*}{m}_3\left(t-\tau \right),$$

(80)

$${\displaystyle \frac{d{m}_3}{dt}}\left(t\right)=f{y}^{*}{m}_2\left(t-\tau \right),$$

(81)

and every solution of (79)–(81) oscillates if and only if the characteristic equation has no real roots (see Theorem 5.1.1 in [21]), i.e.,

$$det\left[\lambda I-A\right]\ne 0$$

(82)

for all $\lambda \in \mathrm{R}.$ Equation (82) is equivalent to the equation

$${\lambda }^3+\alpha {\lambda }^2{e}^{-\lambda \tau }-\left(d\sigma +\beta \gamma \right)\lambda e^{-2\lambda \tau }+d\sigma \alpha e^{-3\lambda \tau }=0,$$

where $\alpha =\frac{x^{*}}{k}$, $\beta =-\frac{ad}{f\left(1+x^{*}\right)},$ $\gamma =\frac{a{x}^{*}}{1+x^{*}}$ and $\sigma =d{z}^{*}.$ In fact,

$$\underset{\lambda \to +\infty }{lim}\left({\lambda }^3+\alpha {\lambda }^2{e}^{-\lambda \tau }-\left(d\sigma +\beta \gamma \right)\lambda e^{-2\lambda \tau }+d\sigma \alpha e^{-3\lambda \tau }\right)=+\infty$$

and

$$\underset{\lambda \to -\infty }{lim}\left({\lambda }^3+\alpha {\lambda }^2{e}^{-\lambda \tau }-\left(d\sigma +\beta \gamma \right)\lambda e^{-2\lambda \tau }+d\sigma \alpha e^{-3\lambda \tau }\right)=+\infty ,\mathrm{since}d\sigma \alpha >0.$$

Then, by (74) and (78), systems (79)–(81) will start oscillating and then all the solutions of systems (75)–(77) will also oscillate. □

Example 1.

Let the parameters of systems (70)–(72) be $a=1.5,$ $d=0.2,$ $b=0.3,$ $f=0.65,$ $h_1=0.5,$ $h_2=0.3,$ $K=9$ and $\tau =2.$ In this case, condition (74) becomes

$$\underset{\lambda \in \mathbb{R}}{min}\left({\lambda }^3+0.3995837778{\lambda }^2{e}^{-\lambda \tau }+0.234339529\lambda e^{-2\lambda \tau }+0.046546037e^{-3\lambda \tau }\right)>0.045$$

and consequently all solutions of the system oscillate around the equilibrium point $\left(3.596254004,0.307692308,2.912157597\right).$

10. Numerical Simulations

Analytical studies can never complete without numerical verification of the results. Moreover, it may be noted that, as the parameters of the model are not based on the real world observation, the main features described by the simulations presented in this section should be considered from a qualitative rather than a quantitative point of view. We choose the parameters of system (1) as $a=1.5,$ $b=0.3,$ $f=0.65,$ $h_1=0.5,$ $h_2=0.3,$ $d=0.2,$ $K=8$ with the initial densities $x(0)=2,$ $y(0)=1.8,$ $z(0)=1$ and observe the dynamical behaviour of system (1). Figure 1a shows that the equilibrium point $E_1$ is locally asymptotically stable and the corresponding phase-portrait is shown in Figure 1b. Figure 2a shows that the equilibrium point $E_2$ is locally asymptotically stable and the corresponding phase-portrait is shown in Figure 2b. Figure 3a shows that the co-existence equilibrium point $E_3$ is locally asymptotically stable and the corresponding phase portrait is shown in Figure 3b. From Figure 4a–d, if all other parameters are fixed and varying transmission rate $a=1.5$ to $a=2$, we observe that oscillation settles down into a stable situation for all three of the species. Stability around this state implies extinction of the infected prey. This study interestingly suggests that the harvesting of both prey prevent limit cycle oscillations and the combined effect of both harvests also prevent disease propagation in the system. We also conclude that the inclusion of stochastic perturbation creates a significant change in the intensity of populations because change of responsive noise parameters causes chaotic dynamics with low, medium and high variances of oscillations (see Figure 5, Figure 6 and Figure 7).

11. Conclusions

In this paper, we have studied stability of a diseased model of susceptible, infected prey and predators around an interior steady state. The positivity of the solutions and boundedness of the system together with stability analysis of boundary equilibrium providing all the necessary information to establish persistence of the system. The deterministic situation theoretical epidemiologists are usually guided by an implicit assumption that most epidemic models (we observe in nature) correspond to stable equilibrium of models. In Theorem 3, we gave the condition for stable co-existence. In Theorem 4, we proved the global stability by using a Lyapunov function. We also worked out the condition for which all three species will persist. The controllability conditions and the conditions for global asymptotic stability have been obtained by using the adaptive back-stepping control approach by using a suitable Lyapunov function. We also studied the stochastic perturbation of model (1), which generates a significant change in the intensity of populations due to low, medium and high variances of oscillations.