# Deterministic Sudden Changes and Stochastic Fluctuation Effects on Stability and Persistence Dynamics of Two-Predator One-Prey Model

Authors: Jawdat Alebraheem; Nasser S. Elazab; Mogtaba Mohammed; Anis Riahi; Ahmed Elmoasry
Journal: Journal of Mathematics
 (2021)

---

## Abstract

In this paper, we present new results on deterministic sudden changes and stochastic fluctuations’ effects on the dynamics of a two-predator one-prey model. We purpose to study the dynamics of the model with some impacting factors as the problem statement. The methodology depends on investigating the seasonality and stochastic terms which make the predator-prey interactions more realistic. A theoretical analysis is introduced for studying the effects of sudden deterministic changes, using three different cases of sudden changes. We show that the system in a good situation presents persistence dynamics only as a stable dynamical behavior. However, the system in a bad situation leads to three main outcomes as follows: first, constancy at the initial conditions of the prey and predators; second, extinction of the whole system; third, extinction of both predators, resulting in the growth of the prey population until it reaches a peak carrying capacity. We perform numerical simulations to study effects of stochastic fluctuations, which show that noise strength leads to an increase in the oscillations in the dynamical behavior and became more complex and finally leads to extinction when the strength of the noise is high. The random noises transfer the dynamical behavior from the equilibrium case to the oscillation case, which describes some unstable environments.


---

## Body

## 1. Introduction

### 1.1. Preface

Theoretical ecology has motivated many mathematicians to discuss different ideas and models from a purely mathematical standpoint; see for example [1–5]. Mathematical modeling is a useful tool to determine how a process works and to predict what may follow [3]. Many problems taken into consideration in mathematical ecology seem simple, but are considered complicated problems due to the difficulty of determining the underlying ecological principles [6]. Nonlinear differential equations are used to mathematically describe predator-prey interactions. However, it is typically difficult to find a suitable mathematical analysis, especially when using nonlinear terms.The Lotka–Volterra model, which is a system of nonlinear coupled first-order ordinary differential equations, has been deemed the basic model for describing predator-prey interactions [7, 8]. Two-predator one-prey models have the form of three species interactions, and thus, these systems are described by a system of three equations. Their dynamical behavior has been studied by some researchers [9–11].Seasonality is an important factor, which plays a vital role in describing the changes and fluctuations in ecological systems with predator-prey interactions [12–18]. Additionally, there are many ecological factors, such as hunting and climate, which have varied effects (positive and/or negative) on the dynamical behaviors of the species. In the literature, a number of studies have investigated the effect of seasonality on the dynamical behavior of predator-prey systems, but most of these studies have focused on the search of chaotic cases in predator-prey systems [12–14]. Several researchers [13, 15] have used impulsive differential equations to describe steep changes, where they studied the systems over a long period. However, in this paper, we will use the novel tool over describing steep changes for a long time or as a new situation is introduced into the system.Deterministic models have been widely used to describe predator-prey interactions and their dynamics. Deterministic models are useful, due to their ability to follow them through mathematical analysis, and they are an important mechanism for describing stable environments. However, random fluctuations appear in unstable environments, so deterministic models are difficult to describe these environments. In addition, the random noises are an important tool to conclude some unexpected dynamical behaviors of predator-prey interactions. Stochastic models play an important role for describing more realistic dynamical modeling of ecosystems. May [19] introduced an important contribution when he investigated stochastic differential equations for describing the limits of niche overlap in a randomly fluctuating environment. Recently, stochastic predator-prey models and their dynamics have been studied by some researchers [20–25].The study of the dynamical behavior of predator-prey interactions has been considered to be an important subject in applied mathematics and mathematical ecology, due to its universal existence and importance [26]. Stability is one of the main important dynamics of predator-prey systems, which is typically the first property considered when studying dynamical behavior.The persistence and extinction dynamics have also been discussed by many researchers [27–30], due to their importance. The analytical definitions of persistence and extinction are as follows: for a population p(t), if p(0) >0 and limt⟶∞pt>0, then p(t) persists, while if p (0) >0 and limt⟶∞pt=0, then p(t) becomes extinct. The geometric meaning of persistence is defined that each trajectory of a system of differential equations is bounded away from the coordinate axes, but the geometric meaning of extinction is that the trajectory of the system of differential equations touches the coordinate axes.The novelty of our work is on consideration of the deterministic and stochastic models taken in such a way; we are to get several results through our analysis. It should be noted that, we transfer the nonautonomous model to the autonomous model(s) by using a novel tool that approximates the model to particular cases.In this paper, we aim to investigate a cosinusoidal function in a Holling type I two-predator one-prey model, in order to study how sudden changes of the dynamics will effect on the dynamical behavior of the model. Investigating the cosinusoidal function and stochastic terms make our assumptions more realistic by concluding new cases of the model. We transfer the nonautonomous model to the autonomous model(s) by using a novel tool which approximates the model to particular cases.The paper is arranged as follows: we introduce inSection 1, the preface and methodology of the paper. In Section 2, we present the mathematical model of the two-predator one-prey system and the seasonality function. In Section 3, we introduce forced deterministic models by sudden changes, divided to two situations: bad and good. In Section 4, we present a mathematical analysis of the deterministic sudden changes. In Section 5, we study the equilibrium points and conduct a stability analysis of these situations. In Section 6, we introduce the stochastic model of the two-predator one-prey system and the numerical simulations. In Section 7, we summarize our conclusions.


