Chaos detection in predator-prey dynamics with delayed interactions and Ivlev-type functional response

1.   Introduction

Investigating the dynamic interaction and interplay between species is essential in mathematical ecology ^[1,2]. Modeling such systems and analyzing their dynamical behavior may give a prediction on the evolution of populations. Particularly, two-species predator-prey models have led to enthusiasm among many scholars ^[3,4,5]. The Gause type two-species predator-prey model is given by

where $x(t)$ and $y(t)$ denote the density of the prey and predator species at time $t$, respectively. Parameters $r, K, s,$ and $Y$ are positive constants, which denote the intrinsic reproduction rate of the prey, the carrying capacity for the prey species, the death rate of the predator species, and the growth yield constant for the conversion of prey to predator density, respectively ^[6,7].

The initial values are $x(0)\geq 0, y(0)\geq 0$ due to their biological meanings. In the absence of $y(t)$ in the model (1.1), the prey increases according to the logistic growth $\dot{x}(t) = rx(t)(1-\frac{x(t)}{K})$. The predator $y(t)$ declined exponentially as $\dot{y}(t) = -sy(t)$ and will eventually die in the long-term if the model lacks prey $x(t)$.

The function $f(x)$ represents the prey-dependent functional response, which is the Ivlev-type functional response ^[8,9], and takes the form

where $\alpha > 0$ is the consumption rate and $\beta > 0$ is the physiological rate at which saturation is achieved. It is a monotone increasing function that saturates, that is, it has a finite positive limit $\alpha$ as $x$ approaches infinity ^[1] (see Figure 1).

This type functional response is both monotonically increasing and uniformly bounded, which was classified to the Holling-Ⅱ functional response by Garvie ^[10,11]. Biologically, it was first proposed to describe the increase of the fish. Hence, our results are useful in designing fishing policies for the fishery industry. Other forms of functional responses can be seen in ^[12,13].

Considering the Ivlev-type trophic response (1.2) in system (1.1), it changes to

The dynamical behaviors of the model (1.3) have been investigated extensively. It has the only coexistence equilibrium $(x^{\ast}, y^{\ast})$ if

where

After applying rescalings in $x, y$, and $t$, it can be assumed that $K = \alpha = Y = 1$. Under the assumption (1.4), if

then there exists a unique stable limit cycle. Otherwise, it has no limit cycles. If (1.4) fails, then system (1.3) has no existence equilibrium. Therefore, no limit cycles of system (1.3) exist.

In 1925, to investigate fish population under harvesting, the predator-prey model with delay was proposed by Volterra. It is described by an integro-differential equation as

The above delayed equation is called an integro-differential equation, and such delays are called distributed delays. We can use the linear chain trick to convert systems into systems with discrete delay. Since then, delayed differential equations (DDEs) have been extensively used to model population dynamics ^[14], neural network ^[15,16], engineering, the life sciences, etc., including predator-prey interactions.

By ^[17], assume the growth rate of the predator species $y(t)$ is proportional to the number of individuals in the population $t-\tau$ time units in the past that manage to survive until time $t$. In order to obtain an expression that describes how many predator individuals alive at time $t-\tau$ are still alive at time $t$, where $\tau$ is the delay due to the gestation of the $y(t)$ ^[18], we need to solve the following first-order ordinary differential equation for $y(t)$,

It implies that

hence

where the factor $e^{-s\tau}$ denotes the survival rate of the predator $y(t)$ which was born at time $t-\tau$ and still remains alive at the time $t$. When the time delay $\tau = 0$, the right side of (1.5) reduces to its prototype $y(t)$. The main difference between $y(t-\tau)$ and $y(t-\tau)e^{-s\tau}$ is shown in Figure 2.

When the time delay $\tau$ is small, they are close ^[19]. However, as the delay increases, the expression $y(t-\tau)e^{-s\tau}$ could describe practical problems better than the expression $y(t-\tau)$. Although the predator-prey interaction with $y(t-\tau)$ has been extensively investigated since its proposal, such systems with the delay term $y(t-\tau)e^{-s\tau}$ are scarce and are not frequently reported. Therefore, compared to existing studies on the predator-prey system ^[20,21,22], this is the main contribution and the novelty of this paper in the aspect of establishing the model.

Similarly, with the well-known Wangersky-Cunningham model, we assume that the change rate of predators depends on the number of prey and of predators present at $\tau$ previous time, that is, the delay $\tau$ in the interaction term $y(t)(1-e^{-\beta x(t)})$ of the second equation. Therefore, we replace $y(t)$ with the Eq (1.5) in model (1.3), and system (1.3) is reduced to

