Three-Species Predator-Prey Stochastic Delayed Model Driven by L\'{e}vy Jumps and with Cooperation Among Prey Species

Our study focuses on analyzing the behavior of a stochastic predator-prey model with a time delay and logistic growth of prey, influenced by L\'{e}vy noise. Initially, we establish the existence, uniqueness, and boundedness of a positive solution that spans globally. Subsequently, we explore the conditions under which extinction occurs, and identify adequate criteria for persistence. Finally, we validate our theoretical findings through numerical simulations, which also helps illustrate the dynamics of the stochastic delayed predator-prey model based on different criteria.


Introduction
Recently, the prey-predator problem has attracted the interest of many researchers in ecology [1][2][3]. The first prey-predator system, trying to describe the evolution of two subpopulations, the predator and the prey, was suggested in 1956 by Lotka and Volterra [4]. This basic model has played a crucial starting view in different studies of prey-predator dynamics. Understanding the different reactions between prey-predator components becomes then an important issue to control the extinction of predator or prey. For example, several models are used to better understand the behavior of viral infections, e.g., in lions, wildebeest and zebra [5], and fishing [6].
To better portray the growth of a population under constrained environments, Verhulst [7] introduced a special sort of growth equation, now known as the logistic equation. The logistic equation is very simple and very useful for one type of animal or one group of populations, as well as for more than one population. Without a predator, the growth rate of the prey population fulfills a logistic growth equation. Much research has been done on modeling by distinguishing the parameters for various species [8][9][10][11][12]. We should note that in some applications, such as the growth of population, the interaction of invulnerability, physiology, etc., means that the response growth rate of the species does not occur immediately, so it is necessary to include some kind of time delay [13]. Delays have a significant effect on the stability of the system. In some examples, the system is stable for a small time-delay value and turns into an unstable system for a higher delay [14]. The first to introduce a delay into the logistic differential equation was Hutchinson, with his paper from 1948 [15].
In the classical prey-predator model, it is assumed that the interaction is only within one prey and one predator. We can, however, take a model representing three interacting sub-populations, say one predator and two prey, where the prey being modeled do not have arXiv:2304.12860v1 [math.DS] 25 Mar 2023 competition among them, while mutual cooperation can be established between them against any predator. A work of 2018 by Kundu and Maitra proposes a prey-predator model with mutualism for two prey and one predator, also taking into consideration the time delays, as follows [16]: where x and y are the densities of the prey, the predator is represented by z, and the time delays τ 1 ≥ 0, τ 2 ≥ 0 and τ 3 ≥ 0 are incorporated in the growth component for each species, subject to the initial data with φ i ∈ C([−τ, 0], R + ) for i = 1, 2, 3. The description of the parameters of model (1) is given in Table 1. Table 1. The description of the parameters for model (1).

