Front. Appl. Math. Stat.Frontiers in Applied Mathematics and StatisticsFront. Appl. Math. Stat.2297-4687Frontiers Media S.A.10.3389/fams.2025.1523276Applied Mathematics and StatisticsHypothesis and TheoryDynamic complexity of a delayed spatiotemporal predator-prey modelHafdaneMohamed1*BabaNossaiba1El FoutayeniYoussef2AchtaichNaceur11Analysis, Modeling and Simulation Laboratory, Hassan II University, Casablanca, Morocco2Cadi Ayyad University, École Nationale des Sciences Appliquées, Marrakech, Morocco
Edited by: Preeti Dubey, University of Washington, United States
Reviewed by: Appanah Rao Appadu, Nelson Mandela University, South Africa
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
This study investigates a delayed spatiotemporal predator-prey model that incorporates key ecological mechanisms, including the Allee effect, fear-induced prey behavior, Holling type II predation with cooperative hunting, toxicity with delayed effects, and both nonlinear (for prey) and linear (for predators) fishing pressures. Using tools from the theory of partial differential equations, stability analysis, and Hopf bifurcation theory, we derive the conditions under which stable coexistence or instability emerges. Our results reveal that system stability is maintained below a critical delay threshold, beyond which oscillatory dynamics arise. In the spatial domain, diffusion can either stabilize populations or lead to heterogeneous patterns such as Turing structures and predator-prey segregation, particularly when diffusion is low and delays are significant. Numerical simulations support and illustrate the analytical findings, showing a variety of dynamic behaviors consistent with observed ecological patterns. This work highlights how the interplay between ecological processes, time delays, and spatial effects governs predator-prey dynamics and offers insights relevant to ecosystem management.
The study of predator-prey interactions has long been fundamental in ecological modeling. Understanding how species coexist and how their populations fluctuate over time and space is essential for both theoretical ecology and natural resource management [1–6]. Natural ecosystems are often subject to various biological and environmental factors that influence species dynamics. In particular, resource competition [7], predation [8], fear effects [9, 10], and anthropogenic disturbances such as fishing play a major role in shaping population dynamics. Numerous studies have explored these dynamics using ordinary and partial differential equations [6, 11–14], allowing for a comprehensive representation of these complex interactions. For example, in Chakraborty et al. [15], the authors investigated a predator-prey model incorporating prey refuge and additional food sources for predators. Supplemental feeding is an effective strategy in integrated pest management and conservation programs. Their findings show that a high level of prey refuge can negate the benefits of supplemental food, making prey control difficult. However, when prey refuge is limited and predators have access to an optimal level of additional food, prey (pest) populations can be effectively regulated.
In this study, we propose a spatiotemporal predator-prey model that incorporates several key ecological factors. First, we integrate the Allee effect [16], which reflects the difficulty prey face in reproducing at low densities, when cooperation among individuals becomes insufficient. The Allee effect, represented by F(b,N)=aNb+N, where N is the prey density and b is the density threshold below which reproduction is strongly reduced, is crucial for many animal and plant populations. This effect is particularly significant in situations where reproduction depends directly on the presence of a sufficient number of individuals [46–50].
Furthermore, the fear effect alters prey behavior in response to predator presence [14, 17, 18]. This mechanism is modeled by a scaling term R(c,P)=11+cP, which reduces prey activity and limits their access to resources, significantly affecting their reproduction rate. Moreover, the fear effect can amplify the Allee effect under certain conditions, thus altering the system dynamics and reducing the prey population in a more complex manner than a simple density effect. Recent studies also emphasize the combined roles of fear of predation, supplemental food availability, and selective predation in shaping ecological stability [19]. Fear can temporarily protect prey by reducing predator encounters, but simultaneously limits prey access to essential resources. Supplemental food supports predator persistence, while selective predation, where predators avoid infected prey, can significantly alter disease transmission and population dynamics in both temporal and spatial settings.
Another fundamental aspect of the model is the incorporation of a Holling type II functional response combined with cooperative predation. The Holling type II functional response describes a saturable consumption rate, reflecting a predator's physiological limitation. However, certain predator species hunt cooperatively, which enhances their efficiency when in groups. This cooperative hunting behavior is integrated into the model through a term that modifies the capture probability as a function of predator density.
Introducing toxicity [14] and its delayed effects adds another critical dimension to the model. Prey accumulates toxins that do not act immediately but instead cause delayed mortality, introducing a memory effect into the system. In contrast, predators experience the effects of toxicity instantaneously when consuming intoxicated prey. This temporal lag in toxin effects can induce instabilities and influence coexistence patterns and population regulation. In Shukla et al. [20], the authors studied the spatial dynamics of a nutrient–phytoplankton system under toxic effects and found that toxicity can lead to spatially inhomogeneous distributions, producing diverse patterns such as stripes, spots, or mixtures thereof. Their findings revealed that certain levels of toxicity could drive spatio-temporal oscillations, further emphasizing the potential for toxicity to generate complex dynamical behaviors in ecological systems.
Incorporating spatial diffusion into our predator-prey models is essential for realistically capturing the movement and distribution of biological populations. Diffusion terms model the natural tendency of individuals to migrate across space, which can interact with local dynamics to generate spatial heterogeneity. For instance, Shukla et al. [20] investigated the effects of cross-diffusion in an algal bloom model and demonstrated that spatial interactions could lead to the formation of complex patterns depending on environmental parameters. Their analysis highlights how spatial diffusion can either stabilize or destabilize homogeneous equilibria, depending on system conditions, and emphasizes the necessity of including such mechanisms in ecological models.
Finally, we account for resource exploitation through fishing. Prey are subject to a nonlinear Michaelis–Menten-type fishing pressure [21–23], which reflects capture saturation at high abundance levels. Predators, on the other hand, experience linear fishing pressure [24–26], corresponding to an extraction rate proportional to their density. These two forms of harvesting allow us to assess the effects of fisheries management on population stability and resilience to anthropogenic pressures [27–30].
The main objective of this study is to analyze the stability of the coexistence equilibrium of both species and explore the spatiotemporal dynamics of the system using advanced mathematical tools. We rely on eigenvalue analysis to identify equilibrium stability conditions, Hopf bifurcation theory [31–33] to examine the emergence of oscillatory behaviors due to time delays, and parabolic and elliptic partial differential equations to study spatial patterns emerging in population distribution. Numerical simulations are also conducted to illustrate the effects of various parameters on species coexistence and spatial population structuring.
This study makes a significant contribution to predator-prey system modeling by simultaneously incorporating multiple ecological and bioeconomic factors. The results will provide a deeper understanding of how spatial diffusion [34–36], delayed toxicity [12, 37–44]), and predator cooperation influence population stability and persistence. Furthermore, the inclusion of fishing pressure in the model makes our study particularly relevant to marine resource management policies and the conservation of exploited species. By combining theoretical analysis with numerical simulations, we aim to characterize the system's dynamic transitions and identify the conditions that promote species coexistence.
The predator-prey system studied is described by the following partial differential equations:
Ω is a smooth and bounded open region in ℝN (where N ≥ 1), and ∂Ω represents its smooth boundary. The operator Δ corresponds to the Laplacian in ℝN. The unit outward normal vector on ∂Ω is denoted as S. The usage of homogeneous Neumann boundary conditions implies that the population under consideration cannot migrate or move across the boundaries of the given domain. Furthermore, we make the assumptions that N0(x, t) and P0(x, t) belong to the space C([−τ, 0], X), where X is defined by:
X={v,w∈W2,2(Ω):∂v(t,x)∂ν=∂w(t,x)∂ν=0,x∈∂Ω}
and the inner product 〈·, ·〉. For the sake of convenience, we limit our study to the one-dimensional spatial domain Ω = (0, lπ) throughout this paper. The parameters used in the model are defined in Table 1.
The meaning of bioeconomical parameters.
Parameter
Meaning
a
Maximum filtering per individual of the population
b
Intensity of the Allee effect
d
Mortality rate of prey
e
Intensity of intraspecific competition
c
Fear induced by predator population
f
The attack rate per predator to prey
g
Predator cooperation in hunting
h
The predators handing time of a prey
m
Death rate for predators
k
the conversion efficiency
τ
The delay for toxins to affect prey
m1
Impact of fishing effort on capture saturation
m2
Impact of fish abundance on capture saturation
θ1, θ2
Toxicity effect
d1, d2
Diffusive coefficients
q1, q2
Capturability coefficients
E1, E2
The level of fishing effort employed to exploit the targeted species
The organization of this article is as follows. In the Section 1, we analyze the temporal model without considering the spatial dimension to better understand its intrinsic dynamics. We study the existence, positivity, and boundedness of the solutions, followed by an analysis of local stability. The Section 2 is dedicated to the study of the spatiotemporal model, where we establish the existence and boundedness of solutions, derive an a priori estimate for positive solutions, and formulate conditions for the non-existence of non-constant positive solutions. A stability analysis of the system is also conducted. In the Section 3, we examine the direction of the Hopf bifurcation and analyze the stability of the emerging periodic solutions. Finally, the Section 4 is devoted to numerical simulations, which illustrate the theoretical results on the effect of delay on stability and examine the influence of diffusion coefficients on the model's dynamics. These analyses are followed by a discussion aimed at interpreting the obtained results and highlighting their biological implications.
With specified initial conditions x(θ) = φ(θ) > 0 and y(θ) = ψ(θ) > 0 for all θ ∈ [−τ, 0], where φ and ψ are continuous functions.
2.1 Existence, positivity, and boundedness of the solution2.1.1 Positivity
Theorem 1. The set {(N, P) ∈ ℝ2 : N ≥ 0, P ≥ 0} is positively invariant for Equation 2.
Proof. Note that the planes N = 0 and P = 0 are invariant. Indeed,
dN(t)dt|N=0=0and dP(t)dt|P=0=0.
This means that if solutions reach N = 0 or P = 0, they do not become negative. Thus, starting with strictly positive initial conditions N(0) > 0 and P(0) > 0, solutions remain in the positive domain for all t > 0.
In conclusion, the set {(N, P) ∈ ℝ2:N ≥ 0, P ≥ 0} is positively invariant for Equation 2.
2.1.2 Boundedness
Theorem 2. All solutions of system (Equation 2), with positive initial values, are bounded.
Proof. Examine the inequality below:
dNdt≤N(a-d-eN).
This implies that N is bounded above by a-de.
Now define U=N+1kP where k is a positive constant, and let Φ be a positive constant. We get:
Theorem 3. Equation 2 has only one possible solution.
Proof. The function g : ℝ2 × ℝ2 → ℝ2 is well-defined and continuous. Moreover, for each gi (i = 1, 2), the partial derivatives exist and are assumed to be continuous and bounded. As a result, g satisfies the Lipschitz condition locally, ensuring the existence and uniqueness of a local solution X(t) according to the Cauchy-Lipschitz theorem for functional differential equations with delay [11].
2.2 Local stability
In this section, we initially identify and characterize all the equilibrium points of Equation 2, followed by an analysis of their local stability.
2.2.1 Equilibrium points
We concentrate exclusively on the dynamic study of internal equilibrium points. To find these, we need to solve the following system:
(H2)β02-γ02⩾0, (α02-2β0-θ22)2-4(β02-γ02)⩾0 and α02-2β0-(θ1N*)2<0
(H3)β02-γ02⩾0, (α02-2β0-(θ1N*)2)2-4(β02-γ02)⩽0 or α02-2β0-(θ1N*)2>0
(H4)β02-γ02⩾0, (α02-2β0-(θ1N*)2)2-4(β02-γ02)⩾0 and α02-2β0-(θ1N*)2>0
Theorem 5. Suppose that the conditions (H0) hold. If one of the conditions (H3) or (H4) is satisfied, then Equation 2 is locally asymptotically stable for all τ ≥ 0.
Next, we demonstrate under what conditions Equation 2 experiences a Hopf bifurcation by considering the delay τ as a bifurcation parameter. The necessary condition for a shift in stability of the interior equilibrium E* is that the characteristic Equation 6 possesses purely imaginary roots. Thus, to derive the stability criterion, we substitute τ=τ^ and ω=ω^ into Equations 7 and 8, and by solving these equations for cos(ω^τ^) or sin(ω^τ^), we obtain:
where n ∈ N. The transversality condition for Hopf bifurcation at τ=τ^ is [dμdτ]τ=τ^>0. Let λ = μ + iω be the root of Equation 6 satisfying μ(τ^)=0 and ω(τ^)=ω^. Differentiating both sides of this equation with respect to τ, we get
The transversality condition [dμdτ]τ=τ^>0 for the occurrence of Hopf bifurcation at τ=τ^ is properly satisfied as long as Q3Q1 − Q4Q2 > 0. Consequently, we obtain the following result:
Theorem 6. Suppose that condition (H0) holds. If (H1) or (H2) is satisfied, then Equation 2 is locally asymptotically stable for τ<τ^ and becomes unstable when τ>τ^. Furthermore, when τ=τ^, Equation 2 undergoes a Hopf bifurcation at (N*, P*) provided that M3M1 − M4M2 > 0.
3 Spatiotemporal model study3.1 Existence and boundedness of the solution
Theorem 7. For Equation 1, we have the following results:
If N0(x, t) ≥ 0 and P0(x, t) ≥ 0, then Equation 1 has a unique positive solution (N(x, t), P(x, t)) for x ∈ Ω and t ∈ (0, ∞).
If (N(x, t), P(x, t)) is a solution of Equation 1, then
lim supt→+∞N(x,t)⩽a-de.
Moreover, there exist constants C1 and C3 such that
Where N0¯=supΩ¯N0(x,t) and P0¯=supΩ¯P0(x,t), ∀t ∈ [−τ, 0]. Let (N~,P~) be the unique solution of the Equation 11. Then (N, P) = (0, 0) and (N¯,P¯)=(N~,P~), are respectively the lower and upper solutions of the system 1. Thus, the Equation 1 has a unique globally defined solution (N(x, t), P(x, t)), which satisfies
0⩽N(x,t)⩽N~(t),0⩽P(x,t)⩽P~(t)
The strong maximum principle ensures that N(x, t) > 0 and P(x, t) > 0.
Proof. Suppose that (N(x), P(x)) is a solution of Equation 13 satisfying N(x), P(x) ⩾ 0. According to the strong maximum principle, we have N(x) > 0 and P(x) > 0. From Lemma 1, we obtain that N(x, t) ⩽ (a − d)/(e + θ1). Multiplying the first equation of 13 by k and adding it to the second equation of 13, we obtain
P⩽kd2(a-d+md1d2)24d2m(e+θ1).3.3 Non-existence of non-constant positive solutions
Theorem 9. For any fixed a, b, c, d, e, E1, E2, f, g, h, k, m, m1, m2, q1, q2, θ1 and θ2, there exists a positive constant d* such that if min{d1,d2}>d*, then Equation 13 has no nonconstant solutions.
Proof. Let (N(x), P(x)) be non-negative solution of Equation 13, We define N̄=|Ω|-1∫ΩN(x)dx and P̄=|Ω|-1∫ΩP(x)dx. It's clear that ∫N-N̄dx=0 and ∫P-P̄dx=0. To facilitate the discussion, let χ(N, P) = (f + gP)NP/(1 + h(f + gP)N). According to the mean value theorem for bivariate functions, we have:
If no delay is present, the characteristic equation will take the following form:
λ2+(αn+θ1N*)λ+βn+γn=0,
We make the following assumption:
(H5)αn+θ1N*>0, βn + γn > 0, for n ∈ ℕ0,
(H6), α0+θ1N*>0,αk+θ1N*<0, (or βk + γk < 0), for some k ∈ ℕ.
Theorem 10. For the Equation 1, suppose that τ = 0. The point (N*, P*) is locally asymptotically stable under (H5) and is Turing unstable under (H6).
Proof. If (H5) holds, we can determine that the characteristic roots of Equation 15 all have negative real parts. Hence, (N*, P*) is locally asymptotically stable. If (H6) holds, then the characteristic roots for k ∈ ℕ have at least one positive real part, but with n = 0, they all have negative real parts. This implies that (N*, P*) is Turing unstable.
Therefore, Re(dλdτ)|τ=τnj,+>0, and Re(dλdτ)|τ=τnj,-<0.
Let τ*=min{τn0∣n∈𝕎1∪𝕎2}. We have the following theorem:
Theorem 11. Suppose that (H5) is satisfied. Then, the following results hold:
The positive equilibrium (N*, P*) of the Equation 2 is asymptotically stable for τ ∈ [0, τ*).
The Equation 2 undergoes a Hopf bifurcation at the positive equilibrium (N*, P*) when τ=τnj,± for n ∈ 𝕎1 ∪ 𝕎2 and j ∈ ℕ0.
4 Hopf bifurcation
In this section, our goal is to obtain the normal form of Hopf bifurcation at the interior equilibrium. Let N̄(x,t)=N(x,τt)-N* and P̄(x,t)=P(x,τt)-P*. In this context, we've omitted the bar for simplicity. Thus, the resulting system is as follows
Respectively, for φ=(φ1,φ2)T∈C1. We know that Λn:={iωnτ~,-iωnτ~} are characteristic roots of
dz(t)dt=-τ~Dn2l2z(t)+Lτ~(zt)
The application of the Riesz representation theorem allows us to establish the existence of a 2 × 2 matrix function ηn(s,τ~),(-1≤s≤0), whose elements are of bounded variation functions such that
for φ∈C([-1,0],ℝ2),ψ∈C([0,1],ℝ2).A(τ~) has a pair of simple purely imaginary eigenvalues ±iωnτ~, and they are also eigenvalues of A*. Define p1(θ)=(1,ζ)Teiωnτ~s(s∈[-1,0]),q1(r)=(1,ϑ)e-iωnτ~r(r∈[0,1]), where
Theorem 12. The subsequent conclusions apply to any critical value τnj,+ (or τnj,-).
μ2 dictates the directions of the Hopf bifurcation: if μ2 > 0 (or μ2 < 0), the Hopf bifurcation is forward (or backward), indicating that the resulting periodic solutions exist for μ > 0 (or μ < 0).
β2 determines the stability of the bifurcating periodic solutions on the center manifold: if β2 < 0 (or β2 > 0), then the bifurcating periodic solutions are orbitally asymptotically stable (or unstable).
T2 dictates the period of bifurcating periodic solutions: if T2 > 0 (or T2 < 0), then the period increases (or decreases).
5 Simulation
In this part, we provide numerical results that explore the influence of delay and diffusion parameters on the dynamics of our model. To achieve this goal, we employ the parameter values provided in the Table 2.
The values of bioeconomical parameters.
Parameter
a
b
c
d
e
f
g
h
q1
Value
2.2
0.1
0.8
0.8
0.1
0.6
0.04
0.6
0.07
Parameter
m
k
m1
m2
E1
E2
q2
θ1
θ2
Value
0.2
0.25
0.5
0.4
0.8
0.9
0.05
0.08
0.05
A straightforward calculation confirms that the Equation 1 has (2.35, 0.61) as its only strictly positive equilibrium point.
5.1 Impact of delay
In the initial segment of this discussion, we designate the diffusion parameters as d1 = 0.1 and d2 = 0.1. In the absence of delay, for the temporal case, we establish the positivity of α0 and β0, and similarly, for the spatiotemporal case, αn and βn are both positive. This ensures that the Routh–Hurwitz conditions, referenced in Theorems 4 and 11, are met, implying local asymptotic stability of the Equation 1 around the internal equilibrium point. Additionally, based on Theorems 6 and 11, we identify the stability interval for our model as [0, 3.217]. By deliberately selecting various delay values inside and outside this stability interval, we conduct an analysis of the model's behavior to validate our theoretical findings. This examination allows us to draw conclusions regarding the model's stability under different delay scenarios. We provide a bifurcation diagram to visually illustrate the stability and instability intervals, along with the nature of bifurcations. Moving to Figures 1, 2, these diagrams depict the phase portrait and the trajectories of solutions over time for both species. Initiated from the point (1, 1), it's evident that the solutions converge toward the internal equilibrium point. This observation highlights the convergence behavior of the model's solutions under these specified conditions.
Phase portrait when τ = 0.
Time series of prey and predators when τ = 0.
The Figures 3, 4 represent the respectively the prey and the predator solution over time and space in domain [0, 100] × [0, π].
Spatiotemporal trajectory of prey when τ = 0.
Spatiotemporal trajectory of predators when τ = 0.
For τ = 2.3 ∈ [0, τ*], the system maintains its stability, which translates in Figures 5, 6.
Phase portrait when τ = 2.3.
Time series of prey and predators when τ = 2.3.
Similarly, for the spatiotemporal case, the prey and the predator populations converge to the equilibrium point, as shown in Figures 7, 8.
Spatiotemporal trajectory of prey when τ = 2.3.
Spatiotemporal trajectory of predators when τ = 2.3.
For τ = 3.5, that is to say outside the stability interval associated with the model, we notice that the Equation 1 experiences a Hopf bifurcation, it loses its stability with the appearance of periodic solutions (see Figures 9, 10).
Phase portrait of Equation 2 when τ = 3.5.
Time series of prey and predators when τ = 3.5.
Of the same for the spatiotemporal solution, which oscillates around the equilibrium point without converging. which is clear in the Figures 11, 12.
Spatiotemporal trajectory of prey when τ = 3.5.
Spatiotemporal trajectory of predators when τ = 3.5.
In concluding this section, we present Figures 13, 14, showcasing the bifurcation diagrams associated with the delay for the prey and predators, respectively. These diagrams are constructed based on the minimum and maximum amplitudes of N (prey) and P (predators). The color differentiation within the diagrams signifies the nature of stability within the model: areas marked in blue delineate the stability interval, indicating regions where the model remains stable. A red star located at τ* = 3.217 designates the bifurcation point, signifying a critical value where a qualitative change in the system's behavior occurs due to variations in delay. The green curve represents the interval of instability, indicating regions where the model exhibits instability. These bifurcation diagrams serve as visual representations, effectively capturing the complex dynamics associated with changes in delay. They provide an intuitive understanding of how the stability properties evolve concerning variations in delay, offering crucial insights into the system's behavior and transitions between stability and instability.
Bifurcation diagram of prey population.
Bifurcation diagram of predator population.
5.2 Impact of diffusion coefficient
In the subsequent segment, our focus shifts toward understanding how alterations in the diffusion coefficient d2 influence the dynamics of these populations. Employing identical parameters listed in the table and initiating the system from the point (2.358+0.001 × sin(πx), 0.609+0.001 × sin(πx)), we set τ = 0 and d2 = 0.1. This deliberate selection adheres to verifying the Routh–Hurwitz condition, crucial for assessing the system's stability. Figures 15, 16 visually represent the outcomes of this analysis, demonstrating the model's stability under these specific conditions. These figures, presenting the phase portraits and solution trajectories over time for both prey and predators, reveal the system's behavior when subjected to variations in the diffusion coefficient d2. The stability observed in these diagrams indicates the model's robustness and predictable dynamics under these prescribed parameters, allowing us to draw conclusions about the impact of d2 on the overall system behavior.
Spatiotemporal trajectory of prey when τ = 0 and d2 = 0.1.
Spatiotemporal trajectory of predators when τ = 0 and d2 = 0.1.
When altering the diffusion coefficient to d2 = 0.01, while maintaining consistency in the values of other parameters, the positive equilibrium E* remains stable for the ordinary differential equation (ODE) system. This means that the system exhibits stability and remains in an equilibrium state under these adjusted conditions. However, as per Theorem 11, the hypothesis H2 of the theorem is validated, indicating a Turing instability in the spatiotemporal domain. This instability arises because the system attains a stable, nonconstant steady-state solution, causing the previously stable positive equilibrium E* to become unstable. This transformation in stability properties is depicted in Figures 17, 18. The consequence of this instability is visually evident in the diagrams, where uneven spatial dispersion of populations is observed. This uneven distribution implies that the populations no longer maintain a homogeneous spread across space; instead, localized patterns emerge, indicative of spatial segregation or clustering within the ecosystem. This phenomenon showcases the intricate relationship between diffusion coefficients, stability properties, and spatial distribution, highlighting the system's propensity for spatially varied population distributions when subjected to specific parameter alterations.
Spatiotemporal trajectory of prey when τ = 0 and d2 = 0.01.
Spatiotemporal trajectory of predators when τ = 0 and d2 = 0.01.
When τ = 3.22 and d2 = 0.001, the system exhibits instability, fostering the existence of periodic inhomogeneous solutions. Under these conditions, a remarkable phenomenon emerges: the prey and predators coexist through spatially inhomogeneous oscillations. Notably, their densities manifest in opposing spatial distributions, a characteristic vividly visible in Figures 19, 20. These figures distinctly illustrate the spatial patterns where the densities of prey and predators display contrasting distributions within the ecosystem, signifying a spatial separation or differentiation that sustains their coexistence.
Spatiotemporal trajectory of prey when τ = 3.22 and d2 = 0.001.
Spatiotemporal trajectory of predators when τ = 3.22 and d2 = 0.001.
6 Discussion
The analysis of the delayed spatiotemporal predator-prey system provides valuable insights into the complex dynamics governing the interactions between predators and prey, incorporating various biological phenomena such as the Allee effect, fear effects on prey, cooperative hunting, the impact of toxicity, and differential fishing pressures on both species. Through this study, we have gained a better understanding of how these factors influence population dynamics by examining the stability and bifurcations of the system within both temporal and spatiotemporal frameworks.
In the temporal case, where diffusion coefficients are set to d1 = d2 = 0.1, we show that the positivity of the parameters α0 and β0 ensures the local asymptotic stability of the internal equilibrium point, as demonstrated by the model without delay. By analyzing the eigenvalue structure and verifying the Routh–Hurwitz conditions, we confirm that the equilibrium is stable for small perturbations, provided the delay remains within a specific stability interval [0, 3.217], as indicated by Theorems 4 and 11. The bifurcation diagrams (Figures 13, 14) offer a visual representation of the system's dynamics, highlighting the critical bifurcation point at τ* = 3.217, beyond which a Hopf bifurcation leads to periodic oscillations in the populations.
In the spatiotemporal case, we observe that diffusion plays a significant role in the spatial distribution of predator and prey populations. When the system operates within the stability interval, the populations converge to the equilibrium point. However, as the delay exceeds this interval, the system undergoes periodic oscillations, and spatially inhomogeneous patterns emerge. These oscillations reflect the system's inability to maintain equilibrium, resulting in non-converging, oscillatory population densities. The behavior of the system also depends on the diffusion coefficient d2, as variations in this parameter influence the stability of the system and the emergence of spatial patterns. When d2 is reduced to 0.01, the positive equilibrium remains stable within the ODE system, but in the spatiotemporal case, this reduction induces a Turing instability, leading to the formation of spatial patterns. This instability arises from the heterogeneity induced by diffusion, even though the temporal dynamics remain stable, illustrating how diffusion can destabilize spatial equilibrium.
Reducing d2 to 0.001 and increasing the delay to τ = 3.22 exacerbates the instability of the system, supporting spatially inhomogeneous periodic solutions. These oscillations lead to spatial segregation between predator and prey, with high prey densities corresponding to low predator densities, and vice versa. This spatial separation is characteristic of ecological systems where environmental factors, habitat fragmentation, or localized resource availability cause spatial differentiation between prey and predator populations.
These results highlight the resilience of the predator-prey system in maintaining coexistence, even under complex conditions of spatial segregation. It demonstrates the model's ability to predict complex ecological behaviors, such as patchiness and species clustering in shared environments, where diffusion-induced instability and temporal delays govern population dynamics. Our findings align with recent studies on diffusive predator-prey systems incorporating time delays and spatial heterogeneity [4, 34, 37], yet they also introduce new perspectives on the role of fishing pressures and toxicity. For example, in Yang et al. [4], nonlocal competition was shown to induce spatially inhomogeneous bifurcating periodic solutions. Similarly, our model reveals that reducing the diffusion coefficient d2 triggers Turing instabilities, leading to spatial pattern formation, suggesting that diffusion-driven instabilities can arise in different ecological scenarios, whether due to competition mechanisms or differential movement rates.
In Song et al. [34], the impact of time delays in a system with a generalist predator was studied, showing that delays induce oscillatory behavior. Our results corroborate this, particularly with the Hopf bifurcation observed at τ = 3.217, marking the transition from stability to periodic oscillations. However, our study extends this analysis by incorporating the combined effects of spatial diffusion and delay, demonstrating that these factors together drive more complex spatial dynamics, such as predator-prey segregation. Moreover, in Yang et al. [37], the authors examined how habitat complexity influences predator-prey interactions. While their model highlights the role of environmental structure, our approach distinguishes itself by explicitly considering how spatial heterogeneity, coupled with diffusion and delay, can generate inhomogeneous patterns without assuming pre-existing habitat constraints. This distinction is crucial for understanding how spatial structures emerge from intrinsic system dynamics, rather than being shaped by external environmental factors. This perspective is particularly relevant for understanding how real-world ecosystems, where predator and prey distributions are influenced by both internal dynamics and external pressures, function and evolve.
The results obtained from this delayed spatiotemporal predator-prey model emphasize the importance of spatial diffusion and time delays in generating complex predator-prey dynamics. These factors play a key role in determining population persistence, stability, and spatial organization, providing a more comprehensive framework for analyzing predator-prey interactions in both managed and natural ecosystems. By considering the interplay between biological and environmental factors, including toxicity and fishing pressures, we gain a deeper understanding of the ecological behaviors that emerge from complex predator-prey systems.
7 Conclusion
This study provides a detailed analysis of a delayed spatiotemporal predator-prey system, integrating various biological and environmental factors, such as the Allee effect, fear effects on prey, cooperative hunting, toxicity, and fishing pressures. Our results reveal that time delays significantly impact the system's stability, leading to periodic oscillations when delays cross a critical threshold. This finding aligns with previous studies, such as those Chakrabortyet al. [45], which observed similar oscillatory behavior in predator-prey models with delays. However, our model extends these findings by incorporating spatial diffusion, which, when altered, induces Turing instabilities and the formation of spatial patterns. While Shukla et al. [20] demonstrated how diffusion can create spatially inhomogeneous patterns in algal blooms, our work applies this concept to predator-prey interactions, showing how diffusion, combined with delays, leads to spatial segregation between predator and prey populations. Furthermore, our study introduces the impact of external pressures such as fishing and toxicity, which are absent in previous models like those of Maity et al. [19], who focused on eco-epidemic dynamics. By incorporating these factors, we provide a more comprehensive understanding of how human interventions destabilize predator-prey systems, a nuance not captured in earlier works. Lastly, our study confirms the spatial patchiness observed in ecological systems, a feature explored in Maity et al. [19], but extends it by showing that such patterns can emerge purely from the internal dynamics of the system, without relying on pre-existing habitat structures, unlike studies that emphasized habitat complexity. Overall, this work builds upon and extends previous research by integrating both internal and external factors, offering a more holistic understanding of predator-prey dynamics in complex ecosystems. In future research, this model could be extended by incorporating stochastic environmental effects and age-structured populations, which would allow for a more realistic representation of ecological uncertainty. Additionally, validating the model using empirical data from real ecosystems would enhance its practical applicability and support the development of more effective management strategies for predator-prey systems under human pressure.
Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.
Author contributions
MH: Writing – original draft. NB: Writing – original draft. YE: Supervision, Writing – original draft. NA: Supervision, Writing – original draft.
Funding
The author(s) declare that no financial support was received for the research and/or publication of this article.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Gen AI was used in the creation of this manuscript.
Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
ReferencesSongYShiQ. Stability and bifurcation analysis in a diffusive predator-prey model with delay and spatial average. Math Methods Appl Sci. (2022) 46:5561–84. 10.1002/mma.8853XiangCHuangJWangH. Bifurcations in Holling-Tanner model with generalist predator and prey refuge. J Differ Equ. (2023) 343:495–529. 10.1016/j.jde.2022.10.018VatsalaASPageniG. Series solution method for solving sequential Caputo fractional differential equations. AppliedMath. (2023) 3:730–40. 10.3390/appliedmath3040039YangRNieCJinD. Spatiotemporal dynamics induced by nonlocal competition in a diffusive predator-prey system with habitat complexity. Nonlinear Dyn. (2022) 110:879–900. 10.1007/s11071-022-07625-xYangRWangFJinD. Spatially inhomogeneous bifurcating periodic solutions induced by nonlocal competition in a predator-prey system with additional food. Math Methods Appl Sci. (2022) 45:9967–78. 10.1002/mma.8349HafdaneMBoutayebHAgmourIEl FoutayeniYAchtaichN. The role of cleaner fish in a predator-prey model: dynamics and optimal harvesting. Commun Math Biol Neurosci. (2023) 2023:92.BhattacharyyaJPalS. Stage-structured cannibalism with delay in maturation and harvesting of an adult predator. J Biol Phys. (2013) 39:37–65. 10.1007/s10867-012-9284-623860833BhattacharyyaJPalS. Hysteresis in coral reefs under macroalgal toxicity and overfishing. J Biol Phys. (2015) 41:151–72. 10.1007/s10867-014-9371-y25708511CresswellW. Predation in bird populations. J Ornithol. (2011) 152:251–63. 10.1007/s10336-010-0638-1XieB. Impact of the fear and Allee effect on a Holling type II prey-predator model. Adv Differ Equ. (2021) 2021:1–15. 10.1186/s13662-021-03592-6HaleJK. Theory of Functional Differential Equations. New York, NY: Springer (1977). 10.1007/978-1-4612-9892-2HafdaneMColleraJAAgmourIEl FoutayeniY. Hopf bifurcation for delayed prey-predator system with Allee effect. Commun Math Biol Neurosci. (2023) 2023:36. Available online at: https://scik.org/index.php/cmbn/article/view/7921BabaNIdmbarekABendahouFE. Navigating the Allee effect: unraveling the influence on marine ecosystems. J Coast Conserv. (2023) 27:65. 10.1007/s11852-023-00989-1BabaNAgmourIEl FoutayeniYAchtaichN. Toxicity impacts on bioeconomic models of phytoplankton and zooplankton interactions. Ocean Dyn. (2023). 74:53–66. 10.1007/s10236-023-01588-2ChakrabortySTiwariPKSasmalSKBiswasSBhattacharyaSChattopadhyayJ. Interactive effects of prey refuge and additional food for predator in a diffusive predator-prey system. Appl Math Model. (2017) 47:128–40. 10.1016/j.apm.2017.03.028AlleeWC. Animal aggregations. In: A Study in General Sociology. Chicago, IL: University of Chicago Press (1931). 10.5962/bhl.title.731333311142LiuJ. Qualitative analysis of a diffusive predator-prey model with Allee and fear effects. Int J Biomath. (2021) 14:2150037. 10.1142/S1793524521500376LaiLZhuZChenF. Stability and bifurcation in a predator-prey model with the additive Allee effect and the fear effect. Mathematics. (2020) 8:1280. 10.3390/math8081280MaitySSTiwariPKShuaiZPalS. Role of space in an eco-epidemic predator-prey system with the effect of fear and selective predation. J Biol Syst. (2023) 31:883–920. 10.1142/S0218339023500316ShuklaJBMisraAKSundarS. Effect of cross-diffusion on the patterns of algal bloom in a lake: a nonlinear analysis. Nonlinear Stud. (2014) 21.ChakrabortyKDasSKarTK. Optimal control of effort of a stage-structured prey-predator fishery model with harvesting. Nonlinear Anal Real World Appl. (2011) 12:3452–67. 10.1016/j.nonrwa.2011.06.007AngTKSafuanHM. Dynamical behaviors and optimal harvesting of an intraguild prey-predator fishery model with Michaelis-Menten type predator harvesting. Biosystems. (2021) 202:104357. 10.1016/j.biosystems.2021.10435733486091HuPCaoHJ. Stability and bifurcation analysis in a predator-prey system with Michaelis-Menten type predator harvesting. Nonlinear Anal Real World Appl. (2017) 33:58–82. 10.1016/j.nonrwa.2016.05.010GuoHJChenLS. The effects of impulsive harvest on a predator-prey system with distributed time delay. Commun Nonlinear Sci Numer Simul. (2009) 14:2301–9. 10.1016/j.cnsns.2008.05.010DouJLiS. Optimal impulsive harvesting policies for single-species populations. Appl Math Comput. (2017) 292:145–55. 10.1016/j.amc.2016.07.027LiWJiJHuangL. Global dynamic behavior of a predator-prey model under ratio-dependent state impulsive control. Appl Math Model. (2020) 77:1842–59. 10.1016/j.apm.2019.09.033TangSYPangWHChekeRAWuJH. Global dynamics of a state-dependent feedback control system. Adv Diff Equ. (2015) 2015:322. 10.1186/s13662-015-0661-xZhangTMaWMengXZhangT. Periodic solution of a prey-predator model with nonlinear state feedback control. Appl Math Comput. (2015) 266:95–107. 10.1016/j.amc.2015.05.016HuangMZLiuSZSongXYChenLS. Dynamics of unilateral and bilateral control systems with state feedback for renewable resource management. Complexity. (2020) 2020:9453941. 10.1155/2020/9453941HockingLM. Optimal Control: An Introduction to the Theory with Applications. Oxford: Oxford University Press (1991). 10.1093/oso/9780198596752.001.0001XiaoZLiZ. Stability analysis of a mutual interference predator-prey model with the fear effect. J Appl Sci Eng. (2019) 22:205–11.SasmalSKTakeuchiY. Dynamics of a predator-prey system with fear and group defense. J Math Anal Appl. (2020) 481:123471. 10.1016/j.jmaa.2019.123471BendahouFEBabaNEl FoutayeniYAchtaichN. Impact of pollution on sardine, sardinella, and mackerel fishery: a bioeconomic approach. Commun Math Biol Neurosci. (2023) 2023:55.SongYPengYZhangT. The spatially inhomogeneous Hopf bifurcation induced by memory delay in a memory-based diffusion system. J Differ Equ. (2021) 300:597–624. 10.1016/j.jde.2021.08.01040439727AkhmetMUBekliogluMErgencTTkachenkoVI. An impulsive ratio-dependent predator-prey system with diffusion. Nonlinear Anal Real World Appl. (2006) 7:1255–67. 10.1016/j.nonrwa.2005.11.00729962293LiuYWeiJ. Double Hopf bifurcation of a diffusive predator-prey system with strong Allee effect and two delays. Nonlinear Anal Model Control. (2021) 26:72–92. 10.15388/namc.2021.26.20561YangRJinDWangW. A diffusive predator-prey model with generalist predator and time delay. AIMS Math. (2022) 7:4574–91. 10.3934/math.2022255AlikakosND. An application of the invariance principle to reaction-diffusion equations. J Differ Equ. (1979) 33:201–25. 10.1016/0022-0396(79)90088-3NiWMTangM. Turing patterns in the Lengyel-Epstein system for the CIMA reaction. Trans Am Math Soc. (2005) 357:3953–69. 10.1090/S0002-9947-05-04010-930656504BanerjeeMTakeuchiY. Maturation delay for the predators can enhance stable coexistence for a class of prey-predator models. J Theor Biol. (2017) 412:154–71. 10.1016/j.jtbi.2016.10.01627825813MengXJiaoJChenL. The dynamics of an age-structured predator-prey model with disturbing pulse and time delays. Nonlinear Anal Real World Appl. (2008) 9:547–61. 10.1016/j.nonrwa.2006.12.001HafdaneMAgmourIEl FoutayeniY. Study of Hopf bifurcation of delayed tritrophic system: dinoflagellates, mussels, and crabs. Math Model Comput. (2023) 10:66–79. 10.23939/mmc2023.01.066JoyRAHawladerMSRahmanMSHossainMRShamimSIMahmudH. Improving quality, productivity, and cost aspects of a sewing line of apparel industry using TQM approach. Math Probl Eng. (2024) 2024:6697213. 10.1155/2024/6697213La TorreDKunzeHRuiz-GalanMMalikTMarsiglioS. Optimal control: theory and application to science, engineering, and social sciences. Abstr Appl Anal. (2015) 2015:890527. 10.1155/2015/890527ChakrabortySTiwariPKMisraAKChattopadhyayJ. Spatial dynamics of a nutrient–phytoplankton system with toxic effect on phytoplankton. Math Biosci. (2015) 264:94–100. 10.1016/j.mbs.2015.03.01025843351StephensPASutherlandWJ. Consequences of the Allee effect for behavior, ecology, and conservation. Trends Ecol Evol. (1999) 14:401–5. 10.1016/S0169-5347(99)01684-510481204BerecLAnguloECourchampF. Multiple Allee effects and population management. Trends Ecol Evol. (2007) 22:185–91. 10.1016/j.tree.2006.12.00217175060CantrellRSCosnerCHutsonV. Permanence in ecological systems with spatial heterogeneity. Proc R Soc Edinb A Math. (1993) 123:533–59. 10.1017/S0308210500025877ZanetteLYWhiteAFAllenMCClinchyM. Perceived predation risk reduces the number of offspring songbirds produce per year. Science. (2011) 334:1398–401. 10.1126/science.121090822158817CourchampFClutton-BrockTGrenfellB. Inverse density dependence and the Allee effect. Trends Ecol Evol. (1999) 14:405–10. 10.1016/S0169-5347(99)01683-310481205