Deterministic and Stochastic Prey–Predator Model for Three Predators and a Single Prey

: In this paper, a deterministic prey–predator model is proposed and analyzed. The interaction between three predators and a single prey was investigated. The impact of harvesting on the three predators was studied, and we concluded that the dynamics of the population can be controlled by harvesting. Some sufﬁcient conditions were obtained to ensure the local and global stability of equilibrium points. The transcritical bifurcation was investigated using Sotomayor’s theorem. We performed a stochastic extension of the deterministic model to study the ﬂuctuation environmental factors. The existence of a unique global positive solution for the stochastic model was investigated. The exponential–mean–squared stability of the resulting stochastic differential equation model was examined, and it was found to be dependent on the harvesting effort. Theoretical results are illustrated using numerical simulations.


Introduction
Lotka and Volterra independently created two types of prey-predatory models known as the "Lotka-Volterra model" [1,2]. Since then, many scientists have modified and developed the Lotka-Volterra model to accurately describe the ecosystem. Numerous studies have examined the case of the presence of more than one predator [2][3][4][5][6][7][8][9][10][11]. Mukhopadhyay and Bhattacharyya [12] formulated a mathematical model of two predators living on a single biotic prey. They assumed that the predation function for the first predator follows the mass action kinetics, while the functional response for the second predator obeys the Holling type-II functional response. They also assumed that one of the predators is economically viable and undergoes harvesting at a rate proportional to its density. According to [13], in the northern Alaskan forest community, moose are the only large herbivores, constituting the primary prey for each of the three predators: black bears, gray wolves, and brown bears. Black bears have been known to attack and consume wolves if the opportunity presents itself. The main feature of this paper was to modify the interference of the predators in the system investigated in [12] by adding an extra predator y(t) where the first predator (black bear) preys on the second predator (gray wolves) in addition to the prey. The focus was on the harvesting rates and carrying capacity parameters of the model. The paper is organized as follows: The mathematical model is given in Section 2. The existence, uniqueness, non-negativity, and boundedness of the system are all verified in Section 3. Section 4 investigates the local and global stability of the system's equilibrium points. The stochastic extension of the deterministic model is conducted in Section 5. The numerical simulations presented in Section 6 are used to verify the theoretical results. Finally, in Section 7, the conclusions are presented.

Mathematical Model
In this paper, we considered a four-species prey-predator model with one prey and three predators as follows: where x(t) is the population size of the single prey species. We assumed that x(t) grows logistically in the absence of predators with intrinsic growth rate r and carrying capacity k. The first predator y(t) has the ability to consume both the prey and second predator z(t) with the Holling type I (linear) functional response. Let the interaction between the third predator w(t) and prey follow the Holling type II functional response. Assume β i (i = 1, 2, 3) denote the predation rates of the first, second, and third predators on the prey, respectively, and a is the half-saturation constant. Furthermore, let m i (i = 1, 2, 3) denote the efficiency of the first, second, and third predators in the presence of the prey. Moreover, δ represents the predation rate of the first predator on the second predator. We assumed that the ecological efficiency of the second predator's biomass z in the first predator's biomass y is unity. We also assumed that the predators economically undergo harvesting at a rate proportional to their density. The constants q i (i = 1, 2, 3) denote the catchability constants, while E represents the harvesting effort. The density of the first, second, and third predator populations decreases due to natural death at constant rates µ 1 , µ 2 , and µ 3 , respectively.

Existence and Uniqueness
In this section, we investigate the existence and uniqueness of the solutions of the prey-predator system (1) in the region Θ 1 × (0, T] where: for sufficiently large ϕ.
Proof. Define a mapping F(X) = (F 1 (X), F 2 (X), F 3 (X), F 4 (X)), in which: For any X,X ∈ Θ 1 , it follows from (1) that: where: Hence, F(X) satisfies the Lipschitz condition with respect to X. According to [14], as F(X) is locally Lipschitz, then there exists a unique local solution to the three-predator-one-prey system (1).

