# Analysis of Stochastic Predator-Prey Model with Disease in the Prey and Holling Type II Functional Response

Authors: C. Gokila; M. Sambath; K. Balachandran; Yong-Ki Ma
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

A stochastic predator-prey model with disease in the prey and Holling type II functional response is proposed and its dynamics is analyzed. We discuss the boundedness of the dynamical system and find all feasible equilibrium solutions. For the stochastic systems, we obtain the conditions for the existence of the global unique solution, boundedness, and uniform continuity. We derive the conditions for extinction and permanence of species. Moreover, we construct appropriate Lyapunov functions and discuss the asymptotic stability of equilibria. To illustrate our theoretical findings, we have performed numerical simulations and presented the results.


---

## Body

## 1. Introduction

Mathematical models are used to study the interrelationship among species and their environment. The study of disease transmission has turned out to be a valuable field of research after the fundamental work of Kermac and McKendric [1] on susceptible-infected framework. Hadeler and Freedman [2] first proposed a disease spread model within interacting populations. Initially, epidemics are created if there are some people susceptible to the infection and some infected people in the population. It is especially essential to view the ecosystem with the influence of epidemiological factors to control the disease in the species. From the ecological point of view, the spread of disease can not be disregarded because its effects are serious. So, various authors have paid attention to the study of transmissible disease in ecology, see for example [3–6] and the references therein. Mondal [7] has examined the disease model with two species and analyzed the dynamical properties of the fractional order system. Haque and Venturino [8] investigated the stability behavior of the deterministic Holling-Tanner predator-prey model. In this paper, we propose the predator-prey model and consider the Holling type II response for predation.In an ecological model, the interactions between two or more species and their dynamics are influenced by each other. So, the growth of one species depends on another and is described by the prey-predator system. Three primary kinds of interaction between the species are: predator-prey, mutualism, and competition. In all predator-prey interactions, Holling functions do not allow the growth of predators to very large extent even if the density of the prey is more. Specifically, Holling type II functional response is defined by a decelerating intake rate which follows from the assumption that the consumer is limited by its capacity to process food. In other words, Holling type II represents the fact that when prey density is small, the predator can take less time for handling prey and if the prey density increases, more prey are attacked so that the handling time also increases. In this article, we have used Holling type II response for both infected and susceptible prey interactions with the predator. This kind of functional response has been widely utilized as a part of biological systems, see few epidemic models [9–11] and chemostat model [12].During the past decades, a study of dynamical behavior of the population species with stochastic impacts has been growing steadily. The interesting situation occurs at the global stability of all feasible equilibria. Pitchaimani and Rajaji [13] constructed the stochastic Nowak-May model and investigated the asymptotic stability. In addition to stability, for every population model, the problem of permanence and boundedness property is also important. Solutions of the population model are called ultimately bounded if they satisfy the following condition: if we find the existence of bounded region in the solution space of our system such that each solution enters the bounded region in limited time and remains within the region forever. The permanence gives a guarantee that if initially the density of all species is positive, then after a specific time the density of each species will be present in some sizeable amount. Ghosh et al. [14] illustrated a seasonally perturbed stochastic model and analyzed the persistence for three species. In the literature, many results that study stability, boundedness, and persistence have been presented for some ecological models with stochastic effect [6, 10, 15–18].Nonetheless, parameters associated with the system are not fully constant, and they always change with time around some average values. These fluctuations occur because of sudden changes in the environment [19] or often created by human interference and natural events in the ecosystem by disturbing the environment. Environmental changes are described as natural disasters, human intervention, or animal or bird contact or infestation of invasive species. So these environmental changes can be outlined as noise. These changes are extreme and produce more effect on the population size in a particular time. Bringing environmental fluctuations into the predator-prey model is the correct way to deal with this situation. May [20] has uncovered the reality that the birth rates, death rates, carrying capacities, and other parameters which describe the remaining factors involved in the ecological process carry the randomness to a large or low extent due to ecological fluctuations. Accordingly, as time tends to be large, every equilibrium solution does not achieve a steady-state value accurately but it fluctuates continuously around the steady state. Recently, Liu et al. [21] developed and analyzed a population model with Holling II response and random effect. To study the model with fluctuations, several authors have introduced ecological fluctuations into every population model to accentuate the reality [11, 13, 16–18, 22–25]. The predator-prey model with two species and ratio dependence is discussed to examine its stability of equilibrium solutions in [26]. Ji et al. [9] introduced two types of functional response and stochastic perturbation into the system. Zhang et al. [27] found the critical value for the stochastic predator-prey system which can be used to determine the extinction and persistence in the mean of the predator population. Zhang and Meng [6] developed the nonautonomous SIRI epidemic model with random disturbance. The above researchers used various noises and different types of functional response depending on the population model. By the above motivation, we consider the predator-prey model with environmental changes in this article.The article is arranged as follows. In Section 2, we present few definitions, lemmas, and theorems which are utilized in further analysis. In Section 3, we discuss the detailed explanation about the formulation and the condition that solution of the deterministic model is bounded. For the stochastic system, we derive the existence of positive solution of the system and its uniqueness and also explore the conditions for stochastic boundedness in Section 4. In addition, we prove that the solution is uniformly continuous. In Section 5, stochastic permanence and extinction under certain parametric restriction are established. Using the corresponding Lyapunov function, we have examined the conditions on global asymptotic stability in Section 6. Next, we have obtained some figures to justify the results in Section 7. Finally, the conclusion based on our results is presented in Section 8.


## 2. Preliminaries