### 1.2. The Methodology

I summarize the mechanism that is followed in this paper through Figure1.Figure 1
Summary of methodology.The methodology of arrays:Array 1: adding the stochastic termArray 2: adding the seasonality functionArray 3: using the approximation methodArray 4: theoretical analysisArray 5: numerical simulation


## 1.1. Preface

Theoretical ecology has motivated many mathematicians to discuss different ideas and models from a purely mathematical standpoint; see for example [1–5]. Mathematical modeling is a useful tool to determine how a process works and to predict what may follow [3]. Many problems taken into consideration in mathematical ecology seem simple, but are considered complicated problems due to the difficulty of determining the underlying ecological principles [6]. Nonlinear differential equations are used to mathematically describe predator-prey interactions. However, it is typically difficult to find a suitable mathematical analysis, especially when using nonlinear terms.The Lotka–Volterra model, which is a system of nonlinear coupled first-order ordinary differential equations, has been deemed the basic model for describing predator-prey interactions [7, 8]. Two-predator one-prey models have the form of three species interactions, and thus, these systems are described by a system of three equations. Their dynamical behavior has been studied by some researchers [9–11].Seasonality is an important factor, which plays a vital role in describing the changes and fluctuations in ecological systems with predator-prey interactions [12–18]. Additionally, there are many ecological factors, such as hunting and climate, which have varied effects (positive and/or negative) on the dynamical behaviors of the species. In the literature, a number of studies have investigated the effect of seasonality on the dynamical behavior of predator-prey systems, but most of these studies have focused on the search of chaotic cases in predator-prey systems [12–14]. Several researchers [13, 15] have used impulsive differential equations to describe steep changes, where they studied the systems over a long period. However, in this paper, we will use the novel tool over describing steep changes for a long time or as a new situation is introduced into the system.Deterministic models have been widely used to describe predator-prey interactions and their dynamics. Deterministic models are useful, due to their ability to follow them through mathematical analysis, and they are an important mechanism for describing stable environments. However, random fluctuations appear in unstable environments, so deterministic models are difficult to describe these environments. In addition, the random noises are an important tool to conclude some unexpected dynamical behaviors of predator-prey interactions. Stochastic models play an important role for describing more realistic dynamical modeling of ecosystems. May [19] introduced an important contribution when he investigated stochastic differential equations for describing the limits of niche overlap in a randomly fluctuating environment. Recently, stochastic predator-prey models and their dynamics have been studied by some researchers [20–25].The study of the dynamical behavior of predator-prey interactions has been considered to be an important subject in applied mathematics and mathematical ecology, due to its universal existence and importance [26]. Stability is one of the main important dynamics of predator-prey systems, which is typically the first property considered when studying dynamical behavior.The persistence and extinction dynamics have also been discussed by many researchers [27–30], due to their importance. The analytical definitions of persistence and extinction are as follows: for a population p(t), if p(0) >0 and limt⟶∞pt>0, then p(t) persists, while if p (0) >0 and limt⟶∞pt=0, then p(t) becomes extinct. The geometric meaning of persistence is defined that each trajectory of a system of differential equations is bounded away from the coordinate axes, but the geometric meaning of extinction is that the trajectory of the system of differential equations touches the coordinate axes.The novelty of our work is on consideration of the deterministic and stochastic models taken in such a way; we are to get several results through our analysis. It should be noted that, we transfer the nonautonomous model to the autonomous model(s) by using a novel tool that approximates the model to particular cases.In this paper, we aim to investigate a cosinusoidal function in a Holling type I two-predator one-prey model, in order to study how sudden changes of the dynamics will effect on the dynamical behavior of the model. Investigating the cosinusoidal function and stochastic terms make our assumptions more realistic by concluding new cases of the model. We transfer the nonautonomous model to the autonomous model(s) by using a novel tool which approximates the model to particular cases.The paper is arranged as follows: we introduce inSection 1, the preface and methodology of the paper. In Section 2, we present the mathematical model of the two-predator one-prey system and the seasonality function. In Section 3, we introduce forced deterministic models by sudden changes, divided to two situations: bad and good. In Section 4, we present a mathematical analysis of the deterministic sudden changes. In Section 5, we study the equilibrium points and conduct a stability analysis of these situations. In Section 6, we introduce the stochastic model of the two-predator one-prey system and the numerical simulations. In Section 7, we summarize our conclusions.


