Stationary distribution and global stability of stochastic predator-prey model with disease in prey population

In this paper, a new stochastic four-species predator-prey model with disease in the first prey is proposed and studied. First, we present the stochastic model with some biological assumptions and establish the existence of globally positive solutions. Moreover, a condition for species to be permanent and extinction is provided. The above properties can help to save the dangered population in the ecosystem. Through Lyapunov functions, we discuss the asymptotic stability of a positive equilibrium solution for our model. Furthermore, it is also shown that the system has a stationary distribution and indicating the existence of a stable biotic community. Finally, our results of the proposed model have revealed the effect of random fluctuations on the four species ecosystem when adding the alternative food sources for the predator population. To illustrate our theoretical findings, some numerical simulations are presented.


Introduction
Mathematical modelling in ecology and epidemiology has become an important tool in analysing the spread and control of infectious diseases.Anderson and May [3] were the first to merge the above two fields and formulate a predator-prey model where some disease infected the prey species.Most infected disease models are created from the classical SIR model of Kermack and McKendrick [18].Subsequently, many researchers have proposed and studied different predator-prey models in the presence of a disease.Studies in the field of eco-epidemiology have increased enormously in the last two decades, and we only mention a few of them [16,20,21,31,32,37].Qin et al. performed Holling III predator-prey model with an infectious predator.A functional response plays a key role in the predator-prey model and describes the number of prey consumed per predator per unit of time for some given time.Here, we used the Beddington-DeAngelis functional response for the predator and the infected first prey interaction.Mathematically, we may consider both prey-dependent and ratio-dependent models as limiting cases of the general Beddington-DeAngelis type predator-prey model.Much progress [10,12,13,15,35,36,39,40] have been seen in the study of the predator-prey model with Beddington-DeAngelis functional response.Recently, Arora et al. [1] proposed a predator-prey system highlighting the significance of additional food for predators with a Beddington-DeAngelis functional response.
Another important factor in ecological and eco-epidemiological species interaction is alternative food sources for the predator population.Generally, predators do not depend on a single prey species for food.Predators will change to alternative food when the density of their preferred prey is low.Fryxell and Lundberg [11] has argued that alternative food plays an essential role in promoting the persistence of predator-prey systems.Many authors [1,9,45] have considered one predator-two prey systems to show the importance of alternative food for the coexistence of species.Thus, a better understanding of the mechanisms, assumptions and biological hypotheses of predator-prey dynamic models involving alternative food is explicitly needed when the disease is present in the prey population.
Here, we consider a model with one predator and two prey populations with the availability of alternative food sources for the predator population when a disease has affected only the first prey population.Srinivasu et al. [38] showed that additional food serves as a tool for the biological conservation of predator-prey systems involving type III functional response.
In reality, a population model can be affected by environmental noise; hence, stochastic population systems play a major role in various branches of applied sciences, including population dynamics and biology.In real life, continuous fluctuations in the environment change the population to a lower or higher extent level.May [30] discovered that all the parameters involved in the population system exhibit random fluctuation as the factors controlling them are not constant.Therefore most authors [5,14,25,41,46] introduce stochastic perturbation into the deterministic models and make the model more realistic for examination; the resulting dynamics solely depend on the intensity of environmental fluctuations.Precisely, enormous efforts have been taken to analyse the dynamical behaviour with environmental fluctuations [8,21,28,34,39,40,42,44,45].Cao et al. [5] transformed the stochastic system with a weak kernel case into an equivalent degenerate system and studied its behaviour.[14] et al. proposed the predator-prey model with Holling II increasing function in the predator.Lu [27] presented Crowley-Martin predator-prey model in a stochastic environment and analysed asymptotical stability in probability by combining the Khasminskii theory of stability with the Lyapunov method.Zhang et al. [43] investigated the influence of anti-predator behaviour due to the fear induced by the predators for the stochastic eco-epidemiological predator-prey model.Recently, Liu et al. [24] developed and analysed a Leslie-Gower stochastic predator-prey model with impulsive toxicant input in polluted environments.Liu [26] carried out several two-species predatorprey models and showed the impact of randomness.A small amount of work has been carried out with the stochastic four-species predator-prey model and disease in the prey.
The present article is organized as follows.Section 2 presents some known definitions, theorems, and lemmas that will be used in the following analysis.In Section 3, we discuss the detailed explanation about the formulation of the stochastic model and analyse the existence and uniqueness of a global positive solution of the system.In Section 4, stochastic permanence and extinction are established under certain parametric restriction.We construct a suitable Lyapunov function to examine the global asymptotic stability for positive equilibrium in Section 5.In Section 6, we establish the conditions for unique stationary distribution.Additionally, we generate some figures to illustrate our main results in Section 7. Finally, we end the paper with a discussion in Section 8.