Here, we give certain notations, definitions, theorems, and lemmas which are used in the following analysis. For more details, see [28–31].Consider the stochastic model (SM) ofd-dimension of the form
(1)dZt=gZt,tdt+hZt,tdBt,t0≤t≤T<∞,with Zt0=Z0. The functions g:Rd×t0,T⟶Rd and h:Rd×t0,T⟶Rd×m are Borel measurable, B=Btt≥t0 is an Rm-valued Wiener process, and Z0 is an Rd-valued random variable.The differential operatorL corresponding to the SM (1) is defined as
(2)L=∂∂t+∑i=1dgiz,t∂∂zi+12∑i,j=1mhz,tgTz,tij∂2∂zi∂zj.Along with the existence and uniqueness assumptions, we make the assumption thatg and h satisfy gz∗,t=0 and hz∗,t=0 for an equilibrium solution z∗, for t≥t0.Definition 1.
The equilibrium solutionz∗ of the SM (1) is stochastically stable if it satisfies for every ε>0 and s≥t0,
(3)limz0⟶z∗ℙsups≤t∥Zs,z0t−z∗∥≥ε=0,where Zs,z0t represents the solution of (1) with Zs=z0 at time t≥s.Definition 2.
The equilibrium solutionz∗ of the SM (1) is said to be stochastically asymptotically stable if it satisfies the stochastic stability condition and
(4)limz0⟶z∗ℙlimt⟶∞Zs,z0t=z∗=1.Definition 3.
The equilibrium solutionz∗ of the SM (1) is said to be globally stochastically asymptotically stable if it satisfies the stochastic stability condition and for every z0 and every s,
(5)ℙlimt⟶∞Zs,z0t=z∗=1.Theorem 4 (see [28]).
Let the functionsg and h have continuous coefficients with respect to t and satisfy the existence and uniqueness properties.
(i)
Suppose that a positive definite functionV∈C2,1Uk×t0,∞ exists, where Uk=z∈Rd:∥z−z∗∥<k, for k>0, such that for allt≥t0,z∈Uk:LVz,t≤0, then, the solution z∗ of (1) is stochastically stable(ii)
Additionally ifV is decresing and a positive definite function V1 exists such that(6)Vz,t≤V1zforallt≥t0,z∈Uk:LVz,t≤0,then the equilibrium solution z∗ is stochastically asymptotically stable
(iii)
If the assumption (ii) holds for a radially unbounded functionV∈C2,1Rd×t0,∞ defined everywhere, then the equilibrium solution z∗ is globally stochastically asymptotically stableLemma 5 (see [29, 31].
Suppose that a stochastic processZt on t≥0 of n-dimension satisfies
(7)EZt−Zsα1≤ct−s1+α2,0≤s,t<∞,where α1, α2, and c are arbitrarily nonnegative constants and a continuous modification Z~t of Zt exists having the property that, for every υ∈0,α2/α1, there exists a random variable ψω>0 such that
(8)ℙω:sup0<t−s<ψω,0≤s,t<∞Z~t,ω−Zt,ωt−sυ≤21−2−υ=1,that is, each sample path of Z~t is locally but uniformly Hölder continuous with υ.Definition 6 (see [30]).
The solutionZt of model (1) is said to be stochastically ultimately bounded, if, for any ε∈0,1, there is a constant δ=δε>0, such that for any initial value Z0∈R+3, the solution Zt of (1) satisfies
(9)limsupt⟶∞ℙZt>δ<ε.Definition 7 (see [30]).
The solutionZt of (1) possesses stochastic permanent property, if there exists a pair of constants φ=φν>0 and χ=χν>0 for any ν∈0,1 such that the solution Zt of (1) for any initial value Z0∈R+3 satisfies the property
(10)liminft⟶∞ℙ∣Zt∣≥φ≥1−ν,liminft⟶∞P∣Zt∣≤χ≥1−ν.


## 3. Deterministic Model

In this section, we propose a predator-prey model with disease among the prey population. Chattopadhyay and Bairagi [32] framed the ecoepidemiological model with two species dividing into three compartments in the Salton sea and analyzed the stability of the positive equilibrium. Because of the disease, susceptible prey and infected prey are there as two groups in the prey population. The predator mostly eats infected prey because they are easy to catch. So these infected preys become more attractive to the predator. We have assumed that both the preys are subject to predation by the predator. In our article, we considered the population model as in [32] with the inclusion of the susceptible prey and predator interaction and functional response as Holling type II for interaction in the following form:
(11)dSdt=rS1−S+IK−λIS−αPSm+S,dIdt=−μI+λIS−βIPa+I,dPdt=−dP+βIPa+I+αPSm+S.Here,St, It, and Pt denote the population densities of susceptible prey, infected prey, and predator at any time t with S0=S0≥0, I0=I0≥0, and P0=P0≥0. r, K, and λ represent the growth rate of S, carrying capacity of susceptible prey, and disease transmission coefficient. α is the search rate of the predator towards susceptible prey and β is the search rate of predator towards infected prey, μ and d are the natural death rates of infected prey and predator. Parameters m and a are half saturation constants. System (11) can have at most five equilibrium solutions:
(i)
The trivial equilibrium solutionE0=0,0,0(ii)
The equilibrium solutionE1=K,0,0 lying on the boundary(iii)
The planar equilibrium solutionE2=S2,I2,0 on the S−I plane where S2=μ/λ and I2=rKλ−c/λr+Kλ(iv)
Another planar equilibrium solutionE3=S3,0,P3 on the S−P plane where S3=dm/α−d and P3=rKmα−dm2−dKm/Kd−α2(v)
The positive equilibrium solutionE4=S∗,I∗,P∗ which is obtained as follows:LetS∗ be a nonnegative root of the following equation
(12)gS∗=S∗2+CS∗+D=0,where C=1/r−d+α+βadr+dKr−dmr−arα−Krα−Krβ+mrβ+adKλ and D=1/r−d+α+βadmr+dKmr−aKαμ−Kmrβ+adKmλ.The roots of the above quadratic equation are(13)S∗=−C±C2−4D2.When any one of the following cases is satisfied, the equilibrium solutionS∗ can have one or two positive values.
(i)
C<0 and D<0(ii)
C<0,D>0 and C2−4D>0(iii)
C>0,D<0where I∗=adm+S∗−αS∗/β−dm+S∗+αS∗,P∗=λS∗−μa+I∗/β.The following relations must hold for the positiveness ofI∗ and P∗:
(14)αS∗<dm+S∗<βm+S∗+αS∗andλS∗>μ.The positive equilibrium solution plays a major role in changing the dynamical behavior. It is the only solution where all the species exist. All other equilibria are the subcases of the coexisting equilibrium solution. Therefore, it is essential to analyze the dynamical properties of positive equilibrium and also it gives the behavior of each species exactly.Now, we provide certain conditions to bound the solutions of the system through the boundedness of the model equation (11).Theorem 8.
All the solutions of system (11) in R+3 with positive initial conditions are uniformly bounded.Proof.
To get the boundedness of solutions of given system (11), we consider the function
(15)U=S+I+P.
Differentiate the above equation with respect to timet to obtain
(16)dUdt=dSdt+dIdt+dPdt,=rS1−S+IK−λIS−αSPm+S−μI+λIS−βIPa+I−dP+βIPa+I+αPSm+S,=rS1−S+IK−μI−dP.
For each0≤η≤minμ,d, the following inequality holds
(17)dUdt+ηU≤r1−SK+ηS+η−μI+η−dP,≤K4rr+η2.
The maximum value of the quadratic functionax2+bx+c is c−b2/4a when a<0. In this way, we get the maxR+r1−S/K+ηS as K/4rr+η2 (refer [33]). Assume that L=K/4rr+η2>0, this implies
(18)dUdt+ηU<L.
By the theory of differential inequalities, we get(19)0≤US,I,P<Lη1−e−ηt+US0,I0,P0e−ηt,and letting t tend to infinity, the above solution is of the form
(20)0≤US,I,P<Lη.
From the above discussion, we conclude that the solution space of system (11) lies within
(21)D=S,I,P∈R+3;U=Lη+ε,foreveryε>0.
Hence, the theorem is proved.


## 4. Stochastic Model

In the natural world, each population in an ecosystem is greatly affected by environmental noises which play a major role in population dynamics. By considering the effect of random environment fluctuations, we have included environmental noise in every equation of our deterministic system (11). In our system, the randomness in the environment will directly affect themselves as fluctuations in the growth rate of the susceptible prey, death rate of the infected prey population, and predator population like
(22)r⟶r+σ1Ḃ1t,−μ⟶−μ+σ2Ḃ2t,−d⟶−d+σ3Ḃ3t,where Bit,i=1,2,3 are independent Brownian motions and σi2i=1,2,3 denote the intensities of the environmental fluctuations and σii=1,2,3 represent the standard deviation. With this fact, we have framed the stochastic system by using the Itô equations as follows:
(23)dS=Sr1−S+IK−λI−αPm+Sdt+σ1SdB1t,dI=I−μ+λS−βPa+Idt+σ2IdB2t,dP=P−d+βIa+I+αSm+Sdt+σ3PdB3t.During the past several years, no work has been reported on the above stochastic model (23). Our aim is to find the dynamics of the stochastic system (23) and show how each population varies with respect to environmental fluctuations.Now, we discuss some important properties like positiveness, boundedness, and continuity of solution of the stochastic model (23).Theorem 9.
ForS0,I0,P0∈IntR+3, system (23) has a unique positive local solution St,It,Pt for t∈0,τe almost surely, where τe is the explosion time.Proof.
Consider the transformation of variables(24)x=logS,y=logI,z=logP.
Using the Itô formula,(25)LV=Vtt,x+Vxt,xft,x+12tracegTt,xVxxt,xgt,x,we get
(26)Lx=r1−ex+eyK−λey−αezm+ex−12σ12,dx=Lxdt+σ1dB1t.
Similarly, we obtain, from system (23),
(27)dx=r1−ex+eyK−λey−αezm+ex−12σ12dt+σ1dB1t,(28)dy=λex−βeza+ey−μ−12σ22dt+σ2dB2t,(29)dz=βeya+ey+αexm+ex−d−12σ32dt+σ3dB3t,with x0=logS0,y0=logI0, and z0=logP0. Now, the functions corresponding to system (28) have initial growth and they satisfy the local Lipchitz property. Hence, a unique local solution xt,yt,zt exists and it is defined in 0,τe. Consequently, there exists a unique positive local solution of (23) as St=ext, It=eyt, and Pt=ezt.Theorem 10.
System (23) has a unique solution St,It,Pt for t∈0,τe and for any initial condition S0,I0,P0∈IntR+3 and the solution remains in IntR+3 with probability one. Therefore, for all t≥0, St,It,Pt∈IntR+3 almost surely.Proof.
To show that the global solution exists, it is enough to prove thatτ=∞ almost surely. Assume that κ0 is a large nonnegative integer such that the closed ball Bκ0∈R+3 contains S0,I0,P0. We choose for any κ≥κ0 and define the stop-time as
(30)τκ=inft∈0,τe:S∉1κ,κorI∉1κ,κorP∉1κ,κ.
Here,inf∅=∞ (∅ is the empty set). Therefore, τκ is increasing as κ⟶∞.
Letτ∞=limκ⟶∞τκ; then, τ∞≤τe almost surely. If τ∞=∞ almost surely is true, then τe=∞ almost surely. If this statement fails, that is, if τ∞≠∞, then the two constants T>0 and ε∈0,1 exist with
(31)ℙτ∞≤T>ε.
Thus, by denotingΩκ=τκ≤T, then κ1≥κ0 is an integer such that, for all κ≥κ1,
(32)ℙτκ≤T≥ε.
DefineV:IntR+3⟶IntR+ by
(33)VS,I,P=S−1−lnS+I−1−lnI+P−1−lnP,where the function VS,I,P>0 for all S,I,P∈IntR+3. Using Itô’s formula, we get
(34)LV=1−1SrS1−S+IK−λSI−αPSm+S+12σ12+1−1I−μI+λSI−βIPa+I+12σ22+1−1P−dP+βIPa+I+αPSm+S+12σ32,=S−1r1−S+IK−λI−αPm+S+I−1−μ+λS−βPa+I+P−1−d+βIa+I+αSm+S+12σ12+12σ22+12σ32,=S−1r−rKS+I−λI+I−1λS−μ−dP−1+αSP−1−αPS−1m+S+βIP−1−βPI−1a+I+σ12+σ22+σ322.
Taking the differential ofVS,I,P, one gets
(35)dVS,I,P=fS,I,Pdt+gS,I,PdBt,where
(36)gS,I,P=σ1S−1+σ2I−1+σ3P−1,S,I,P=σ12+σ22+σ322−S−1rKS+I+λI−r−I−1μ−λS−dP−1+αSP−1−αPS−1m+S+βIP−1−βPI−1a+I.
From [29], it is easy to show thatfS,I,P is bounded above, say by N, in R+3 (that is to say ∣fS,I,P∣≤N, for all S,I,P∈R+3). From equation (35), we have
(37)∫0τm∧TdVS,I,P≤∫0τm∧TNdt+∫0τm∧TgS,I,PdBt,where τκ∧T=minτκ,T. Taking the expectations of the above inequality, one gets
(38)EVSτκ∧T,Iτκ∧T,Pτκ∧T≤VS0,I0,P0+NT.
Note that no less than one ofSτκ,ω, Iτκ,ω, and Pτκ,ω belongs to the set κ,1/κ, for every ω∈Ωκ; therefore, we get
(39)VSτκ∧T,Iτκ∧T,Pτκ∧T≥κ−1−lnκ∧1κ−1−ln1κ.
Hence, from (32),
(40)EVSτκ∧T,Iτκ∧T,Pτκ∧T≥EIΩκωVSτκ∧T,Iτκ∧T,Pτκ∧T,≥εκ−1−lnκ∧1κ−1−ln1κ,in which IΩκω denotes the indicator function of Ωκ. It follows from (39) that
(41)VS0,I0,P0+NT≥εκ−1−lnκ∧1κ−1−ln1κ.
κ⟶∞ leads to a contradiction: ∞>VS0,I0,P0+NT=∞. Therefore, τ∞=∞ almost surely. Thus, τ=∞ then St,It,Pt∈R+3 almost surely.With the existence of solution, next, we analyze how the solution changes inR+3.Theorem 11.
For any initial valueW0=S0,I0,P0∈R+3, the solutions of system (23) are stochastically ultimately bounded.Proof.
By Theorem10, the solution Wt remains in R+3 for all t≥0. Consider the function V1t,S=etSθ for θ>0. Using the Itô formula, we compute
(42)LV1=etSθ1+r1−S+IK−λI−αPm+Sθ+θ−1θ2σ12,≤et1+rθ+θ−1θ2σ12Sθ−rθkSθ+1≤M1θet.
By considering the integral and expectation on two sides of the above equation, we getetESθt−ES0θ≤M1θet. So, we have limsupt⟶∞ESθt≤M1θ<+∞.
Define the functionV2t,I=etIθ for θ>0; using the Itô formula, we get
(43)LV2=etIθ1+−μ+λS−βPa+Iθ+θ−1θ2σ22,≤etIθ1+λθ−μθ+θ−1θ2σ22≤M2θet.
Then,etEIθt−EI0θ≤M2θet. So, we have limsupt⟶∞EIθt≤M2θ<+∞. Similarly, defining the function V3t,P=etPθ for θ>0 and applying the Itô formula, we get
(44)LV3=etPθ1+−d+βIa+I+αSm+Sθ+θ−1θ2σ32,≤etPθ1+βaθ+αmθ−dθ+θ−1θ2σ32≤M3θet.
Then,etEPθt−EP0θ≤M3θet. So, we have limsupt⟶∞EPθt≤M3θ<+∞.
ForWt=St,It,Pt∈R+3, we may get
(45)Wtθ≤3maxS3t,I3t,P3tθ/3≤3θ/3Sθt+Iθt+Pθt.
Consequently,(46)limsupt⟶∞EWtp≤M4θ<+∞,where M4θ=3θ/3M1θ+M2θ+M3θ. Applying Chebyshev inequality, we get that all solutions are stochastically bounded.Using fundamental properties and suitable Lyapunov functions, we continue to show that the positive solutionWt=St,It,Pt is uniformly Hölder continuous.Theorem 12.
Every sample path ofSt,It,Pt is uniformly continuous, where St,It,Pt is a solution of system (23) on t≥0 with S0,I0,P0∈R+3.Proof.
The modified form of the first equation of system (23) is
(47)St=S0+∫0tSur1−S+IK−λI−αPm+Sdu+∫0tσ1SudB1u.
Assume thatf1u=Sur1−Su+Iu/K−λIu−αPu/m+Su andf2u=σ1Su.
From Theorem (28), we deduce that
(48)Ef1tθ=ESr1−S+IK−λI−αPm+Sθ,=ESθr−rS+IK−λI−αPm+Sθ,≤12ES2θ+12Er+rS+r+λI+αP2θ,≤12ES2θ+44θ−3/2r2θ+r2θES2θ+r+λ2θEI2θ+α2θEP2θ,≤12M12θ+44θ−3/2r2θ+r2θM12θ+r+λ2θM22θ+α2θM32θ,=F1θ,Ef2tθ=Eσ1Suθ=σ1θESθ≤σ1θM1θ≤F2θ.
For stochastic integrals, we observe the moment inequality and apply for0≤t1≤t2 and θ>2, to get
(49)E∫t1t2f2udB1uθ≤θ−1θ2θ/2t2−t1θ−2/2∫t1t2Ef2uθdu,≤θ−1θ2θ/2t2−t1θ/2F2θ.
Then, for0<t1<t2<∞, t2−t1≤1,1/θ+1/κ=1, we have
(50)ESt2−St1θ=E∫t1t2f1udu+∫t1t2f2udB1uθ,≤2θ−1E∫t1t2f1uduθ+2θ−1E∫t1t2f2udB1uθ,≤2θ−1t2−t1θ/κE∫t1t2f1uθdu+2θ−1θ−1θ2θ/2t2−t1θ/2F2θ,≤2θ−1t2−t1θ/κF1θt2−t1+2θ−1θ−1θ2θ/2t2−t1θ/2F2θ,=2θ−1t2−t1θF1θ+2θ−1θ−1θ2θ/2t2−t1θ/2F2θ,=2θ−1t2−t1θ/2t2−t1θ/2F1θ+θ−1θ2θ/2F2θ,≤2θ−1t2−t1θ/21+θ−1θ2θ/2Fθ,where Fθ=maxF1θ,F2θ. By Lemma 5, for each exponent υ∈0,θ−2/2θ, we get that each sample path of St is uniformly and locally Hölder continuous and which shows the uniform continuity of each sample path of St on t∈R+3. Similarly, the uniform continuity of It and Pt is proved on R+3. Therefore, we get the uniform continuity of each sample path of St,It,Pt to system (23) on t≥0.


## 5. Long Time Behavior of System

Here, we look at the solution behaviour of system (23) as time becomes very large. For that, we define the hypotheses which are useful in further analysis.
(51)H1:LηmaxrK,αm+12maxσ12,σ22,σ32,<minr−μ,r−d,H2:r−σ122<0,λ−μ−σ222<0,α+β−d−σ322<0.First, we will prove stochastic permanence which plays an essential part in population dynamics. We discuss this property as follows:Theorem 13.
If the assumptionH1 holds, then system (23) is stochastically permanent.Proof.
For any initial valueW0=S0,I0,P0∈R+3, we show that there exists a solution Wt=St,It,Pt such that
(52)limsupt⟶∞E1Wtγ≤M,where arbitrary nonnegative constant γ satisfies
(53)LηmaxrK,αm+γ+12maxσ12,σ22,σ32,<minr−μ,r−d.
By (53), there is a positive constant ρ>0 such that
(54)minr−μ,r−d−ρ−LηmaxrK,αm−12maxσ12,σ22,σ32,>0.
DefineVS,I,P=S+I+P for S,I,P∈R+3 and ZS,I,P=1/VS,I,P; from the Itô formula, we have
(55)dVS,I,P=Sr1−S+IK−λI−αPm+S+I−μ+λS−βPa+I+P−d+βIa+I+αSm+Sdt+σ1SdB1t+σ2IdB2t+σ3PdB3t,dZW=−Z2WSr1−S+IK−λI−αPm+S+I−μ+λS−βPa+I+P−d+βIa+I+αSm+Sdt+Z3Wσ1S2+σ2I2+σ3P2dt−Z2Wσ1SdB1t+σ2IdB2t+σ3PdB3t,=LZWdt−Z2Wσ1SdB1t+σ2IdB2t+σ3PdB3t.
Under the hypothesisH1, we introduce γ as a positive constant such that condition (53) holds. By Itô’s formula, we obtain
(56)L1+ZWγ=γ1+ZWγ−1LZW+12γγ−11+ZWγ−2Z4W×σ1S2+σ2I2+σ3P2.
Next, we chooseρ>0 to be small such that condition (54) holds. Then,
(57)Leρt1+ZWγ=ρeρt1+ZWγ+eρtL1+ZWγ,=eρt1+ZWγ−2ρ1+ZW2+B,where
(58)B=−γZ2WSr1−S+IK−λI−αPm+S+I−μ+λS−βPa+I+P−d+βIa+I+αSm+S−γZ3WSr1−S+IK−λI−αPm+S+I−μ+λS−βPa+I+P−d+βIa+I+αSm+S+γZ3Wσ1S2+σ2I2+σ3P2+γγ+12Z4Wσ1S2+σ2I2+σ3P2.
The upper bound of the function1+ZWγ−2ρ1+ZW2+B is defined in the following way:
(59)γZ3Wσ1S2+σ2I2+σ3P2≤γZWmaxσ12,σ22,σ32,γγ+12Z4Wσ1S2+σ2I2+σ3P2≤γγ+12Z2Wmaxσ12,σ22,σ32.
Hence,(60)Leρt1+ZWγ=eρt1+ZWγ−2ρ1+ZW2+B,≤eρt1+ZWγ−2ρ+2ρ−γminr−μ,r−d+γLηmaxrK,αm+γmaxσ12,σ22,σ32ZW+ρ−γminr−μ,r−d+γLηmaxrK,αm+γγ+12maxσ12,σ22,σ32Z2W.
From (53) and (54), we get a nonnegative constant Q satisfying Leρt1+ZWγ≤Qeρt. This implies that
(61)Eeρt1+ZWγ≤1+ZWγ+Qeρt−1ρ.
Therefore,(62)limsupt⟶∞EZγWt≤limsupt⟶∞E1+ZWtγ≤Qρ.
Note that,(63)S+I+Pγ≤3γS3+I3+P3γ/3=3γWγ,where W=S,I,P∈R+3. Accordingly,
(64)limsupt⟶∞E1Wtγ≤3γlimsupt⟶∞EZγWt≤3γQρ≔R.
Assume thatφ=v/R1/γ for any ν>0; then, by Chebyshev’s inequality, we get
(65)P∣Wt∣<φ=PWt−γ>φ−γ≤EWt−γ/φ−γ=φγEWt−γ,that is,
(66)liminft⟶∞P∣Wt∣≥φ≥1−ν.
Similarly, we can getχ>0 for any ε>0 such that liminft⟶∞P∣Wt∣≤χ≥1−ν. Hence, by Definition 7, system (23) is stochastically permanent.In population dynamics, there is a chance to lose the species population fully. So, the study of disappearance of species is also much important in the ecosystem.Theorem 14.
If the assumptionH2 holds, then the solution Wt=St,It,Pt of system (23) will be extinct with probability one for any given initial value W0=S0,I0,P0∈R+3.Proof.
DefineV4=lnS. By Itô’s formula, we obtain
(67)dlnS=r1−S+IK−λI−αPm+S−12σ12dt+σ1dB1t.
Integrating it from 0 tot gives
(68)lnSt=lnS0+r−σ122t−rK∫0tSudu−rK+λ∫0tIudu−α∫0tPum+Sudu+σ1∫0tdB1u.
Then,(69)lnSt≤lnS0+r−σ122t−rK∫0tSudu+σ1B1t.
Divide the above inequality on two sides byt and taking t⟶∞, we get
(70)limsupt⟶∞lnStt≤r−σ122<0,almost surely.
Also, define the Lyapunov functionV5=lnI; use Itô’s formula to obtain
(71)dlnI=−μ+λS−βPa+I−12σ22dt+σ2dB2t.
Integrating this from 0 tot, we have
(72)lnIt=lnI0−μ+σ222t−λ∫0tSudu−β∫0tPua+Iudu+σ2∫0tdB2u.
Consequently,(73)lnIt≤lnI0+λ−μ−σ222t+σ2B2t.
Dividing above inequality byt and taking t⟶∞, we obtain
(74)limsupt⟶∞lnItt≤λ−μ−σ222<0,almost surely.
Similarly, we define the Lyapunov functionV6=lnP and by Itô’s formula, we get
(75)dlnP=−d+βIa+I+αSm+S−12σ32dt+σ3dB3t.
Integrating it, we get(76)lnPt=lnP0−d+σ322t+β∫0tIua+Iudu+α∫0tSum+Sudu+σ3∫0tdB3u.
Accordingly,(77)lnPt≤lnP0+α+β−d−σ322t+σ3B3t.
Then, dividing this byt and taking t⟶∞, we get
(78)limsupt⟶∞lnPtt≤α+β−d−σ322<0,almost surely.
Hence, the desired assertion is proved.


## 6. Stochastic Asymptotic Stability

In this section, we prove that both planar and coexistence equilibrium of system (23) are stochastically asymptotically stable under certain assumptions.Theorem 15.
The equilibrium solutionE1=K,0,0 of system (23) is stochastically asymptotically stable on D, if
(79)R0≤1,R1≤1,12σ12KS2≤η1S,I,P,where η1S,I,P=r/KS−K2−λK−μI−αK/m+K−dP.Proof.
We use Theorem (1) and construct the suitable Lyapunov function as follows:
(80)VS,I,P=S−K−KlnSK+I+P.
We define the infinitesimal generatorL on the Lyapunov function and get
(81)LVS,I,P=S−Kr1−S+IK−λI−αPm+S+−μI+λIS−βIPa+I+−dP+βIPa+I+αPSm+S+12σ12KS2,=S−Kr1−S+IK−λIS−K−αPSm+S+αPKm+K+−μI+λIS−βIPa+I+−dP+βIPa+I+αPSm+S+12σ12KS2,≤S−Kr1−SK+λK−μI+αKm+K−dP+12σ12KS2,=−η1S,I,P+12σ12KS2.
IfR0≤1,R1≤1, it follows that η1S,I,P≥0. By assumptions, LVS,I,P<0 on D. Therefore, E1 is stochastically asymptotically stable.Theorem 16.
The equilibrium solutionE2=S2,I2,0 of system (23) is stochastically asymptotically stable on D, if
(82)d>d1,12σ12S2S2+12σ22I2I2≤η2S,I,P,where η2S,I,P=3r/KS−S22+r/KI−I22−αS2/m+βI2/a−d and d1=mrβKλ−μ+aαμKλ+r/amλKλ+r.Proof.
Consider the Lyapunov function for the equilibrium solutionE2(83)VS,I,P=S−S2−S2lnSS2+I−I2−I2lnII2+P.
Define the infinitesimal generatorL on the above function as follows:
(84)LVS,I,P=S−S2r1−S+IK−I−αPm+S+I−I2−μ+λS−βPa+I+−dP+βIPa+I+αPSm+S+12σ12S2S2+12σ22I2I2,=S−S2rS2+I2K−S+IK−λI−I2−αPm+S+I−I2λS−S2−βPa+I+−dP+βIPa+I+αPSm+S+12σ12S2S2+12σ22I2I2,=−rKS−S22−rKS−S2I−I2+αS2m+S+βI2a+I−dP+12σ12S2S2+12σ22I2I2,≤−rKS−S22−rKS−S22+I−I222+αS2m+S+βI2a+I−dP+12σ12S2S2+12σ22I2I2,≤−3r2KS−S22−r2KI−I22+αS2m+S+βI2a+I−dP+12σ12S2S2+12σ22I2I2,=−η2S,I,P+12σ12S2S2+12σ22I2I2.
Ifd>d1 where d1=mrβKλ−μ+aαμKλ+r/amλKλ+r, we obtain η2S,I,P≥0. By assumptions, LVS,I,P<0 on D. Therefore, E2 is stochastically asymptotically stable.Theorem 17.
The equilibrium solutionE3=S3,0,P3 of system (23) is stochastically asymptotically stable on D, if
(85)m1<m<m2and12σ12S3S2+12σ32P3P2≤η3S,I,P,wherem1=d−αdK+2r−2Kr/2dr, m2=μα−d/dλ, and η3S,I,P=r/K−αP3/m+S3+αS3/2m+S3S−S32+αS3/2m+S3P−P32−λS3−μI.Proof.
The Lyapunov function for the equilibrium solutionE3 can be written as
(86)VS,I,P=S−S3−S3lnSS3+I+P−P3−P3lnPP3.
Define the infinitesimal generatorL on VS,I,P as
(87)LVS,I,P=S−S3r1−S+IK−λI−αPm+S+−μI+λIS−βIPa+I+P−P3−d+βIa+I+αSm+S+12σ12S3S2+12σ32P3P2,≤S−S3r1−SK+S−S3−λI−αPm+S+−μI+λIS−βIPa+I+P−P3−d+βIa+I+αSm+S+12σ12S3S2+12σ32P3P2,=−rKS−S32+S−S3−λI+αP3m+S3−αPm+S+−μI+λIS−βIPa+I+P−P3βIa+I+αSm+S−αS3m+S3+12σ12S3S2+12σ32P3P2,≤−rKS−S32+λS3−μI−αS−S3PS3−SP3m+Sm+S3+12σ12S3S2+12σ32P3P2,=−rKS−S32+λS3−μI+αP3S−S32m+Sm+S3−αS3S−S3P−P3m+Sm+S3+12σ12S3S2+12σ32P3P2,≤−rKS−S32+λS3−μI+αP3S−S32m+Sm+S3−αS3m+Sm+S3×S−S32+P−P322+12σ12S3S2+12σ32P3P2,=−rK−αP3m+Sm+S3+αS32m+Sm+S3S−S32−αS3P−P322m+Sm+S3+λS3−μI+12σ12S3S2+12σ32P3P2,≤−rK−αP3m+S3+αS32m+S3S−S32−αS3P−P322m+S3+λS3−μI+12σ12S3S2+12σ32P3P2,=−η3S,I,P+12σ12S3S2+12σ32P3P2.
One can see thatη3S,I,P≥0, if m1<m<m2 where m1=d−αdK+2r−2Kr/2dr and m2=μα−d/dλ. By assumptions, LVS,I,P<0 on D. Therefore, E3 of system (23) is stochastically asymptotically stable.Theorem 18.
The equilibrium solutionE4=S∗,I∗,P∗ of system (23) is stochastically asymptotically stable on D, if
(88)r1<r<r2,12σ12S∗S2+n12σ22I∗I2+n22σ32P∗P2≤η4S,I,P,where r1=1/2KαP∗/m+S∗−λn1−1/2, r2=1/2Kλn1−1/2−βP∗/a+I∗, and η4S,I,P=r/2K−αP∗/m+S∗+λn1−1/2S−S∗2+λn1−1/2−r/2K−n1βP∗/a+I∗I−I∗2,
forn1=am+S∗/ma+I∗ and n2=m+S∗/m.Proof.
Construct the suitable Lyapunov function for the coexistence equilibrium as follows(89)VS,I,P=S−S∗−S∗lnSS∗+n1I−I∗−I∗lnII∗+n2P−P∗−P∗lnPP∗,
whereni>0,i=1,2 are suitable constants to be determined in further analysis. Define the infinitesimal generator L on VS,I,P and we get
(90)LVS,I,P=S−S∗r1−S+IK−λI−αPm+S+n1I−I∗−μ+λS+βPa+I+n2P−P∗−d+βIa+I+αSm+S+12σ12S∗S2+n12σ22I∗I2+n22σ32P∗P2,=S−S∗−rKS−S∗−rKI−I∗−λI−I∗+αP∗m+S∗−αPm+S+n1I−I∗λS−S∗+βP∗a+I∗−βPa+I+n2P−P∗βIa+I−βI∗a+I∗+αSm+S−αS∗m+S∗+12σ12S∗S2+n12σ22I∗I2+n22σ32P∗P2,=−rK−αP∗m+Sm+S∗S−S∗2+n1βP∗a+Ia+I∗I−I∗2−λn1−1−rKI−I∗S−S∗+−n1β+n2βaa+I∗I−I∗P−P∗a+I+−α+n2mαm+S∗P−P∗S−S∗m+S+12σ12S∗S2+n12σ22I∗I2+n22σ32P∗P2,≤−rK−αP∗m+S∗S−S∗2+n1βP∗a+I∗I−I∗2−λn1−1−rK×I−I∗2+S−S∗22+−n1β+n2βaa+I∗I−I∗P−P∗a+I+−α+n2mαm+S∗P−P∗S−S∗m+S+12σ12S∗S2+n12σ22I∗I2+n22σ32P∗P2,=r2K−αP∗m+S∗+λn1−12S−S∗2+λn1−12−r2K−n1βP∗a+I∗×I−I∗2+−n1β+n2βaa+I∗I−I∗P−P∗a+I+−α+n2mαm+S∗P−P∗S−S∗m+S+12σ12S∗S2+n12σ22I∗I2+n22σ32P∗P2,=−η4S,I,P+12σ12S∗S2+n12σ22I∗I2+n22σ32P∗P2.
One can observe thatη4S,I,P≥0, if r1<r<r2 where r1=1/2KαP∗/m+S∗−λn1−1/2 and r2=1/2Kλn1−1/2−βP∗/a+I∗. By assumptions, LVS,I,P<0 on D. Therefore, E4 of system (23) is stochastically asymptotically stable.


## 7. Numerical Simulation

From the practical point of view, numerical results play a major role in population dynamics. Here, we have provided graphical solution curves to carry out and check our theoretical results of system (23). To find the strong solution of system (23) with given initial value, Higham [34] introduced Milstein method. By using this method, we have derived the following discretization equations from our given system (23):
(91)Sk+1=Sk+Skr1−Sk+IkK−λIk−αPkm+SkΔt+σ1SkΔtξk+σ122Skξk2−1Δt,Ik+1=Ik+IkλSk−βPka+Ik−μΔt+σ2IkΔtηk+σ222Ikηk2−1Δt,Pk+1=Pk+PkβIka+Ik+αSkm+Sk−dΔt+σ3PkΔtζk+σ322Pkζk2−1Δt,where Δt is a positive time increment parameter and ξk,ηk,ζk,k=1,2,3,⋯,n represent independent random Gaussian variables N0,1, it is given with the help of pseudorandom number generators. Using MATLAB software and possible parameters, we have generated the solution curves of the deterministic system (11) and corresponding stochastic system (23).In Figures1–3, we take the initial conditions S0=0.1, I0=0.1, and P0=0.1 and parameters r=1, K=3, β=1, α=1.15, c=0.7, d=1, a=1, m=1, and λ=0.5. The only difference between Figures 1–3 is the intensities of environmental changes. Specially, we choose σ1=σ2=σ3=0.1 in Figure 1; in Figure 2, we have chosen σ1=σ2=σ3=0.04 and in Figure 3σ1=σ2=σ3=0.01. From Figures 1–3, we have observed that the coexistence equilibrium S∗,I∗,P∗ solution curves of stochastic model (23) always oscillates with respect to the curves of the deterministic model (11). From Figures 1–3, we have arrived at the conclusion that with the decrease in the values of σ1,σ2,σ3 the fluctuations of the solution curves of the stochastic system are reduced and coincide with that of the solution curves of the deterministic system.Figure 1
Simulation showing the solution curves of deterministic system (11) and corresponding stochastic system (23) with large noise σ1=σ2=σ3=0.1.Figure 2
Simulation showing the solution curves of deterministic system (11) and corresponding stochastic system (23) with noise σ1=σ2=σ3=0.04.Figure 3
Simulation showing the solution curves of deterministic system (11) and corresponding stochastic system (23) with small noise σ1=σ2=σ3=0.01.In Figure4, the parametric values are r=0.09, K=4, β=0.09, α=1.1, c=0.02, d=0.8, a=1, m=1, and λ=0.005 with initial values S0=1.5,I0=0.8,P0=0.6. By Theorem 14, we observe that all the species in the stochastic case (23) disappear completely under the hypothesis H1 and we have seen this graphically via Figure 4. Further, we have simulated the corresponding deterministic case in Figure 4 which satisfies condition (14) and leading to the permanence of each species. Under parametric values r=0.001, K=2.5, β=0.01, α=0.04, c=0.02, d=0.01, a=1, m=1, and λ=0.003 and initial parametric condition S0=2,I0=1.3,P0=0.5, by Theorem 14, all the species of the stochastic system (23) as well as corresponding deterministic system (11) disappear fully under the hypothesis H1 and the solution curves are shown in Figure 5. If the deterministic system goes to extinction, the stochastic system also goes to extinction, in which case solution is not dependent on the intensity of noise value (Compare Figure 4 with Figure 5).Figure 4
Simulation showing that the solutions of the stochastic system (23) become extinct with σ1=0.5,σ2=0.2,σ3=0.3 and the solutions of the deterministic system (11) are persistent.Figure 5
Simulation showing that the solution curves of both deterministic system (11) and corresponding stochastic system (23) are extinct with σ1=0.5,σ2=0.2,σ3=0.3.Furthermore, in Figure6, parameters chosen are r=3, K=3, β=1.5, α=1, c=1, d=1, a=1, m=1, and λ=1.7 with initial values S0=1.5,I0=0.8,P0=0.1 and simulations also confirm the result that all three species of system (11) and (23) are permanent by Theorem 13.Figure 6
Simulation showing that the solution curves of both deterministic system (11) and corresponding stochastic system (23) are permanent with σ1=0.05,σ2=0.04,σ3=0.06.Finally, Figure7 illustrates the Theorem 18. We choose the parameters r=2.5, K=3, α=1.5, β=0.9, c=0.8, d=1, m=1, a=1, and λ=1.2 with different initial values satisfying condition (88). Theorem 18 says that S∗,I∗,P∗ of (23) is stochastically asymptotically stable.Figure 7
Simulations showing the global stability convergence of the stochastic system (23) around the coexistence equilibrium solution S∗,I∗,P∗ with different initial values and σ1=σ2=σ3=0.01


## 8. Concluding Remarks

In reality, predator-prey system reacts to more consequences in the environmental changes in the ecosystem. So, the study of dynamical behavior of the predator-prey model has been receiving more attention, and many efforts have been taken in the field of population dynamics by several authors. To overcome the effect of random fluctuations and control the disease in an ecosystem, we have proposed the stochastic Holling type II predator-prey model (23) with disease in the prey. We are interested to know the changes of dynamical behavior of model (23) with respect to the intensities of environmental fluctuations. There is also need to reveal the relationship between the intensities of the environmental fluctuations and parameters associated with the system since the predator-prey model is disturbed by environmental noises.First, we explain the basic assumptions of our model (11) with boundedness of solutions. We have introduced the stochastic term to the deterministic model (11) and proved that for any positive initial value there exists a positive global solution of the stochastic system (23). For a stochastic predator-prey model, the boundedness of a solution is also verified because it gives a guarantee of system validity. In Theorem 12, we have checked the property of uniform continuity of the positive solution of system (23). Permanence and extinction property are provided for system (23) since as time tends to be large it validates the long time and short time survival in a ecosystem. Under assumption H1, our stochastic system (23) attains the permanence. The stochastic system goes to extinction state if the assumption H2 holds. Furthermore, with the help ofthe Lyapunov function theory we have examined the stochastic asymptotic stability of all feasible equilibrium solutions of (23). The parametric restriction for stochastic stability is established. (see the inequalities in Theorem 15, 16, 17, and 18).In a stochastic system, we have observed that if the intensities of the environmental changes increase, the fluctuations of the solutions increase. The numerical result shows that if the strong environmental changes happen in an ecosystem making each species to disappear it can occur in reality while the stochastic system (23) can maintain permanence under sufficiently small environmental noise. In permanence, if we decrease the level of environmental fluctuations, solutions become stable. By the above fact, the stochastic results give more accurate results compared to the deterministic results.



---

*Source: 3632091-2020-09-04.xml*