System (1.6) has the initial data

where $(\phi(\eta), \psi(\eta)) \in C\left([-\tau, 0], \mathbb{R}_{+0}^{2}\right)$ is the Banach space of continuous functions mapping the interval $[-\tau, 0]$ into $\mathbb{R}_{+0}^{2}$, where $\mathbb{R}_{+0}^{2} = \{(x, y): x \geq 0, y \geq 0\}$ ^[23]. By the fundamental theory of DDEs, system (1.6) has a unique solution $x(t), y(t)$ satisfying initial data (1.7).

The main goal of this paper is to show how the delay $\tau$ affects the dynamics of model (1.6). This paper is organized as follows. In Section 2, we prove the positivity and boundedness for the solution of system (1.6). In Section 3, when time delay is accounted as a bifurcation parameter, the stability analysis is given for the coexistence equilibrium for model (1.6). We analytically prove that the local Hopf bifurcation critical values are neatly paired. In Section 4, numerical explorations using the numerical continuation software XPPAUT and DDE-Biftool are carried out in order to substantiate the obtained theoretical results. Simulations indicated that as the delay increases, the positive equilibrium loses its stability and bifurcates a family of orbitally asymptotically stable periodic solutions. The coexistence equilibrium undergoes stability switches. For large enough delay, the predator will die out. Before the extinction of the predator, rich dynamics such as Hopf bifurcation, period doubling bifurcation and strange attractor have been demonstrated when time delay is accounted as bifurcation parameter, and the abundance of steady-state chaotic solutions appears via a cascade of period-doubling bifurcations is also detected. The coexistence equilibrium undergos transcritical bifurcation at the die out critical value. We summarize our conclusions in Section 5, especially on the impact of time delay from the biological aspect.

2.   Positivity and boundedness

For the system (1.6) with positive initial data (1.7), we first prove the following two theorems concerning the positivity and boundedness of the solution ^[24].

Theorem 2.1. Solutions of system (1.6) with positive initial data (1.7) remain positive for $t > 0$.

Proof. Assume $(x(t), y(t))$ is a solution of system (1.6) satisfying initial data (1.7). By the Theorem 2.1 of ^[25], we solve the following ordinary differential equation (ODE):

Integrating between the limit from $0$ to $t$, the solution is

Obviously, the exponential function is always positive, regardless of the integrand. It implies that $x(t)$ is positive for $t > 0$ and $\phi(0) > 0$.

Next, we show that $y(t)$ is positive on $t\in[0, +\infty)$. Based on the theory of Hale ^[26], it is obvious that $y(t)$ is well-defined on $[-\tau, +\infty)$ and

Since $\varphi(0) > 0$ and initial data (1.7), we have $y(t) > 0$ when $t\in[0, \tau]$, therefore $y(t) > 0$ for $t\in[0, +\infty]$. Positivity implies that the cone of the solutions is invariant in the system. □

Theorem 2.2. Solutions of system (1.6) with positive initial data (1.7) are uniformly ultimately bounded.

Proof. Define the following function:

The derivative of $\omega(t)$ with respect to $t$ is

Substituting $\dot{x}(t)$ and $\dot{y}(t+\tau)$ into the above expression, we obtain

where $x_{0}$ is the upper bound of $x(t)$. By the Lemma 2.1 of ^[27,28], we obtain

when $t$ is sufficiently big. It implies that $x(t)$ and $y(t)$ are ultimately bounded. The boundedness of the model ensures that there is a restriction on the growth of the species due to limited resources in nature. This completes the proof. □

3.   Local stability and bifurcation analysis

To begin, we consider the possible equilibria of model (1.6).

Proposition 3.1. (i) System (1.6) has two distinct equilibria, the trivial equilibrium $E_{0} = (0, 0)$ and the semi-trivial equilibrium $\bar{E} = (K, 0)$.

(ii) If

holds, the $\tau_{c}$ is the critical value, then system (1.6) has a coexistence equilibrium $E^{\ast} = (x^{\ast}, y^{\ast})$, where

The existence of $E^{\ast}$ ensures that $1-e^{-\beta K} > 0$. Note that the equilibrium value depends on $\tau$: $x^{\ast}$ is an increasing function with respect to the delay $\tau$, while $y^{\ast}$ is a decreasing function when $x^{\ast} > \frac{K}{2}$, that is: The larger the delay, the higher the number of the prey population, and the lower the number of predators at the equilibrium.