Preliminaries
Consider the d-dimensional stochastic model (SM) of the form with an initial value The differential operator L associated with the SM (1) is given by In addition to the existence and uniqueness assumptions, we assume that f and g satisfy f (x * , t) = 0 and g(x * , t) = 0 for an equilibrium solution x * , for t ≥ t 0 .Moreover we assume that x 0 is a non-random constant with probability 1.

Definition 2.1 (see [22]):
The solution X(t) of ( 1) is said to possess stochastic permanent if for any ε ∈ (0, 1), there exists a pair of positive constants δ = δ(ε) and χ = χ(ε) such that for any initial value X 0 ∈ R 4 + , the solution X(t) to (1) has the properties Theorem 2.1 (see [2]): Assume that f and g satisfy the existence and uniqueness assumptions and that they have continuous coefficients with respect to t.
(i) Suppose that there exists a positive definite function Then, the equilibrium solution x * of (1) is stochastically stable.(ii) If, in addition, V is decreasing and there exists a positive definite function then the equilibrium solution x * is stochastically asymptotically stable.
Let X(t) be a homogeneous Markov process defined in the l dimensional Euclidean space E l and described by the following system of stochastic differential equations: The diffusion matrix is defined by Hasminskii [19], We assume there exists a bounded domain U ⊂ E l with regular boundary , having the following properties: P1: In some neighbourhood of a domain U, the smallest eigenvalue of the diffusion matrix A 4 (x) is bounded away from zero.P2: If x ∈ E l \U, the mean time τ at which a path emerging from x reaches the set U is finite, and sup x∈S E x τ < ∞ for every compact subset S ⊂ E l .
Lemma 2.1: If above assumptions hold, then the Markov process X(t) has a stationary distribution ϕ(•).Let g(•) be a function integrable with respect to the measure ϕ.Then for all x ∈ E l .