Parameters Description
The intrinsic growth rate for prey x r 2 The intrinsic growth rate for prey y K 1 The carrying capacity for x K 2 The carrying capacity for y α 1 The predation rate of prey x α 2 The predation rate of prey y β The rate of cooperation of preys x and y against predator z δ The predator death rate α 3 The rate of intra-species competition within the predators a 1 The transformation rate of predator to preys x a 2 The transformation rate of predator to preys y Many stochastic predator-prey problems have been deployed in order to describe the effect of white noise on predator-prey dynamics [3,[17][18][19][20]. In this context, in 2020, Rihan and Alsakaji studied a stochastic prey-predator model by introducing a Brownian perturbation to the three variables of the model and obtaining different conditions of extinction and persistence [17].
Lévy jumps present an important tool to model many real dynamic phenomena [21][22][23][24][25][26][27]. Because of the inherent stochastic properties of the prey-predator processes, it is natural to assume that the dynamic model may experience sudden strong perturbations in the predator-prey progression [28].
The authors of the previous delay model in [16] studied the effect of time delays on the stability of an equilibrium. Here we make model (1) more realistic and improve the results of [16] by considering that the predator and prey are exposed to uncertainties and randomness in the progress of a natural conflict. As a result, we consider here a stochastic system, which is more advantageous than the deterministic one because it takes into consideration not only the mean tendency but also the variance aspect surrounding it. In contrast to the deterministic model of [16], which always generates the same results for fixed initial values, our stochastic model may project different realities. In reality, here we add both random walks and Lévy jumps to (1) and, in contrast to [16], we study the questions of persistence and extinction, while in [16] their interests are different, being focused on bifurcation classification. Motivated by these facts and by the previous studies, here, we introduce jumps into model (1) as follows: where W i (t), i = 1, 2, 3, are standard Brownian motions defined on a complete probability space (Ω, F , (F t ) t≥0 , P) with the filtration (F t ) t≥0 and satisfying the usual conditions. We denote by x(t−), y(t−) and z(t−) the left limits of x(t), y(t) and z(t), respectively. Here N(dt, du) is a Poisson counting measure with the stationary compensator ν(du)dt,Ñ(dt, du) = N(dt, du) − ν(du)dt, where ν is defined on a measurable subset U of the nonnegative half-line with ν(U) < ∞ and the intensity of W i (t) is σ i , i = 1, 2, 3. The jump intensities are represented by q i (u), i = 1, 2, 3. Our model and work is new and original because it introduces both white and Lévy noise into the system (1), to show the effect on the sub-population of the prey under their cooperation and by taking their recruitment rate in logistic growth, and considers the effect of the stochastic perturbation of the predator. On the other hand, we also deal for the first time in the literature with the effect of the delays on the three considered species, describing the intra-production of the new prey and the generation of the predator after reaction with the prey.
The paper is organized as follows. In Section 2, we give some properties of the solution of model (3). The stochastic extinction of both the prey and the predator is established in Section 3. In Section 4, we prove the persistence of the prey and the extinction of the predator. The stochastic persistence of both the predator and the prey is obtained in Section 5 and Section 6 is devoted to some numerical simulations that illustrate the theoretical findings. We end with Section 7 of conclusions, final discussions, and some possible future directions of research.

Properties of the Solution
In this section, we show the well-posedness of our stochastic predator-prey delayed model driven by Lévy jumps. In reality, we prove the existence and uniqueness of a global positive solution (Theorem 1) as well as its boundedness (Theorem 2).

Existence and Uniqueness of a Global Positive Solution
Our first theorem guarantees the existence and uniqueness of a global positive solution for (3). Please note that the unique global solution (x(t), y(t), z(t)) predicted by Theorem 1 satisfies x(t) ≤ K 1 and y(t) ≤ K 2 since K 1 is the carrying capacity for x and K 2 is the carrying capacity for y. Theorem 1. Let δ > α 3 . For any given initial condition (2), the system (3) has a unique global solution (x(t), y(t), z(t)) ∈ R 3 + , t ≥ −τ, almost surely (a.s.).
Proof. The drift and the diffusion are locally Lipschitz. Thus, for any initial condition (2), there exists a unique local solution (x(t), y(t), z(t)) for t ∈ [−τ, τ e ), where τ e is the explosion time.
To demonstrate that this solution is global, we need to prove that τ e = ∞ a.s. First, we will prove that (x(t), y(t), z(t)) do not tend to infinity for a finite time. Let k 0 > 0 be a sufficiently large number, in such a way that (φ 1 , . Let us define, for each integer k ≥ k 0 , the stopping time We have established that τ k is an increasing number when k ↑ ∞. Let τ ∞ = lim k→∞ τ k , where τ ∞ ≤ τ e a.s. We need to show that τ ∞ = ∞, which means that τ e = ∞ and (x(t), y(t), z(t)) ∈ R 3 + a.s. Assume the opposite is verified, i.e., assume that τ ∞ < ∞ a.s. Then, there exist two constants 0 < < 1 and where Therefore, we have Integrating both sides of equation (4) between 0 and τ k ∧ T, we obtain For each h > 0, let us define Consequently, letting k → ∞, we have which represents a contradiction of the previous assumption. Therefore, τ ∞ = ∞ and the model has a unique global solution (x(t), y(t), z(t)) a.s.

Stochastic Boundedness
Theorem 1 shows that the solution of system (3) remains in the positive cone R 3 + . However, this nonexplosion property in a population dynamical system is often not good enough: the property of ultimate boundedness is more desired. Next, we prove stochastically ultimate boundedness.

Theorem 2. Under the conditions
the solution of system (3) given by Theorem 1 is stochastically ultimately bounded for any initial condition (2). where Then, Using conditions (5), we find that function f (x, y, z) has an upper bound. Denote f (x, y, z) and N = M + 1.
Since f (0, 0, 0) = 0, then N > 0. According to formula (6), Then, using Itô's formula, we obtain that Integrating both sides of (7) from 0 to t, and then taking expectations, we have where X = (x, y, z). This fact implies that For any > 0, let A = √ N √ . Using Chebyshev's inequality, we obtain that The proof is complete.
We remark that our condition (5) assures the boundedness of the solution, which is crucial in biological models. Roughly speaking, Theorem 2 tells us that we need to control the intensities σ i of the Brownian motions and the jump intensities q i .

Stochastic Extinction of the Prey and Predator
In this section, we show conditions under which the population becomes extinct with probability of one (Theorem 3). In the follow-up, we use the following notation: x(s)ds.

Theorem 3.
Let If max{c 1 , c 2 , c 3 } < 0, then the solution of system (3) satisfies Proof. Let us define Using Itô's formula, we have Observe that x → log(1 + x) − x is a nonpositive function. Since and so In addition, Using the same method, we prove that lim sup Now, consider Using Itô's formula, it follows that Then, Using the same technique we have used for x , we obtain that lim sup t→+∞ z ≤ 1 This completes the proof.

Stochastic Extinction of Predator
We now show conditions for which the population of predators becomes extinct with probability of one (Theorem 4).

Proof.
Using the first equation of (3), we have Let us define Then, using Itô's formula, we have Thus, Using the strong law of large numbers for martingales, we obtain With the same method, we prove that lim inf t→+∞ y ≥ r 2 − σ 2 2 2 1 − r 2 + 2r 2 /K 2 .

Now, consider
Using Itô's formula, we have that Then, We conclude that lim sup t→+∞ z ≤ 1 and the proof is complete.

Stochastic Persistence
This section is devoted to proving the persistence of the prey and predator populations.
Using now Itô's formula with Then, which proves the intended result.

Numerical Simulations
This section is devoted to illustrating our mathematical findings using numerical simulations. In the following examples, we apply the algorithm presented in [30] to solve system (3). In our simulations, the time period is in days. The different values for the parameters used in our numerical simulations are given in Table 2.
−0.008 −0.008 −0.008 Figure 1 shows the evolution of the two prey, x and y, and the predator z, in the case that both prey and predator vanish. This means the extinction of the populations, which is in agreement with Section 3. Figure 2 presents the dynamics of x, y and z in the case studied in Section 4, where the predator tends to zero and the prey remains strictly positive: the population of predators will be extinct while the population of prey persists. Figure 3 illustrates the behavior of x, y and z in a situation where both the predator and the preys remain strictly positive, i.e., all populations persist. This agrees with our theoretical result of Section 5. Figure 4 shows the impact of the transformation rates of predator to prey, a 1 and a 2 , based on the dynamics of the predator. We can remark that a 2 has more influence on the number of predators than the parameter a 1 . Figure 5 shows the impact of the carrying capacities K 1 and K 2 on the dynamics of the predator. We observe that K 2 has more influence on the number of predators than the parameter K 1 .
It is known that delays can affect the behavior of the studied population [31,32]. For this reason, now we study the effect of the delay on prey 1, prey 2, and the predator.     Figure 6 represents the dynamics of prey 1 for different values of τ 1 . We observe that when one increases the value of τ 1 from 0.5 to 2 days, that increases the amplitude of the prey 1 with a delay translation.   Figure 7 represents the dynamics of prey 2 for different values of τ 2 . We conclude that if we increase the value of τ 2 from 0.5 to 2 days, then the amplitude of the prey 2 increases with a delayed translation. In Figure 8, we show the behavior of the predator population for different values of τ 3 . We note that an increase of the value of τ 3 from 0.5 to 2 days results in an increase of the amplitude of the predators with a delayed response. Finally, in Figure 9 we illustrate the effect of the three delays. We observe that an increase on the values of τ 1 , τ 2 and τ 3 from 0.5 to 1 day implies an increase in the amplitude of all the populations with a translation.

Conclusions and Discussion
In this paper, we have proposed and analyzed a three-compartment model that depicts the interaction between two prey and one predator. The new stochastic predator-prey model incorporates Lévy noise and considers time delays of the two prey in the logistic functional production and in the transformation rate of predator to prey. Both white noise and Lévy jump perturbations are integrated into all model compartments. We have established the existence and uniqueness of a global positive solution and its boundedness, demonstrating the well-posedness of our stochastic predator-prey mathematical model. We have also provided sufficient conditions for the extinction of both prey and predator as well as for the persistence of the prey with the extinction of the predator. Additionally, we have presented a sufficient condition for the persistence in mean of both prey and the predator. Our theoretical findings are reinforced by different numerical results, demonstrating three possible scenarios: (i) the extinction of both prey and predator populations, (ii) the persistence of the prey and the extinction of the predator, and (iii) biological persistence, which is vital for maintaining the continuity of the reaction between prey and predator.
Several interesting research questions require further investigation in future studies. Specifically, our team intends to examine the spatial diffusion of a delayed prey-predator problem with Lévy jumps. Additionally, we plan to explore the bifurcation and chaotic behavior of the stochastic time-delayed prey-predator model. A comprehensive analysis is nontrivial and requires further inquiry, which we will undertake in another study. Furthermore, our analysis of the model under consideration can be expanded to include fractional-order derivative models, as described in [33,34], and we can explore alternative views of stochasticity, as demonstrated by the authors of [35]. We also note that the parameter values presented in Table 2 for numerical simulations of our theoretical results may not be applicable to real-world problems. Future research should strive to establish a connection with practical objectives and address the challenge of parameter estimation.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.