3.1.   Linear stability of $E_{0}$ and $\bar{E}$

The characteristic equation corresponding to $E_{0} = (0, 0)$ is

whose roots are obtained as $\lambda_{1} = r > 0$.

Lemma 3.1. For all $\tau\geq0$, the trivial equilibrium $E_{0}$ is always unstable.

The characteristic equation with respect to $\bar{E} = (K, 0)$ is

implying that $\lambda_{1} = -r < 0$ and

Let

Therefore,

for any $\tau\geq0$. Thus, if $\tau < \tau_{c}$, $f(0) < 0$, and $f(\lambda) = 0$ has at least one positive root. Thus, when $\tau < \tau_{c}$, the semi-trivial equilibrium is unstable.

Lemma 3.2. When $\tau < \tau_{c}$, $\bar{E}$ is unstable.

3.2.   Linear stability and Hopf bifurcation of $E^{\ast}$

In this part, assume that (3.1) is satisfied, then the interior (coexistence) equilibrium $E^{\ast}$ exists. The linearized system of (1.6) about $E^{\ast}$ is

where $X(t) = (x(t), y(t))^{T}$, $X(t-\tau) = (x(t-\tau), y(t-\tau))^{T}$, $A_{0} = \left[\begin{array}{ccc} r(1-\frac{2x^{\ast}}{K})+\alpha\beta y^{\ast}e^{-\beta x^{\ast}} & -\alpha(1-e^{-\beta x\ast}) \\ 0 & -s \end{array}\right]$, $A_{1} = \left[\begin{array}{ccc} 0 & 0 \\ \alpha\beta Yy^{\ast}e^{-s\tau-\beta x^{\ast}} & \alpha Ye^{-s\tau}(1-e^{-\beta x^{\ast}}) \end{array}\right]$.

The linearization system (3.2) around $E^{\ast}$ has the following characteristic equation:

that is,

where

Notice that if $\tau = 0$, Eq (3.3) reduces to the second-order polynomial equation

and it follows that all eigenvalues of Eq (3.5) have negative real parts if, and only if,

The transcendental Eq (3.3) has infinitely many roots. Note that polynomials $p_{i}(\tau)(i = 1, 2, 3, 4)$ are dependent on $\tau$. The transcendental equation associated with (3.2) at $E^{\ast}$ is

where

For the characteristic equation, before applying the criterion due to Beretta and Kuang ^[29] to evaluate the existence of a purely imaginary root, we first verify the following properties for $\tau\in[0, \tau_{c})$, where $\tau_{c}$ is the maximum value when $E^{\ast}$ exists.

$({\rm{i}})\; P(0, \tau)+Q(0, \tau)\neq0;$

$({\rm{ii}})\; P({\rm{i}}\omega, \tau)+Q({\rm{i}}\omega, \tau)\neq0, \forall\; \omega\in\mathbb{R};$

$({\rm{iii}})\; \lim\limits_{|\lambda|\rightarrow \infty} {\rm{sup}}\big\{|\frac{Q(\lambda, \tau)}{P(\lambda, \tau)}|; {\rm{Re}} \lambda\geq0\big\} < 1;$

$({\rm{iv}})\; F(\omega, \tau): = |P({\rm{i}}\omega, \tau)|^{2}-|Q({\rm{i}}\omega, \tau)|^{2}$ has a finite number of zeros;

$({\rm{v}})$ Each positive root $\omega(\tau)$ of $F(\omega, \tau) = 0$ is continuous and differentiable in $\tau$ whenever it exists ^[18].

Obviously,

(ⅰ) is satisfied. Assumption ensures that $\lambda = 0$ is not the root of Eq (3.7).

Assume that $p_{2}(\tau)+e^{-s\tau}p_{4}(\tau)\neq\omega^{2}, \; p_{1}(\tau)+e^{-s\tau}p_{3}(\tau)\neq0$, then

It follows from (3.7) that

hence (ⅲ) follows.

For the function $F$ defined in (ⅳ), which follows

and

such that

where

$F(\omega, \tau)$ has at most four roots, therefore assumption (iv) is satisfied. Assume that $(\omega_{0}, \tau_{0})$ is a point in its domain such that $F(\omega_{0}, \tau_{0}) = 0$. In a certain neighborhood of $(\omega_{0}, \tau_{0})$, the partial derivatives $F_{\omega}$ and $F_{\tau}$ exist and are continuous, and $F_{\omega}(\omega_{0}, \tau_{0})\neq0$. Assumption (ⅴ) is satisfied by the implicit function theorem. Assumption (ⅳ) guarantees that Eq (3.7) has at most a finite number of purely imaginary roots ^[18], i.e., the roots cross the imaginary axis a finite number of times as $\tau$ varies.