## 1.2. The Methodology

I summarize the mechanism that is followed in this paper through Figure1.Figure 1
Summary of methodology.The methodology of arrays:Array 1: adding the stochastic termArray 2: adding the seasonality functionArray 3: using the approximation methodArray 4: theoretical analysisArray 5: numerical simulation


## 2. Mathematical Model and Seasonality Function

### 2.1. Mathematical Model

We use a nondimensional system of the Holling type I two-predator one-prey model [31] as follows:(1)dxdt=x1−xk−αxy−βxz,dydt=−uy+e1αxy−e1αy2−c1yz,dzdt=−uz+e2βxz−e2βz2−c2yz,subject to initial conditions(2)x0=x0>0,y0=y0>0,z0=z0>0.The biological meaning of the variables and parameters is as follows:x: prey densityy: first predator densityz: second predator densityk: carrying capacity of the systemα and β: searching and capturing efficiency of predators y and zu and w: loss rates of predators y and ze1 and e2: birth rate of the predator for each prey consumedc1 and c2: interspecific competition between the predatorsThe parameters and initial conditions of the model (1) are supposed to be positive values.Theorem 1.
All the solutions of system (1) which initiate in R+3 for t≥ 0 are bounded.Proof.
according to the first equation of system (1), we prove that it is bounded as follows:(3)dxdt≤x1−xk.
The solution of equation (3) is xt=ket+kc/−1+et+kc, where c is the integration constant.
Then,0≤limt⟶∞supxt≤k∀t>0.
Then, we prove thatxt+yt+zt≤Q,∀t≥0.
LetRt=xt+yt+zt.
The derivative ofR with respect to time (t) is as follows:(4)dRdt=dxdt+dydt+dzdt,dRdt=1−xk−αy−βzx+−u+e1αx−e1αy−c1zy+−u+e2βx−e2βz−c2yz.
Since all the parameters are positive and the solutions initiating continue in the nonnegative quadrant inR+3, we can suppose the following:(5)dRdt=1−xkx+−u+e1αx−e1αyy+−u+e2βx−e2βzz.
We have that(6)maxR+x1−xk=k4.
By substituting in (5), it becomes as follows:(7)dRdt≤k4+−u+e1αx−e1αyy+−w+e2βx−e2βzz,(8)dRdt≤k4+−u+e1αx−e1αyy+−w+e2βx−e2βzz+Rt−Rt.
Equation (8) can be written as follows:(9)dRdt+Rt≤K4+x+−u+e1αx−e1αy+1y+−w+e2βx−e2βz+1z.
Sincext≤k, then(10)dRdt+Rt≤5k4+−u+e1α−e1αy+1y+−w+e2βx−e2βz+1z,but(11)maxR+−u+e1αk−e1αy+1y=−1+e1αk2−2e1αuk+u24e1α,maxR+−w+e2βk−e2βz+1z=−1+e2βk2−2e2βwk+w24e2β.
So, equation (10) becomes as follows:(12)dRdt+Rt≤Q,where(13)Q=145k+−1+e1αk2−2e1αuk+u2e1α+−1+e2βk2−2e2βwk+w2e2β.
Thus,Rt≤Q+ρe−t, where ρ is a constant of integration:(14)limsupt⟶∞Rt≤limt⟶∞Q+ρe−t.
Then,Rt≤Q.


### 2.2. Seasonality Function

Cosinusoidal and sinusoidal functions [12, 14, 16] are used for describing the effects of seasonality on the dynamical behavior of the model (1). The cosinusoidal function is as follows:(15)Ct=1+ϵcosμt,where the parameter ϵ indicates the seasonality degree (or strength seasonal degree) and the parameter μ represents the angular frequency of the fluctuations caused by impacts.


## 2.1. Mathematical Model