Mathematical model
In this section, we develop the predator-prey model with two prey and one predator population where the disease infects only the first prey population.Sharma and Samanta [37] pay attention to a similar model and analyse the global stability of the positive equilibrium.We develop the model as follows: with initial densities Here, X S (t), X I (t), Y(t) and Z(t) denote the numbers of susceptible first prey, infected first prey, second prey and predator, respectively, and all parameters in the model are positive.We choose the Beverton-Holt functional response to represent the birth rate of the prey population.In the absence of disease and predator population, susceptible first prey population will follow the growth equation where is the maximum per-capita growth rate of prey and b is per-capita mortality rate; moreover, X S a+X S can be represented as − a a+X S , which approaches as X S goes to infinity.The parameters used in the system modelling with their meaning are listed in the Table 1.
The given model requires the following assumptions: (1) Only the first prey population is infected by the disease and because of the presence of the disease, the first prey population is divided into two sub-classes, namely (i) susceptible first prey population (X S ) and (ii) infected first prey population (X I ).
(2) The infected first prey population dies before having the capability of reproducing and does not recover from the disease.(3) Disease does not spread to predators from the infected first prey.(4) Predators do not consume the susceptible first prey population.It has more possibilities to consume the infected first prey population as they are easily available to be preyed upon owing to weakening due to the infection, i.e. they lose their power to protect themselves from a predator.The predators also consume the second prey because of their significant availability in the total population.Thus, our assumption is reasonable in such cases.This explains why the infected first prey and second prey increase their availability to predation [21].(5) The interaction between the infected first prey and predator is assumed to be governed by a Beddington-DeAngelis type functional response.This functional response was introduced by Beddington [4] and DeAngelis et al. [7].The interaction between the second prey and predator is considered a Lotka-Volterra functional response.There is no interaction between the first prey and second prey population.
(6) The growth of the second prey population is exponential, so the availability of second prey becomes very large when the predator becomes extinct.Because of this, there is no need to consider the search time for the second prey population.
In reality, environmental noises significantly affect population dynamics, which are essential components in an ecosystem.Taking into account the effect of a randomly fluctuating environment, we included white noise in each equation of the system (2).Let us assume that the fluctuations in the environment incorporate themselves mainly as fluctuations in the death rate of the susceptible first prey population, the death rate of the infected first prey population, the growth rate of the second prey population, and the death rate of the predator population, as follows: where B 1 (t), B 2 (t), B 3 (t), B 4 (t) represent the independent Brownian motions and σ 1 , σ 2 , σ 3 , σ 4 represent the intensities of the white noises.Therefore the corresponding stochastic system can be described by the Itô equations In the past several years, no work has been done on the given stochastic predator-prey system (4).Our aim is to analyse the dynamical properties of the stochastic model and show the impact of environmental noise in the population system.
Theorem 3.1: All solutions of the system (2) with initial conditions Proof: The proof of this theorem is similar to the proof of Theorem 1 in [37].
Proof: Consider the following transformation of variables Let us use Itô's formula Using this, we obtain and dp = Lp dt + g(t, p) dB(t).
Similarly, from system (4), we obtain . Now, the functions associated with the above system have an initial growth and satisfy a local Lipchitz condition.Therefore there exists a unique local solution (p(t), q(t), r(t), s(t)) defined in some interval [0, τ e ).Consequently, X S (t) = e p(t) , X I (t) = e q(t) , Y(t) = e r(t) , Z(t) = e s(t) is the unique positive local solution of (4) starting from an interior point of the first quadrant.
Theorem 3.3: For any initial condition (X S 0 , X I 0 , Y 0 , Z 0 ) ∈ Int(R 4 + ), the system (4) possesses a unique solution (X S (t), X I (t), Y(t), Z(t)) for t ∈ [0, τ e ) and the solution remains in Int(R 4  + ) with probability one; therefore, for all t ≥ 0, (X S (t), Proof: To prove that this solution is global, it is enough to show that τ = ∞ almost surely.Let k 0 be a sufficiently large nonnegative integer such that the closed ball B(k 0 ) ∈ R 4 + contains (X S 0 , X I 0 , Y 0 , Z 0 ).We choose for any k ≥ k 0 and define the stop-time by using the ideas in [17] as where inf ∅ = ∞ (∅ being the empty set).Clearly τ k is increasing as k → ∞.
( 6 ) Thus, by denoting k = {τ k ≤ T}, there is an integer k 1 ≥ k 0 such that, for all k ≥ k 1 , Define a C 4 function V : where the function V(X S , X I , Y, Z) is positive definite for all (X S , X I , Y, Z) ∈ Int(R 4 + ).Using Itô's formula, we obtain Taking the differential of V(S, I, P) along the solution trajectories of the system (4) using above equation, we obtain Define three positive constants We utilize Lemma 2.2 (the proof is given in [6]) in the form where where Using the properties of Itô's integral and with the above inequality we obtain that Taking the expectations of both sides of the above inequality and applying Fubini's theorem, we obtain By Grownwall's inequality given in [29] and the above inequality, we obtain where From ( 7), we have P( k ) ≥ .It is noteworthy that there is at least one of X S (τ k , ω), X I (τ k , ω), Y(τ k , ω) and Z(τ k , ω) belongs to the set {k, Hence from ( 6) and ( 8), it follows that + almost surely.The proof is complete.
Theorem 3.3 tells us that the solution of system (4) will remain in R 4 + .With this property, we continue to discuss how the solution varies in R 4  + in more detail.

Theorem 3.4:
The solutions of the system (4) are stochastically ultimately bounded for any initial value W 0 = (X S 0 , Proof: This proof is similar to the proof of Theorem 3 in [44].

Long time behaviour of system
Here, we look at the long-time behaviour of the system (4), such as permanence and extinction to determine how the solution varies in R 4 + in more detail.First, we impose the following hypotheses, which will be useful later on.
(H1): max{b, c, d, e} In population dynamics, stochastic permanence is one of the most essential properties for species survival.We discuss this property with the fruitful ideas of [28] as follows: Theorem 4.1: Suppose that assumption (H1) holds; then system (4) will be stochastically permanent.
Proof: First, we show that for any initial value where κ is an arbitrary positive constant satisfying max{b, c, d, e} By (9), there exists an arbitrary positive constant ρ > 0 satisfying ; by Itô's formula, we have and Under Assumption (H1), we can certainly choose a positive constant κ satisfying (9).Again, applying Itô's formula, we obtain Then, we can choose ρ > 0 sufficiently small such that it satisfies (10); consequently where Let us now discuss the upper boundedness of the function and Hence From ( 9) and ( 10), we know that there exists a positive constant Q such that Le ρt (1 + U(W)) κ ≤ Qe ρt which implies that Note that Thus, for any ε > 0, let ν = ( ε M ) 1 κ .By Chebyshev's inequality, we obtain Similarly, we obtain that for any ε > 0, there exists χ > 0 such that lim inf t→∞ P{|W(t)| ≤ χ } ≥ 1 − ε.Hence the theorem is proved.
In ecology, the worse occurs when a species disappears completely.Here we prove that if the noise is sufficiently large, the solution will become extinct with probability one.Theorem 4.2: Assume that (H2) holds, then for any given initial value W(0) = (X S (0), X I (0), Y(0), Z(0)) ∈ R 4  + , the solution W(t) = (X S (t), X I (t), Y(t), Z(t)) of system (4) will be extinct with probability one.
Proof: Define the Lyapunov function V 1 = ln X S .Then, by Itô's formula, we have Integrating it from 0 to t yields ln X S (t) ≤ ln X S (0) Dividing by t on both sides and letting t → ∞, we obtain lim sup Furthermore, define the Lyapunov function V 2 = ln X I ; by Itô's formula, we have Integrating it from 0 to t, we obtain Dividing by t and letting t → ∞, we obtain lim sup Define the Lyapunov function V 3 = ln Y; by Itô's formula, we arrive at Integrating it, we get ln Then, dividing it by t and letting t → ∞, we obtain Similarly, we define the Lyapunov function V 4 = ln Z and by Itô's formula, we get Integrating it, we have ln Dividing it by t and letting t → ∞, we have lim sup The desired assertion is thus proved.

Stochastic asymptotic stability
Using the ideas of Pitchaimani and Rajaji [33], we shall analyse the stochastic stability, which will lead to certain conditions.Here, the discussion is entirely about the mathematical theory on stochastic asymptotic stability of the positive equilibrium solution (X * S , X * I , Y * , Z * ) of the system (4).The positive equilibrium point E * = (X S * , X I * , Y * , Z * ) where Z * = γ −d δ and others are obtained as follows: Let X I * be a positive root of the following quadratic equation where The roots of the above quadratic equation are When any one of the following cases is satisfied, then the equilibrium points X I * can have one or two positive values.We have used the criteria in Theorem 2.1 and the ideas of [23] to prove the stochastic system stability.

Theorem 5.1:
The equilibrium solution E * = (X S * , X I * , Y * , Z * ) of the system (4) is stochastically asymptotically stable, if Proof: Consider the Lyapunov function The infinitesimal generator L acting on the Lyapunov function can be written as follows One can observe that φ(X S , By assumption, LV(X S , X I , Y, Z) becomes negative definite; therefore the equilibrium solution E * of the system (4) is stochastically asymptotically stable.

Existence of stationary distribution
Here, we shall give the existence of the stationary distribution of prey and predator populations.For this purpose, we must find the stationary distribution for solutions of the system (4), which implies stability in the stochastic sense.We use lemma 2.1, (that is given in [19]) to prove the main theorem on the stationary distribution.Theorem 6.1: ) simultaneously.Then there exists a stationary distribution ϕ(•) for the stochastic system (4).

Proof:
The proof of this theorem is based on the following positive definite function V : Applying Itô's formula, we obtain and where Here, the ellipsoid lies entirely in Int(R + 4 ).Let U to be a neighbourhood of the ellipsoid with Ū ⊆ E 4 = IntR + 4 , where Ū is the closure of U.Then, for every x ∈ U\E 4 , we have LV < 0, which implies condition (P2) in Lemma 2.1 is satisfied.
In addition, there exist for all (X S , X I , Y, Z) ∈ Ū, ξ ∈ R 4 , which implies that condition (P1) in above Lemma 2.1 is satisfied.Therefore we can conclude from the above discussion that the stochastic system (4) has a stationary distribution ϕ(•).

Numerical simulation
In this section, let us consider some examples to visualize our theoretical findings by simulating the stochastic system (4).To find the strong solution of system (4) with a given initial value, we use the Milstein method.Consider the following discretization equations of model ( 4) where time increment t > 0 and ξ k , η k , ζ k , k , k = 1, 2, 3, . . ., n, are independent random Gaussian random variables N(0, 1) that can be generated numerically by pseudo random number generators.We can construct suitable simulations for the stochastic system (4) and the corresponding deterministic system (2) with the help of MATLAB software and possible parameters.Choose the initial value as X S (0) = 0.6, X I (0) = 0.5, Y(0) = 0.5, Z(0) = 0.8 for all figures that are going to be discussed in this platform.
The only difference between Figure 1(a,b) is the intensity of white noise.From the figures, we see that the solution curve of the stochastic model (4) always fluctuates around the solution curves of the deterministic system (2).We observe that with the decreasing of σ 1 , σ 2 , σ 3 , σ 4 values, the dynamics of the stochastic system become increasingly similar to the deterministic system.For the parametric and intensity values σ 1 = 0.7, σ 2 = 0.6, σ 3 = σ 4 = 0.5 condition (H2) in Theorem 4.2 is verified by the observer; hence, Figure 2(a) shows that all the four species of the stochastic system (4), as well as the corresponding deterministic system (2), will die out.Furthermore, in Figure 2(b), the parameters chosen are the same as  those used in (12); moreover, the intensities σ 1 = σ 2 = 0.03, σ 3 = σ 4 = 0.01 satisfy condition (H1) in Theorem 4.1, and simulations also confirm that all three species of both systems are permanent.
Figure 3 illustrates Theorem 5.1.We choose the parameters values as in ( 12) with different initial values (0.3, 0.2, 0.4, 0.5), (0.2, 0.5, 0.4, 0.7), (0.5, 0.2, 0.6, 0.8) and (0.2, 0.3, 0.3, 0.4).Our model admits a stochastically asymptotically stable solution around the positive equilibrium for the above parameters.From this, we can conclude that the positive equilibrium agrees with Theorem 5.1 under the assumed conditions.The stochastic stability shown for each of the four species in the phase-portrait plane is presented in Figure 4.
The parameter values were chosen (12), and the choice of σ 1 , σ 2 , σ 3 and σ 4 was based on the conditions needed for the existence of a stationary distribution (see Theorem 6.1).We kept all parameter values fixed and repeated the simulation 10, 000 times but never observed any extinction scenario appear at any future time.Figure 5 shows that all species are distributed normally around the mean values (0.30, 0.85, 0.07, 1.27).After that, we increased the strengths of environmental forcing to σ 1 = 0.06, σ 2 = 0.06, σ 3 = 0.03, σ 4 = 0.03 and, again, we focussed on the species distribution fluctuate around the deterministic steady-state value; however, the amplitude of fluctuation is more comparable to that of the earlier case given in Figure 6.In this case, the susceptible first prey population is distributed within (0.17, 0.42); the infected first prey population occupies within (0.7, 1.1), the second prey lies within (0, 0.2), and the predator population remains within the range (0.79, 1.78).In the earlier case, the distribution of each species was within (0.22, 0.37), (0.75, 0.99), (0, 0.09) and (1, 1.5).Figures 5 and 6 show that in stochastic stability, all the solution curves will converge to one particular region.In deterministic stability, all the solution curves will converge to a single point.In Figures 5 and 6, the left-hand side shows that the mean solution of the stochastic system coincided with the solution curves of the deterministic system.Figure 7(a,b) present a clear representation of the stationary distribution for all species in a single graph, similar to Figures 5 and 6.

6.
The left panel represents the solution trajectories of both system from one simulation run; the right panel represents the stationary distribution of each species in the system (4) separately from 10,000 simulation run with intensity of noise σ 1 = σ 2 = 0.06, σ 3 = σ 4 = 0.03.

Discussion
From an applications viewpoint, the population model has become more interesting and studied extensively by several researchers.Our work is the first attempt to consider the stochastic four species predator-prey model with disease in the first prey (4).The main objective of the present paper is to study the model regarding how the intensities of environmental noises affect each population in the system (4) where an alternative food source is considered for predator population and disease has affected only the first prey population.Our analytical findings have pointed out that the intensity of noises can play a significant role in the survival of all species.These result can be verified with numerical support.
First, we explain the rigour of our model (2) under several assumptions.We introduce the stochastic term to the system (2) and show that the stochastic system (4) has a local and global positive solution for any positive initial value.The property of permanence and extinction is more desirable since it means the long-time and short time survival in a population dynamics.Under assumption (H1), our stochastic system (4) attains the permanence.The stochastic system goes to the extinction state if assumption (H2) holds.In real life, this may happen when a serious disease or severe weather occurs.Furthermore, we discuss the stochastic asymptotic stability of positive equilibrium with the help of the invariance principle and Lyapunov's second method that are given in Theorem 5.1.Thereafter, we have shown the existence of a stationary distribution for all species under certain parameter restriction discussed in Theorem 6.1.This parametric restriction follows the idea that a large intensity of noise can destabilize the system, and in this particular situation we cannot find a stationary distribution.
In a stochastic system, we see that larger intensities of the white noise lead to increasingly larger fluctuations of the solutions.The numerical result shows that strong noise may make the species extinct, which can occur in reality.In permanence, if we reduce the intensity of noise, solutions become stable.From this aspect, one can see the major difference when comparing the deterministic with stochastic models: that is, in the deterministic system the parameters are assumed to be constant, which is not true biologically in real life.

Figure 2 .
Figure 2. (a) This graph shows that the species of system (2) and (4) goes to extinction and (b) each species in both system goes to permanent.

Figure 4 .
Figure 4.The phase trajectories clearly shows the stochastic stability of individual species in the system (4) with different initial values and with the noise σ 1 = 0.04, σ 2 = 0.03, σ 3 = 0.02, σ 4 = 0.01 converges to the region where the positive equilibrium occur.

Figure 5 .
Figure 5.The left panel shows the solution trajectories of both deterministic and stochastic systems from one simulation run; the right panel shows the stationary distribution of each species in the system (4) separately from 10,000 simulation runs with intensity of noise σ 1 = σ 2 = 0.03, σ 3 = σ 4 = 0.01.

Table 1 .
Parameters and meaning -predator-prey model with disease in the first prey.