Next, we assume that $\lambda = {\rm{i}}\omega(\omega > 0, {\rm{i}} = \sqrt{-1})$ is the pure imaginary root of expression (3.7), then $\lambda = {\rm{i}}\omega$ satisfies

i.e., because $|e^{- {\rm{i}}\omega\tau}| = 1$, $\omega(\tau)$ is the positive zero root of

Define the set

Therefore,

has positive root $\omega = \omega(\tau)$ if $\tau \in I$. Otherwise, $F(\omega, \tau) = 0$ does't have a positive zero point. Furthermore, we obtain

In addition, define the function $\theta(\tau)\in[0, 2\pi]$, which satisfied (3.9) for $\tau \in I$, i.e.,

The $\omega(\tau)\tau$ in (3.9) and $\theta(\tau)$ in (3.10) have the following relationship:

Introduce map $\tau_{n}:I\rightarrow \mathbb{R}_{0}^{+}$:

where $\omega(\tau)$ is the positive root of (3.8). From $I$ to $\mathbb{R}$, define the continuous and differential function $S_{n}(\tau)$ as

Let $\lambda(\tau)$ be the eigenvalues satisfied by $\lambda(\tau^{\ast}) = {\rm{i}}\omega(\tau^{\ast})$, and the transversality condition is obtained as

Theorem 3.1. (i) For model (1.6), if either the set $I$ is empty or the function $S_{n}(\tau)$ has no positive zero in $I$, for $0 < \tau < \tau_{c}$, the positive equilibrium $E^{\ast}$ is asymptotically stable.

(ii) If Eq (3.11) has positive roots in $I$ denoted by $\{\tau_{1}, \tau_{2}, \cdots, \tau_{m}\}$ with $\tau_{j} < \tau_{j+1}$ and $S'_{n}(\tau_{1}) > 0$, the positive equilibrium $E^{\ast}$ is asymptotically stable for $\tau\in[0, \tau_{1})\cup(\tau_{m}, \tau_{c})$ and unstable for $\tau\in(\tau_{1}, \tau_{m})$, with Hopf bifurcations occurring when $\tau = \tau_{j}, j = 1, 2, \cdots, m.$

4.   Numerical simulations

To illustrate the analytical local Hopf Bifurcation results, we shall present some numerical simulations in this section and will extend them further with the help of numerical bifurcation analysis. The graphs are mainly drawn using DDEBifTool ^[30,31]. Similar dynamics have been numerically detected in the discrete delay system with Holling type Ⅱ and Beddington-DeAngelis trophic response ^[19].

4.1.   Hopf bifurcation

Hereafter, parameters are fixed at the following values,

According to the biological meaning of the parameters, because $r$ is relatively large, it indicates that the prey has a high breeding rate. The $s$ is relatively small which indicates that the predator has a low death rate. In addition, $r$ is smaller than $\alpha$ to some extent, and it indicates that the changes in the number of prey are less influenced by their own birth rate and more influenced by their ability to evade natural enemies to prevent predation. $Y\alpha$ is larger than $s$, and it indicates that the changes of the number of predators are more influenced by their ability to prey on food.

Under (4.1), we consider

The $\phi(0)$ is initial prey population, and it represents the number of prey at the start of the model (4.2). This value can affect the food supply available to predators. If the initial number of prey is low, predators may face a food shortage. In contrary, if the initial number of prey is high, predators may have an abundant food supply. Furthermore, the $\varphi(0)$ is initial predator population, and it represents the number of predators at the start. This value can influence the initial pressure that predators exert on prey. If the initial number of predators is low, the pressure on prey may be relatively small; conversely, if the initial number of predators is high, the pressure on prey may be greater. The $\phi(0)$ and $\varphi(0)$ are chosen as 0.1 here, and they are at the median level relatively. The growth and mortality parameters represent the biological characteristics of predators and prey, such as growth rates and mortality rates. These parameters can influence the population dynamics of predators and prey, thereby affecting the stability of the ecosystem.

System (4.2) has three equilibria: $E_{0} = (0, 0)$, $\bar{E} = (1, 0)$, and the coexistence equilibrium $E^{\ast} = (x^{\ast}, y^{\ast})$ exists when $\tau < \tau_{c} = 250$.