We use a nondimensional system of the Holling type I two-predator one-prey model [31] as follows:(1)dxdt=x1−xk−αxy−βxz,dydt=−uy+e1αxy−e1αy2−c1yz,dzdt=−uz+e2βxz−e2βz2−c2yz,subject to initial conditions(2)x0=x0>0,y0=y0>0,z0=z0>0.The biological meaning of the variables and parameters is as follows:x: prey densityy: first predator densityz: second predator densityk: carrying capacity of the systemα and β: searching and capturing efficiency of predators y and zu and w: loss rates of predators y and ze1 and e2: birth rate of the predator for each prey consumedc1 and c2: interspecific competition between the predatorsThe parameters and initial conditions of the model (1) are supposed to be positive values.Theorem 1.
All the solutions of system (1) which initiate in R+3 for t≥ 0 are bounded.Proof.
according to the first equation of system (1), we prove that it is bounded as follows:(3)dxdt≤x1−xk.
The solution of equation (3) is xt=ket+kc/−1+et+kc, where c is the integration constant.
Then,0≤limt⟶∞supxt≤k∀t>0.
Then, we prove thatxt+yt+zt≤Q,∀t≥0.
LetRt=xt+yt+zt.
The derivative ofR with respect to time (t) is as follows:(4)dRdt=dxdt+dydt+dzdt,dRdt=1−xk−αy−βzx+−u+e1αx−e1αy−c1zy+−u+e2βx−e2βz−c2yz.
Since all the parameters are positive and the solutions initiating continue in the nonnegative quadrant inR+3, we can suppose the following:(5)dRdt=1−xkx+−u+e1αx−e1αyy+−u+e2βx−e2βzz.
We have that(6)maxR+x1−xk=k4.
By substituting in (5), it becomes as follows:(7)dRdt≤k4+−u+e1αx−e1αyy+−w+e2βx−e2βzz,(8)dRdt≤k4+−u+e1αx−e1αyy+−w+e2βx−e2βzz+Rt−Rt.
Equation (8) can be written as follows:(9)dRdt+Rt≤K4+x+−u+e1αx−e1αy+1y+−w+e2βx−e2βz+1z.
Sincext≤k, then(10)dRdt+Rt≤5k4+−u+e1α−e1αy+1y+−w+e2βx−e2βz+1z,but(11)maxR+−u+e1αk−e1αy+1y=−1+e1αk2−2e1αuk+u24e1α,maxR+−w+e2βk−e2βz+1z=−1+e2βk2−2e2βwk+w24e2β.
So, equation (10) becomes as follows:(12)dRdt+Rt≤Q,where(13)Q=145k+−1+e1αk2−2e1αuk+u2e1α+−1+e2βk2−2e2βwk+w2e2β.
Thus,Rt≤Q+ρe−t, where ρ is a constant of integration:(14)limsupt⟶∞Rt≤limt⟶∞Q+ρe−t.
Then,Rt≤Q.


## 2.2. Seasonality Function

Cosinusoidal and sinusoidal functions [12, 14, 16] are used for describing the effects of seasonality on the dynamical behavior of the model (1). The cosinusoidal function is as follows:(15)Ct=1+ϵcosμt,where the parameter ϵ indicates the seasonality degree (or strength seasonal degree) and the parameter μ represents the angular frequency of the fluctuations caused by impacts.


## 3. Forced Deterministic Models by Sudden Changes

Events that happen unexpectedly (i.e., as the result of some environmental factors) on predator-prey interactions are called sudden changes. We apply the approximation method to describe such changes, in order to simplify the mathematical analysis of the model and make it biologically sensible. The approximation method has been applied for analyzing SIR models by some researchers [31, 32]. However, Alebraheem [16] has applied this technique to transfer a nonautonomous model containing seasonality terms to a autonomous model(s) by approximating the model to particular cases, in order to study the dynamical behavior of predator-prey systems.We apply the approximation method by taking the smallest and biggest values of the seasonality degree∈, where 0≤∈≤1. Hence, we approximate the cosinusoidal function (equation (15)) by the following two situations:(16)Ct≅Pt=0,bad situation,2,good situation.We interpret the “bad” and “good” situations as indicating surrounding circumstances are bad or good, respectively.We investigate the cosinusoidal function (equation (15)) in system (3) through three different cases, as follows.If sudden changes are forced for the whole system, we have the following:(17)dxdt=x1−xk−αxy−βxzPt,dydt=−uy+e1αxy−e1αy2−c1yzPt,dzdt=−wz+e2βxz−e2βz2−c2yzPt.If sudden changes are forced for the prey species through the growth rate of the prey, we have the following:(18)dxdt=Ptx1−xk−αxy−βxz,dydt=−uy+e1αxy−e1αy2−c1yz,dzdt=−wz+e2βxz−e2βz2−c2yz.If sudden changes are forced for both predator’ species through the birth rate of the predator for each prey consumed, we have the following:(19)dxdt=x1−xk−αxy−βxz,dydt=−uy+e1Ptαxy−e1αy2−c1yz,dzdt=−wz+e2Ptβxz−e2βPtz2−c2yz.