Non-Negativity and Boundedness
Considering the biological significance of the problem, we were only interested in nonnegative and bounded solutions. The prey-predator system (1) can be written as follows: with initial values Thus, the solution of the model (1), with non-negative initial conditions remains nonnegative. Furthermore, the solution satisfies the Lipschitz condition, as stated in Theorem 1. By Theorems 5 and 6 in [14], the solution of the prey-predator model (1) satisfies the non-negativity. The boundedness of the solutions of model (1) is given in the following theorem.
Theorem 2. The solutions of the prey-predator model (1) starting in R 4 + are uniformly bounded.
As a result, all the solutions of the prey-predator model (1) that start in R 4 + are uniformly bounded in the region: In the following, three critical parameters R 1 , R 2 , and R 3 , can be used to classify the dynamics of the prey-predator model (1). The threshold parameter R 1 is defined by . Using the next-generation method, one can obtain the basic reproduction number: .
One can note that the threshold parameter R 1 appears as a result of additional predator y(t) in the system considered in [12].

Equilibria and Stability
The prey-predator model (1) has the following seven equilibrium points: 1.
The predator free equilibrium point E 1 = (k, 0, 0, 0), which always exists; 3. The equilibrium point The equilibrium point E 5 = (x 5 , y 5 , z 5 , 0), where: The coexistence equilibrium point E 6 = (x 6 , y 6 , z 6 , w 6 ), where: One can note that the additional predator y(t) causes two new equilibrium points E 2 and E 5 to be obtained, which were not present in [12]. Now, the local stability of the system (1) is investigated. The Jacobian matrix is given as follows: The eigenvalues of J around the trivial point E 0 are r, −(µ 1 + q 1 E), −(µ 2 + q 2 E) and −(µ 3 + q 3 E); therefore, for all parameters, E 0 is a saddle with three-dimensional stable manifolds and a one-dimensional unstable manifold. The stability of the free predators' equilibrium point E 1 = (k, 0, 0, 0) is studied as follows: Proof. The Jacobian matrix of the model (1) around E 1 (J(E 1 )) is as follows: The eigenvalues of

Theorem 4. If
Proof. One can consider the positive-definite Lyapunov function as follows.
By calculating the time derivative of V 1 , one obtains: In accordance with Lyapunov-Sasalle's invariance principle, E 1 is globally stable when One can note that the global stability of E 1 depends on the parameters β 1 , µ 1 and q 1 of additional predator y(t), which were not present in [12]. The local bifurcation near the equilibrium point E 1 of the system (1) is now investigated using Sotomayor's theorem [15].