Consider Figures 3 and 4: The blue solid line (the red dotted line) represents stable equilibrium (unstable equilibrium) and the filled green circle (open blue circles) indicates stable periodic orbit (unstable periodic orbits).

Figure 3 indicates that the $E_{0} = (0, 0)$ is unstable. When $\tau\in[0, \tau_{3})$, $\bar{E} = (1, 0)$ is unstable. When $\tau > \tau_{3}$, it is locally asymptotically stable, i.e., at $\tau = \tau_{3}$ the predator goes extinct. The coexistence equilibrium $E^{\ast}$ is stable when $\tau\in[0, \tau_{1})\cup(\tau_{2}, \tau_{3})$ while it's unstable when $\tau\in(\tau_{1}, \tau_{2})\cup(\tau_{3}, \tau_{c})$. At $\tau = \tau_{3}$, the semi-trivial equilibrium $\bar{E} = (1, 0)$ and positive equilibrium $E^{\ast}$ exchange their stability, leading to the appearance of transcritical bifurcation. Note that $\tau_{1} = 4.4730$, $\tau_{2} = 69.3061$, and $\tau_{3} = 132.9233$.

We produced a corresponding bifurcation diagram as Figure 4, using $\tau$ as the primary bifurcation parameter. Figure 4 shows the first critical value $\tau_{1}$ with $\omega_{1} = 0.1215$, which the bifurcated Hopf bifurcation and the second critical value is $\tau_{2}$ with $\omega_{2} = 0.0289$, generating the supercritical Hopf bifurcation at $\tau_{1}$ and $\tau_{2}$. Biologically, above phenomenon could be interpreted as there being an interval ($\tau_{1}$, $\tau_{2}$) of survival that may exist even though the positive equilibrium is unstable.

Figures 5–8 show the trajectories and phase graph of system (4.2) with $\tau = 4.3$, $\tau = 4.5$, $\tau = 50$, and $\tau = 70.5$, respectively. Figure 5 illustrates that the coexistence equilibrium $E^{\ast} = (0.0669, 1.8724)$ is locally asymptotically with $\tau = 4.3 < \tau_{1}$. It will lose its stability and a bifurcating periodic solution occurs once $\tau = 4.5 > \tau_{1}$ as the time delay increases, as shown in Figure 6. Figure 7 indicates that a stable periodic solution occurs with $\tau = 50$. Additionally, Figure 8 shows that the coexistence equilibrium $E^{\ast}$ is locally asymptotically stable with $\tau = 70.5 > \tau_{2}$. The results are coincident with Figures 3 and 4.

The above simulations indicate that there exists a unique global Hopf bifurcation connecting $\tau_{1}$ and $\tau_{2}$. The global Hopf bifurcation is bounded, and each global Hopf branch connects a pair of Hopf bifurcation values. In the next subsection, we will detect the global Hopf bifurcation ^[32].

4.2.   Chaotic solutions

By the global Hopf bifurcation result of ^[33,34], it shows that for the following delayed Lotka-Volterra model,

After the second critical value, the local Hopf bifurcation implies a global Hopf bifurcation. However, for the delayed model (1.6), it is not valid. We will show that as follows.

Parameters are fixed at the following values,

According to the biological meaning of the parameters, because $r$ is relatively large, it indicates that the prey has a high breeding rate. The $s$ is relatively small, which indicates that the predator has a low death rate. Because $r$ is larger than $\alpha$, to some extent, and it indicates that the changes in the number of prey are more influenced by their own birth rate. In addition, $Y\alpha$ is larger than $s$, and it indicates that the changes in the number of predator are more influenced by their ability to prey on food ^[35].

Under (4.3), we consider the following model:

Consider Figures 9 and 10, where the open blue circles stand for unstable periodic orbits. Figure 9 indicates that the $E_{0} = (0, 0)$ is always unstable. The $\bar{E} = (1, 0)$ is unstable when $\tau\in[0, \tau_{3})$. When $\tau > \tau_{3}$, it is locally asymptotically stable where $\tau_{3} = 130.2649$. In $\tau\in[0, \tau_{1})\cup(\tau_{2}, \tau_{3})$, the positive equilibrium $E^{\ast}$ is stable, and unstable when $\tau\in(\tau_{1}, \tau_{2})\cup(\tau_{3}, \tau_{c})$ where $\tau_{1} = 2.6330$ and $\tau_{2} = 79.8515$. By Figure 10, the unstable periodic orbits appear between $\tau_{4}$ and $\tau_{5}$, where $\tau_{4} = 67.2228$ and $\tau_{5} = 70.4188$. Biologically, since the appearance of unstable bifurcating periodic solutions, the two species could coexist in a chaotic mode.