## 4. Mathematical Analysis of Deterministic Sudden Changes

In this section, we analyze the sudden changes’ effects on system (3) mathematically, so we substitute the values of Pt (equation (16)) through three cases (i.e., systems (17)–(19)) as follows:The first case: if sudden changes have an effect on the whole system.The bad situation: when we use Pt=0, system (17) becomes as follows:(20a)dxdt=0,(20b)dydt=0,(20c)dzdt=0.The solutions of equations (20a)–(20c) are as follows:(21)limt⟶∞xt=x0,limt⟶∞yt=y0,limt⟶∞zt=z0.We conclude from systems (20a)–(20c), the system will be set at the initial conditions.The good situation: when we use Pt=2, system (17) becomes as follows:(22)dxdt=2x1−xk−2αxy−2βxz,dydt=−2uy+2e1αxy−2e1αy2−2c1yz,dzdt=−2wz+2e2βxz−2e2βz2−2c2yz.The second case: if sudden changes have an effect on the prey species through the growth rate of the prey.The bad situation: when we use Pt=0, system (18) becomes as follows:(23a)dxdt=0−αxy−βxz,(23b)dydt=−uy+e1αxy−e1αy2−c1yz,(23c)dzdt=−wz+e2βxz−e2βz2−c2yz.The solution of equation (23a) becomes as follows.Becausey>0 and z>0, so we can reduce equation (23a) to become the following:(24)dxdt=−αx,xt=e−αt,then limt⟶∞xt=0.Sincey and z follow x, then limt⟶∞yt=0 and limt⟶∞zt=0.The good situation: when we use Pt=2, system (18) becomes as follows:(25)dxdt=2x1−xk−αxy−βxz,dydt=−uy+e1αxy−e1αxy2−c1yz,dzdt=−wz+e2βxz−e2βStz2−c2y.The third case: if sudden changes have an effect on both predator’ species through the birth rate of the predator for each prey consumed.The bad situation: when we use Pt=0, system (19) becomes as follows:(26a)dxdt=x1−xk−αxy−βxz,(26b)dydt=−uy+0−0−c1yz,(26c)dzdt=−wz+0−0−c2yz.For equations (26b) and (26c), we remove the terms −c1yz and −c2yz because they are negative terms and to simplify the mathematical analysis, so we have the following:(27)dydt=−uy−c1yz≅−uy,(28)dzdt=−wz−c2yz≅−wz.The solution of equation (27) is as follows:(29)yt=y0e−ut.Then, the solution leads to the following:(30)limt⟶∞yt=0.The solution of equation (28) is as follows:(31)zt=z0e−ut.The solution of this equation leads to the following:(32)limt⟶∞zt=0.Sincelimt⟶∞yt=0 and limt⟶∞zt=0, then equation (26a) becomes as follows:(33)dxdt=x1−xk.The solution of equation (26a) is as follows:(34)xt=ket+kc−1+et+kc,where c is the integration constant, then limt⟶∞xt=k.The good situation: when we use Pt=2, system (19) becomes as follows:(35)dxdt=x1−xk−αxy−βxz,dydt=−uy+2e1αxy−2e1αy2−c1yz,dzdt=−wz+2e2βxz−2e2βz2−c2yz.


## 5. Equilibrium Points and Stability Analysis