Theorem 5.
The prey-predator model (1) goes through a transcritical bifurcation regarding the bifurcation parameter q 1 around Proof. The Jacobian matrix of the prey-predator model (1) where ν 1 is any non-zero real number. Similarly, the eigenvector corresponding to J(E 1 ) T Q 2 = 0 is given by Q 2 = (0, τ 2 , 0, 0) T , where τ 2 is any non-zero number. Thus: In accordance with Sotomayor's theorem, the prey-predator model (1) has a transcritical bifurcation at q * 1 , which is equivalent to R 1 = 1. Therefore, the proof is complete.
The stability around E 2 = (x 2 , y 2 , 0, 0) is studied as follows: Proof. The eigenvalues of J(E 2 ) are: The eigenvalues λ 3 and λ 4 have negative real parts. Thus, if Proof. One can consider the positive-definite Lyapunov function as follows: By taking the time derivative of V 2 , one obtains, Thus, E 2 is globally stable if m 3 x 2 a(µ 3 +q 3 E) < 1 and The stability of the equilibrium point E 3 = (x 3 , 0, z 3 , 0) is investigated as follows: Proof. The eigenvalues of J(E 3 ) are: The eigenvalues λ 3 and λ 4 have negative real parts. Thus, if (µ 1 +q 1 E) < 1, then the equilibrium point E 3 is locally stable.
Proof. One can consider the positive-definite Lyapunov function as follows: By taking the time derivative of V 3 , one obtains, The stability of the equilibrium point E 4 = (x 4 , 0, 0, w 4 ) is studied as follows: Proof. The eigenvalues of J(E 4 ) are: The other eigenvalues are determined by: The stability around E 5 = (x 5 , y 5 , z 5 , 0) is studied as follows: The other roots are determined by where: when rδ > φ, then c 1 c 2 > c 3 . Hence, due to the Routh-Hurwitz criterion, all the eigenvalues of the Jacobian matrix J(E 5 ) around E 5 have a negative real part. Thus, the proof is complete.
Proof. One can consider the positive-definite Lyapunov function as follows.
By taking the time derivative of V 6 , one obtains, In accordance with Lyapunov-Sasalle's invariance principle, E 6 is globally stable if β 3 w 6 a(a+x 6 ) < r k and β 1 m 2 = β 2 m 1 .
In this section, we show that at the positive equilibrium point E 6 , a Hopf bifurcation arises, by taking the catchability constant q 3 , as a bifurcation parameter. The following lemma is presented first. Lemma 1. The characteristic Equation (5) has a pair of purely imaginary roots, and the remaining roots have negative real parts if and only if Suppose (5) has two eigenvalues, which have negative real parts, and two complex conjugates eigenvalues (call them (5) and separating the real and imaginary parts, one obtains: Following [16,17], substituting (6) into (7), differentiating with respect to q 3 , and utilizing m(q * 3 ) = 0 and n(q * 3 ) = 0, one obtains: Theorem 14. The system around the coexistence E 6 enters into Hopf bifurcation when q 3 passes q * 3 if the coefficients B j (q 3 )(j = 1, 2, 3, 4) at q 3 = q * 3 satisfy the following conditions: According to Theorem 14, there exists a Hopf bifurcation in the model (1), where the Hopf bifurcation is controlled by q 3 .

Stochastic Models
In this section, we perform a stochastic extension of the deterministic model (1) in two ways. Firstly, a randomly fluctuating driving force can be directly added to the deterministic model. Secondly, the catchability constants are replace by random parameters.

Stochastic Perturbations
Considering the effect of environmental noise, one can introduce a stochastic perturbation into the system (1); the stochastic prey-predator model takes the form: where W i (i = 1, 2, 3, 4) are independent standard Brownian motions with W i (0) = 0 and σ i > 0 denote the intensities of the white noise. In many applications, the solution of the Itô stochastic differential equation must preserve the positivity of the solutions [18][19][20]. According to Theorem 2.2 and Corollary 1 in [18], the solutions of (8) emanating from non-negative initial data (almost surely) remain non-negative if they exist. In the next theorem, another approach according to [17] to prove the existence and uniqueness of a positive global solution of model (8) is given.
Theorem 15. For any given initial value x 0 , y 0 , z 0 , w 0 ∈ R 4 + , there exists a unique solution x(t), y(t), z(t), w(t) of the system (8) on t ≥ 0, and the global positive solution will remain in R 4 + with probability one.
Here, we allowed stochastic perturbations of x, y, z, w around the free predators' equilibrium point E 1 . The linearized stochastic system can be written as: where: U(t) = (u 1 (t), u 2 (t), u 3 (t), u 4 (t)) T . One can note that the free predators' equilibrium E 1 of the system (1) corresponds to the trivial solution of the system (10). Following [17,20], let B be the set defined as B = [(t ≥ t 0 ) × R n , t 0 ∈ R + ] and V ∈ C 0 2 (B) be a twice-differential function with respect to U and a continuous function with respect to t. Now, we require the following theorem to prove the asymptotically mean-squared stability of the trivial solution of (10).
Following [21,26,27], the Lyapunov operator LV(t, U) associated with (12) is defined as: Theorem 17. The trivial solution of (10) is asymptotically mean-squared stable if: Proof. Consider the following Lyapunov function: The first condition of Theorem 16 holds for the Lyapunov function defined in (13) with p = 2. Now, Lyapunov operator LV 7 (t, U) becomes: and this leads to LV 7 (t, U) ≤ −K 3 U 2 , where: One can note that the conditions of Theorem 17 indicate that the exponential-meansquared stability of the system (10) depends on the harvesting effort.

Random Harvesting
Here, we studied the effect of random harvesting on the three predators. The stochastic extension of (1) is as follows: The catchability parameters q 1 , q 2 , and q 3 were perturbed by independent Gaussian white noise terms ζ 1 , ζ 2 , and ζ 3 in the system (14) because, usually in the prey-predator system, harvesting is performed randomly, where ζ i , i = 1, 2, 3 are independent Gaussian white noises satisfying: is the Dirac delta function; δ ij is the Kronecker delta; . is the expectation.

Numerical Simulations
In this part, the numerical simulations are compared with the previous theoretical analysis. The numerical simulation was conducted using the following parameters The effect of catchability constants can be shown by drawing the bifurcation diagram regarding q 1 as a bifurcation parameter. The transcritical bifurcation value is centered at q * 1 = 0.2 as indicated in Figure 1. Note that the bifurcations that are presented in Theorem 5 are illustrated because q * 1 = 0.2 is equivalent to R 1 = m 1 k µ 1 +q 1 E = 1. One can draw the bifurcation diagram regarding q 3 to indicate the effect of the harvesting. The supercritical Hopf bifurcation value is centered at q * 3 = 0.254669 as indicated in Figure 2. It can also be noted that for q 3 > 0.254669, the prey-predator model (1) is locally stable as indicated in Figure 3, while for q 3 < 0.254669, the system goes through the limit cycle behavior. One can find that all the conditions of Theorem 14 hold as Φ 1 (0.254669) = 0, Φ 2 (0.254669) = 0 and dΦ 1 (q 3 ) dq 3 | q 3 =0.254669 = 0. This confirms the existence of a Hopf bifurcation at q * 3 = 0.254669. As a result, the harvesting parameter q 3 can break the oscillating behavior of the deterministic system (1) and drive it to the required state. In the same way, the bifurcation of the system can be studied using the parameter q 2 , as shown in Figure 4. x,y,z,w q 3 =0.4 x y z w  To better understand the effect of the caring capacity k, one can draw the bifurcation diagram with respect to k. It can be seen that the supercritical Hopf bifurcation value is localized at k = 0.45 as shown in Figure 5. The supercritical Hopf bifurcation value is centered at k = 0.45, as indicated in Figure 5. When k > 0.45, the prey-predator model (1) goes through limit cycle oscillation, as indicated in Figures 5 and 6. For k < 0.45, E 4 = (0.075, 0, 0, 0.15 − 0.01125/k) is locally stable, as indicated in Figure 6. It can also be noted that the conditions of local stability that were established in Theorem 10 were verified because when k = 0.4, one has R 3 = 2.8571 < 1 + (am 3 ) ((µ 3 +q 3 E)(a+k)) = 3.1429. To give some numerical findings for the prey predator system (8), one can use the Milstein method mentioned in [31,32]. The prey-predator system (8) reduces to the following discrete system.
where h is a positive time increment and ij , (i = 1, 2, 3, 4) are independent random Gaussian variables N(0, 1). Figure 7 represents the dynamical behavior of the model (8) when the noise strength is the law (σ i = 0.05). One can note that for the given parameters, the strength of environmental noise is very close to zero, and the system behaves as a deterministic model. Following [33], one can note that in the deterministic case, if R 0 < 1, then the prey-predator system (1) has a predators' free equilibrium point E 1 = (k, 0, 0, 0). In the stochastic model (14), if one gradually increases the intensities of fluctuation and keeps the remaining parameters unchanged, the fluctuations around E 1 become larger, as seen for the values of σ i = 0.2 and 0.9 shown in Figure 7. The black line in Figure 7 represents the prey when (σ i = 0). From Figure 8, it is seen that increasing the catchability constants has a stabilizing effect on the stochastic model (14).   (14) with (q 1 = q 2 = q 3 = 0.1) and (q 1 = q 2 = q 3 = 0.5).

Conclusions
In this paper, a mathematical prey-predator model was proposed and analyzed. The interference of the predators in the system investigated in [12] was modified by adding an extra predator y(t) where the first predator preys on the second predator in addition to the prey. The interaction between the three predators and single prey was studied. The impact of harvesting on the first and the second predator was investigated. Sufficient conditions were obtained to ensure local stability. It was concluded that the dynamics of the population can be controlled by harvesting. The harvesting rates of the three predator species played an important role in controlling the local and global dynamics of the preypredator system. They can break the oscillating behavior of the deterministic system and drive it to the required state. To investigate the effect of environmental noise, we performed a stochastic extension of the deterministic model to study the fluctuation of the ecological factors. The existence of a unique global positive solution for the stochastic model was investigated. We used stochastic perturbation around the free predators' equilibrium point. Constructing an appropriate Lyapunov function and applying Itô's formula, we note that the deterministic model was robust with respect to stochastic perturbation. The criterion of stochastic stability depends on the intensities of noise σ i , i = 1, 2, 3, 4. The exponential-mean-squared stability of the resulting stochastic differential equation model was examined, and it was found to be dependent on the harvesting effort.