To see the influence of delay $\tau$ on the dynamical behaviors of the model, we detect the complex dynamical behavior when $\tau\in(\tau_{4}, \tau_{5})$ by Figures 11–14. By Figure 11, for the system (4.4), when $\tau = 65.9$, the system has a limit cycle whose period is approximately $250$. The periodic orbits are always stable until $\tau < 66.6$. When $\tau > 66.6$, stable periodic solutions undergo period-doubling bifurcation; as Figure 12 shows, when $\tau = 66.6$, the system bifurcates twice the period. When $\tau = 69.5$, the system bifurcates with a sequence of period-doubling bifurcations. When $\tau$ continues to increase from $65.9$ to $70$, by Figure 14, the system (4.4) achieves chaotic oscillation through period-doubling bifurcation with a chaotic attractor. In a stable equilibrium, the periodic oscillation by $2$, $2^{2}$, $2^{3}\cdots$ cycles eventually lead to chaos. Eventually, a cascade of period doubling bifurcations leads to chaos, which resembles the chaotic attractor of the following Mackey-Glass equation ^[36]

The existence of chaotic solutions implies that even a small environmental or parameter perturbation can disrupt the dynamics of the system ^[37,38]. In current research on the existence of chaos, only the phase diagram or the time course diagram of the system is generally provided, with few discussions on the chaotic path. In fact, the chaotic path can clearly illustrate the dynamic transition process of the system under the influence of parameters. Therefore, compared with existing results on the dynamical behaviors of predator-prey systems ^[21,39,40], the detection of chaos and the analysis of the chaotic path are the main contributions of this paper in terms of dynamical behaviors.

Remark 4.1. Since chaos is sensitive to initial conditions, the strange attractor is not obvious under the parameter values as specified in (4.3), which results in the critical values $\tau_4-\tau_6$ in Figure 10 being not accurate enough. Furthermore, the interval from $\tau_4$ to $\tau_6$ is not long enough, and we did not detect the dynamical behaviors in detail within this interval. We hypothesize that within this interval, as the delay $\tau$ increases from $70$ to $\tau_6$, the system (4.4) undergoes four cycle bifurcations, doubling the period twice, leading to a period doubling bifurcation and a limit cycle. The dynamic behavior of the system does not change substantially.

5.   Conclusions

In this paper, a delayed prey-predator model with an Ivlev-type functional response is investigated, focusing on the effect of the delay on the dynamical behaviors of the model. The supercritical Hopf bifurcation and period doubling types of bifurcations, as well as a strange attractor, can occur at the positive equilibrium when time delay is considered as a bifurcation parameter. The chaotic attractor appears, followed by a sequence of period-doubling bifurcations for small enough of the death rate of the predator species. This study delves into the intricate interplay where time delays and nonlinear responses converge, offering a deeper insight into the chaotic behaviors that may arise within these complex systems.

From a biological perspective, there are intriguing explanations. If delay is minimal, predator and prey populations stabilize. However, as delay escalates, species exhibit asymptotic, periodic, or quasi-periodic fluctuations, suggesting an oscillatory coexistence of predator and prey. As the time delay continues to increase, the system will exhibit a chaotic phenomenon of 'lose a millimeter, miss a thousand miles'. Consequently, short-term observations can be deceptive in forecasting due to the presence of bifurcation and chaos, highlighting the complexity of long-term ecological dynamics. The results could be very essential for biologists who work with delayed prey-predator systems. In conclusion, this paper makes two contributions: the introduction of an exponential delay term and the detection of chaos.

In reality, the processing time delay rarely has the same length at every instance; instead, it follows a distribution with some mean value. Our follow-up work will investigate the dynamical behaviors of the model incorporating distributed delay and compare the dynamics resulting from using various distributions, including discrete delay.

Author contributions

Qinghui Liu: The data analysis and wrote the paper; Xin Zhang: The formal analysis and the validation. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the Basic Science (Natural Science) Research Project of Jiangsu Province (No. 24KJB110010), the Special Reform and Development Project of Nanjing University of Finance and Economics in 2023 (No. XGFB3202311), Teaching Reform Project of Nanjing University of Finance and Economics in 2023 (No. JG23902).

Conflict of interest

All authors declare no conflicts of interest in this paper.