One of the main dynamical behaviors is stability. We find the positive equilibrium points to study the stability. To check the local stability, we compute the variational matrices corresponding to each equilibrium point and using the Routh–Hurwitz criterion for studying the stability. To check the global stability, we do that by constructing the Dulac function and Lyapunov function and using them to prove the global stability. We summarize the results of the equilibrium points of good situations when the sudden changes are forced through three cases in the following table (Table1).Table 1
Description of the equilibrium points, positive equilibrium points’ conditions, and stability condition(s) to each case for the good situation.
The equilibrium pointsPositive equilibrium points’ conditionsThe dynamical behaviorThe first caseE0=0,0,0No conditionsSaddle pointE1=k,0,0No conditionsGlobally stableE2=x,y,0=ku+e1/e1αk+e1,e1αk−u/e1α2k+e1α,0e1αk>uGlobally stableE3=x,0,z=kw+e2/e2βk+e2,0,e2βk−w/e2β2k+e2βe2βk>wGlobally stableE4=x,y,z= (kc1c2−kwc1α+kc2βu−ke1αβw−kue2αβ−ke1e2αβ/c1c2+kc2e1αβ+kc1e2αβ−e1e2αβ−ke1e2α2β−ke1e2αβ2,c1w+ke1αβw−ke2c1β−e2βu−ke2β2u+ke1e2αβ/−c1c2−kc2e1αβ−kc1e2αβ+e1e2αβ+ke1e2α2β+ke1e2αβ2,c2u−ke1αc2−e1αw−ke1α2w+ke2αβu+ke1e2αβ/−c1c2−kc2e1αβ−kc1e2αβ+e1e2αβ+ke1e2α2β+ke1e2αβ2)e1αk>ue2βk>wGlobally stableThe second caseE0=0,0,0No conditionsSaddle pointE1=k,0,0No conditionsGlobally stableE2=x,y,0=ku+2e1/e1αk+2e1,2e1αk−2u/e1α2k+2e1α,0e1αk>uGlobally stableE3=x,0,z=kw+2e2/wk+2e2,0,2e2βk−2w/e2β2k+2e2βe2βk>wGlobally stableE4=x,y,z= (2c1c2k−wc1αk+c2βuk−ke1αβw−kue2αβ−2ke1e2αβ/2c1c2+c2e1αβk+c1e2αβk−2e1e2αβ−ke1e2α2β−ke1e2αβ2,2c1w+ke1αβw−2ke2c1β−2e2βu−ke2β2u+2ke1e2αβ/−2c1c2−kc2e1αβ−kc1e2αβ+2e1e2αβ+ke1e2α2β+ke1e2αβ2,2c2u−2ke1αc2−2e1αw−ke1α2w+ke2αβu+2ke1e2αβ/−c1c2−kc2e1αβ−kc1e2αβ+e1e2αβ+ke1e2α2β+ke1e2αβ2)e1αk>ue2βk>wGlobally stableThe third caseE0=0,0,0No conditionsSaddle pointE1=k,0,0No conditionsGlobally stableE2=x,y,0=ku+2e1/2e1αk+2e1,2e1αk−u/2e1α2k+2e1α,02e1αk>uGlobally stableE3=x,0,z=kw+2e2/2e2βk+2e2,0,2e2βk−w/2e2β2k+2e2β2e2βk>wGlobally stableE4=x,y,z= (c1c2k−wc1αk+c2βuk−2ke1αβw−2kue2αβ−4ke1e2αβ/c1c2+2c2e1αβk+2c1e2αβk−4e1e2αβ−4ke1e2α2β−4ke1e2αβ2,c1w+2ke1αβw+2ke2c1β+2e2βu+2ke2β2u−2ke1e2αβ/−c1c2−2kc2e1αβ−2kc1e2αβ+4e1e2αβ+4ke1e2α2β+4ke1e2αβ2,c2u−2ke1αc2−2e1αw−2ke1α2w+2ke2αβu+4ke1e2αβ/−c1c2−2kc2e1αβ−2kc1e2αβ+4e1e2αβ+4ke1e2α2β+4ke1e2αβ2)e1αk>ue2βk>wGlobally stableWe present only the proof of the first case, and in the same manner, the proofs of the second and third cases will be followed, so the proofs of the second and third cases will be omitted.Theorem 2.
(i)
The trivial equilibrium pointE0=0,0,0 is a saddle point(ii)
The peak equilibrium pointE1=k,0,0 is locally asymptotically stable in x-direction, but it is locally asymptotically stable in y−z plane if it holds the conditions (38) and (39)Proof.
(i) we compute the variational matrix ofE0 which is given as follows:(36)M1=2000−2u000−2w.
Through the variational matrixM1, we see that the eigenvalues of y-direction and z-direction are negative, but the eigenvalue of x-direction is positive; this explains that the manifold is unstable along x-direction, but stable along y-direction and along z-direction. Then, the trivial equilibrium point E0 is the saddle point.
(ii) The variational matrix ofE1 is given as follows:(37)M2=−2−2kα−2kβ0−2u+2e1αk000−2w+2e2βk.
Through the variational matrixM2, we notice that the equilibrium point E1 is locally asymptotically stable, if the following conditions are satisfied:(38)u>e1αk,(39)w>e2βk.Theorem 3.
The peak equilibrium pointE1=k,0,0 is globally asymptotically stable under the following conditions:(40)u>4e1α2k,(41)w>4e2β2k.Proof.
Consider the following Lyapunov function aboutE1:(42)V1=x−k−klnxk+y2e1α+z2e2β.
whereV1 is a continuously differentiable real-valued function defined on R+3. Therefore, we have the following:(43)dV1dt=1−x¯xdxdt+12e1αdydt+12e2βdzdt,dV1dt=B1x−kxx1−xk−2αy−2βz+y2e1α−u+2e1αx−2e1αy−c1z+z2e2β−w+2e2βx−2e2βz−c2y,dV1dt=−1kx−k2−y2−z2−2αxy−2βxz−2e2βc1+2e1αc24e1e2αβyz−u−4e1α2k2e1αy−w−4e2β2k2e2βz.
If the conditions (40) and (41) are satisfied, then we obtain that dV1/dt<0 for any point in R+3.Theorem 4.
(i)
The equilibrium pointE2 is globally asymptotically stable in the interior of the positive quadrant of x−y plane(ii)
The equilibrium pointE3 is globally asymptotically stable in the interior of the positive quadrant of x−z plane
We prove part (i), and in the same manner, part (ii) can be proved.Proof.
letGx,y=1/xy, where G is a Dulac function. It is continuously differentiable in the positive quadrant of the x−y plane A=x,y|x>0,y>0:(44)N1x,y=2x1−xk−2αxy,N2x,y=−2uy+2e1αxy−2e1αy2.
Thus,ΔGN1,GN2=∂GN1/∂x+∂GN2/∂y=−2/yk−2e1α/x.
We find thatΔGN1,GN2<0 for all x>0 and y>0 in the positive quadrant of the x−y plane. By using Bendixson–Dulac criterion, there is no periodic solution in the interior of the positive quadrant of the x−y plane. E2 is globally asymptotically stable in the interior of the positive quadrant of the x−y plane.Theorem 5.
The persistence equilibrium pointE^=x^,y^,z^ of system (22) is globally asymptotically stable.Proof:.
we use the Lyapunov function to prove the global stability of positive equilibrium pointE^ as follows:(45)V=B1x−x^−x^lnxx^¯+B2y−y^−y^lnyy^+B3z−z^−z^lnzz^.
Equation (45) can be expressed as follows:(46)V=B1h1x,x¯+B2h2y,y¯+B3h3z,z¯,where h1x,x¯=x−x¯−x¯lnx/x¯.(47)h2y,y¯=y−y¯−y¯lnyy¯,h3z,z¯=z−z¯−z¯lnzz¯.
System (22) can be written as follows:(48)dxdt=xJx,y,z,dydt=yL1x,y,z,dzdt=zL2x,y,z,where(49)Jx,y,z=2−2xk−2αy−2βz,L1x,y,z=−2u+2e1αx−2e1αy−2c1z,L2x,y,z=−2w+2e2βx−2e2βz−2c2y.
Let(50)h′λ,λ¯=∂h∂λλ,λ¯=1−λ¯λ.
We compute the derivative ofV along the trajectories of system (22):(51)dVdt=B1h1′x,x¯dxdt+B2h2′y,y¯dydt+B3h3′z,z¯dzdt,which is(52)dVdt=B11−x¯xdxdt+B21−y¯ydydt+B31−z¯zdzdt,(53)dVdt=B1x−x¯xxJx,y,z+B2y−y¯yyL1x,y,z+B3z−z¯zzL2x,y,z.
Equation (53) can be expressed as follows:(54)=B1x−x¯xxJx,y,z−Jx¯,y¯,z¯+B2y−y¯yyL1x,y,z−L1x¯,y¯,z¯+B3z−z¯zzL2x,y,z−L2x¯,y¯,z¯.where Jx¯,y¯,z¯=0, L1x¯,y¯,z¯=0, and L2x¯,y¯,z¯=0, so we have(55)dVdt=B1x−x^2−2xk−2αy−2βz−2−2x^k−2αy^−2βz^+B2y−y^−2u+2e1αx−2e1αy−2c1z−−2u+2e1αx^−2e1αy^−2c1z^+B3z−z^−2w+2e2βx−2e2βz−2c2y−−2w+2e2βx^−2e2βz^−2c2y^.
Rearrange the terms of equation (55):(56)dVdt=B1kx−x^−2x−x^−2αy−y^−2βz−z^+B2y−y^2e1αx−x^−2e1αy−y¯−2c1z−z^+B3z−z^2e2βx−x^−2e2βz−z^−2c2y−y^,dVdt=B1k−2x−x^2−2αx−x^y−y^−2βx−x^z−z^+B22e1αy−y^x−x^−2e1αy−y^2−2c1y−y^z−z^+B32e2βz−z^x−x^−2e2βz−z^2−2c2z−z^y−y^.
By selectingB1=1, B2=1/e1, and B3=1/e2, so(57)dVdt=−2kx−x^2−2αy−y^2−2c1e1y−y^z−z^−2βz−z^2−2c2e2z−z^y−y^.
We find thatdV/dt is negative under no condition (i.e., no restrictions on parameters).
From Theorem5, we notice that the persistence dynamical behaviors of system (22) are globally stable.


## 6. Stochastic Model

In this section, we give numerical simulation to the stochastic version of our model. This consideration is due to the prevalence of randomness in almost all wild animal life, which makes the use of stochastic differential equations more realistic and efficient to describe some predictions of dynamical behaviors, see Figures2–5.Figure 2
Time series of system (59) without noiseσ=0.0.Figure 3
Time series of system (59) with low noise strengthσ=0.05.Figure 4
Time series of system (59) with medium noise strengthσ=0.5.Figure 5
Time series of system (59) with high noise strengthσ=1.4.The standard Ito^ stochastic differential equation is written as follows [24]:(58)dxt=Ft,xtdt+Gt,xtdWt,xt0=x0,where the first term represents the drift coefficient and the second term represents the random noise in the environment, which is sometimes called Gaussian white noise.We use a stochastic term in the deterministic model (2) as in reference [24], so we have the following model:(59)dxdt=x1−xk−αxy−βxz+σ1xdW1,dydt=−uy+e1αxy−e1αy2−c1yz+σ2ydW2,dzdt=−wz+e2βxz−e2βz2−c2yz+σ3zdW3,where σi, i=1,2,3 represent the strength of noise, and dWi, i=1,2,3 is a standard Wiener or Brownian motion processes.We have theoretically proven that the dynamical behavior of deterministic sudden changes is globally stable. In this section, we present the effects’ stochastic fluctuations on the dynamical behavior. The MATHEMATICA program was used to perform the numerical simulations. The values of parameters were selected to fulfil the positive values of a nontrivial equilibrium point, called the coexistence point (i.e., to satisfy conditions (38) and (39)) in deterministic models. In addition, the values of σi=σ were set as in [24], to represent three levels of noise strength; that is, low, medium, and high noise strengths. The parameters and initial conditions’ values were taken as follows:k=2.0, α=1.0, β=1.4, e1=0.6, e2=0.65, c1=0.07, c2=0.04, u=0.45, w=0.6, x0=0.6, y0=0.3, and z0=0.25.Figure2 represents the dynamical behavior of the model (59) without noise (i.e., σ=0), which gives the deterministic model. In Figure 2, the dynamical behavior of the species was stable coexistence, which corresponds with the theoretical analysis of the deterministic model. Figure 3 represents the dynamical behavior of the model (59) when the strength of the noise was low. Figure 3 shows that the dynamical behavior of the species was coexistence with smooth oscillations. However, with an increase in the strength of noise, such as the medium-noise situation shown in Figure 4, the dynamical behavior of the species was coexistence with sharp oscillations, which may lead to extinction. Figure 5 represents the dynamical behavior of the model (59) when the noise strength was high. The dynamical behavior of the species became more complex, and they tended to extinction. The random noises transfer the dynamical behavior from the equilibrium case to the oscillation case, which describes some unstable environments.We conclude that increasing the noise strength led to an increase of oscillation in the dynamical behavior, which can be interpreted biologically as increasing the probability of extinction, representing the worst-case scenario of dynamical behavior. This result corresponds with the numerical simulations. These results correspond well with the results of reference [24] with the difference being the mathematical model used, whereas increasing the noise strength led to an increase in oscillations in the dynamical behavior, finally leading to extinction when the noise strength was high.


## 7. Conclusion

We investigated the seasonality effects in a Holling type I two-predator one-prey model, which can more realistically describe the species of interaction more realistic. The nonautonomous models are transferred to autonomous models by approximating the model to particular cases representing sudden changes, so the situations are classified to bad and good situations, according to the surrounding circumstances. A mathematical analysis of sudden changes is introduced, and the equilibrium points and stability are discussed. We made the following conclusions.For the bad situations, we obtained the following outcomes:If sudden fluctuations have an effect on the whole system, then the system will remain at the initial conditionsIf sudden fluctuations have an effect on the prey species, then both predators’ species and the prey species will go extinctIf sudden fluctuations have an effect on both predator’ species, then the prey species will reach carrying capacity, while both of the predator’ species will go extinctThe equilibrium points of each case were obtained and found to be stableFor the good situations, we obtained the following outcomes:The one-prey two-predators system interacted through three different systems (22), (25), and (35) which represented the three casesWe obtained five positive equilibrium points, in each caseWe proved that the general dynamical behavior is globally stable, except for the trivial equilibrium point (which was a saddle point)The dynamical behavior in the case of good situations presented that the persistence dynamics is only a stable dynamical behaviorThrough numerical simulations, we presented effects of stochastic fluctuations on interactions, which showed that noise strength led to an increase in the oscillations in dynamical behavior and became more complex, finally leading to extinction when the noise strength was high.



---

*Source: 6611970-2021-03-03.xml*