Relaxation Oscillations and Dynamical Properties in a Time Delay Slow–Fast Predator–Prey Model with a Piecewise Smooth Functional Response

Abstract

In the past few decades, the predator–prey model has played an important role in the dynamic behavior of populations. Many scholars have studied the stability of the predator–prey system. Due to the complex influence of time delay on the dynamic behavior of systems, time-delay systems have garnered wide interest. In this paper, a classical piecewise smooth slow–fast predator–prey model is considered. The dynamic properties of the system are analyzed by linearization. The existence and uniqueness of the relaxation oscillation are then proven through the geometric singular perturbation theory and entry–exit function. Finally, a stable limit cycle is obtained. A numerical simulation verifies our results for the systems and shows the effectiveness of the method in dealing with time delays.

1. Introduction

The relationship between predator and prey is an important topic both in mathematics and biology. Functional response [1,2,3] refers to a response in which the predation rate of each predator varies with prey density, namely, the predation effect of a predator on a prey. The predator–prey model with piecewise smooth response functions not only has the dynamic characteristics of the smooth system but also has many complex dynamic behaviors. Compared with the smooth system, the non-smooth system shows many bifurcation phenomena. In 2021, Li et al. [4] studied a class of predator–prey models with piecewise smooth response functions: the specific model is as follows:

$$\begin{array}{c}\frac{dx}{dt}=x(1-\frac{x}{K})-p(x,y)y,\hfill \\ \frac{dy}{dt}=\epsilon y(p(x,y)-\alpha ),\hfill \end{array}$$

(1)

where

$$\begin{array}{c}p(x,y)=\frac{1}{2+hy}\frac{x+2-\left|x-2\right|}{2}.\end{array}$$

In System (1), parameters K, $\alpha$, h are positive numbers with specific biological meanings. The parameter $\epsilon$, which represents the ratio of reproductive rate of predator to prey, can be a small parameter in some models, such as a predator lynx and its prey hare.

Ref. [4] mainly considered the case where $\epsilon$ is small, in which case in System (1) is a slow–fast system. We can study it using geometric singular perturbation theory. As for the study of geometric singular perturbation theory and exchange lemma, readers can refer to [5,6,7,8,9,10,11,12,13,14,15,16].

In recent years, the complex dynamic behaviors of fast–slow dynamical systems have been extensively investigated. In 2016, Prohens et al. [17] proved a Fenichel-like theorem, and they derived the explicit expression of the invariant slow manifold. In 2018, Wang et al. [18] considered a planar $C^{\infty }$ differential model through geometric singular perturbation theory and entry–exit function. Detailed research on the entry-e-xit functions can be found in Reference [19]. In 2020, Wang et al. [20] found relaxation oscillations in several classes of predator–prey models with time delays; Karl et al. [21] classified bifurcations of the critical set and investigated the connection between bifurcation of the critical set and general bifurcation of the attractive relaxation oscillations; Shen et al. [22] considered a predator–prey system with evolutionary effects and analyzed the fast–slow dynamic behaviors of the system by geometric singular perturbation theory; Giné et al. [23] discussed the existence of a modified Leslie–Gower model around the positive equilibrium point. In 2021, Dehingia et al. [24] considered a discrete cancer model with time delay. They studied the effect of parameter variation and time delay on dynamic behavior and proved that Hopf bifurcation occurs at positive equilibria. The direction of Hopf bifurcation and the stability of bifurcation periodic solution are obtained using normal theory and the central manifold theorem developed by Hassard et al. [25]; Li et al. [26] studied a Leslie–Gower model with piecewise smooth response functions based on geometric singular perturbation theory. In 2022, Das et al. [27] studied the dynamic behavior of cancer mathematical model with time delay and analyzed the stability of the equilibria of the system with time delay. Due to the addition of time delay, the coexistence equilibrium lost its stability and Hopf bifurcation occurred. Finally, the theoretical results were verified by numerical simulation; Khan et al. [28] investigated the predator–prey model with three species and analyzed the local stability and bifurcation behavior at the equilibrium point; Shang et al. [29] proposed a predator–prey model with a simplified Holling-IV functional response to study the bifurcation behavior of the model.

We added a time delay into System (1) to study the dynamic behavior of the predator–prey system with time delay. Lemma 1 was obtained by calculating the equilibrium of the system (7). Lemmas 2 and 3 were obtained by verifying the conditions of Hopf bifurcation in the system (7). It can be proven that Hopf bifurcation occurs at a positive equilibria. By studying the dynamic behavior of the limit system, it was found that the original system has a unique and stable limit cycle.

The rest of this paper is organized as follows. In Section 2, we give a brief overview of the entry–exit function on the plane. In Section 3, we study the equilibria, Hopf bifurcation and dynamics of limit systems about the time delay system. In Section 4, we give a theorem to explain the existence and uniqueness of the stable relaxation oscillation cycle. In Section 5, a numerical example is given to verify our theoretical results. In Section 6, the conclusions are addressed.

2. Entry–Exit Function

Consider a slow–fast vector field in the form of

$$\begin{array}{c}\frac{dx}{dt}=xf(x,y,\epsilon ),\hfill \\ \frac{dy}{dt}=\epsilon g(x,y,\epsilon ),\hfill \end{array}$$

(2)

where $(x,y)\in {\mathbb{R}}^2$ are state-space variables and the parameter $\epsilon$ represents the ratio of time scales. Let $\epsilon \to 0$. We then get the equilibria of System (2) (See Figure 1a). We define functions f and g satisfying

$$\begin{array}{c}f(0,y,0)<0,fory>0,\hfill \\ f(0,y,0)>0,fory<0,\hfill \\ g(0,y,0)<0.\hfill \end{array}$$

(3)

According to [16], we can define $p_{\epsilon }\left(y_0\right)$, which satisfies $\underset{\epsilon \to 0}{lim}{p}_{\epsilon }\left(y_0\right)=p_0\left(y_0\right)$ (See Figure 1b), where $p_0\left(y_0\right)$ is determined by

$$\begin{array}{c}{\int }_{y_0}^{p_0\left(y_0\right)}\frac{f(0,y,0)}{g(0,y,0)}dy=0.\end{array}$$

(4)

Therefore, we called $p_0\left(y\right)$ an entry–exit function.

3. Dynamical Properties of the Time-Delay System

Consider a system with a time delay $\tau$, which is a sufficiently small constant representing the delay of the prey in functional response.

$$\begin{array}{c}\frac{dx}{dt}=x(1-\frac{x}{K})-p(x(t-\tau ),y)y,\hfill \\ \frac{dy}{dt}=\epsilon y(p(x,y)-\alpha ).\hfill \end{array}$$

(5)

Using Taylor formula,

$$\begin{array}{c}x(t-\tau )\approx x\left(t\right)-x^{\prime }\left(t\right)\tau .\hfill \end{array}$$

(6)

Substituting (6) into System (5) yields

$$\begin{array}{c}\frac{dx}{dt}=(x(1-\frac{x}{K})-p(x,y)y)q\left(y\right),\hfill \\ \frac{dy}{dt}=\epsilon y(p(x,y)-\alpha ),\hfill \end{array}$$

(7)

where

$$q\left(y\right)=\left\{\begin{array}{ccc}& \frac{2+hy}{2+hy-\tau y},\hfill & \hfill x<2\\ & 1.\hfill & \hfill x\ge 2\end{array}\right.$$

3.1. Equilibria

Lemma 1.

For System (7), we have the following results:

(a) Systems (5) and (7) share the same equilibria. Therefore, System (7) has equilibria $E_0(0,0)$ and $E_1(K,0)$.

(b) If $h<\frac{8(1-\alpha )}{K}\le \frac{K(1-\alpha )}{K-2}$, System (7) has a unique positive equilibrium $E_{*}(x_{*},y_{*})$ with $x_{*}=\frac{\sqrt{K^2{(1-h)}^2+8\alpha kh}-K(1-h)}{2h},y_{*}=\frac{x_{*}-2\alpha }{\alpha h}$.

(c) The set $S=\left\{(x,y)|0<x<K,0<y<\frac{2-2\alpha }{\alpha h}\right\}$ is the invariant set of System (7).

Proof. 

Since the time delay does not affect the value of the equilibria of the system, Systems (5) and (7) have the same equilibria.

Let the right end of System (5) be 0, then

$$\begin{array}{cc}\hfill x\left(1-\frac{x}{K}\right)-p(x(t-\tau ),y)y& =0\hfill \\ \hfill \epsilon y\left(p\right(x,y)-\alpha )& =0\hfill \end{array}$$

Then, we have $E_0(0,0),E_1(K,0),E_{*}(x_{*},y_{*})$, where

$$\begin{array}{ccc}\hfill x_{*}& =\frac{\sqrt{K^2{(1-h)}^2+8\alpha kh}-K(1-h)}{2h},y_{*}\hfill & \hfill =\frac{x_{*}-2\alpha }{\alpha h}.\end{array}$$

By the definition of invariant set, we have (c). □

3.2. Hopf Bifurcation

We consider the bifurcation of System (7) at the unique positive equilibrium $E_{*}(x_{*},y_{*})$ about parameter $\tau$. For convenience, we denote

$$\begin{array}{ccc}\hfill & F(x,y)=(x(1-\frac{x}{K})-p(x,y)y)q\left(y\right),\hfill & \hfill G(x,y)=\epsilon y\left(p\right(x,y)-\alpha ).\end{array}$$

(8)

Thus, the linearized system of System (7) at the unique positive equilibrium $E_{*}(x_{*},y_{*})$ can be written as

$$\left(\begin{array}{c}{x}^{\prime }\hfill \\ y^{\prime }\hfill \end{array}\right)=\left(\begin{array}{cc}{F}_x\hfill & F_y\\ G_x\hfill & G_y\end{array}\right)\left(\begin{array}{c}x\hfill \\ y\hfill \end{array}\right).$$

(9)

Then, the characteristic equation of System (9) is

$$\begin{array}{c}{\lambda }^2-(F_x+G_y)\lambda -G_x{F}_y+F_x{G}_y=0.\end{array}$$

(10)

When $x\ge 2$, we have $G_x=0$. Therefore, Equation (10) has no pair of pure imaginary roots.

When $x<2$, we find

$$\begin{array}{cccc}\hfill & F_x=(1-\frac{2x_{*}}{K}-\frac{y_{*}}{2+h{y}_{*}})q\left(y_{*}\right),\hfill & \hfill & F_y=-x_{*}\frac{2q\left(y_{*}\right)+(2y_{*}+h{y}_{*}^2)q^{\prime }\left(y_{*}\right)}{{(2+h{y}_{*})}^2}+x_{*}(1-\frac{x_{*}}{K})q^{\prime }\left(y_{*}\right),\hfill \\ & G_x=\frac{\epsilon y_{*}}{2+h{y}_{*}},\hfill & \hfill & G_y=\epsilon (\frac{2{\alpha }^2}{x_{*}}-\alpha ).\hfill \end{array}$$

(11)

Lemma 2.

For Equation (10), we have

If $\tau =\frac{h{x}_{*}}{x_{*}-2\alpha }-\frac{h{x}_{*}^2}{\epsilon \alpha {(x_{*}-2\alpha )}^2}(1-\frac{2x_{*}}{K}-\frac{x_{*}-2\alpha }{h{x}_{*}})$, then Equation (10) has a pair of pure imaginary roots.

Proof. 

The condition guaranteeing Equation (10) has a pair of pure imaginary roots:

$$\begin{array}{c}{F}_x+G_y=0,G_x{F}_y-F_x{G}_y<0.\end{array}$$

(12)

Consider $\epsilon$ as a sufficiently small constant while $\tau$ is a small constant. We can easily obtain $G_x>0$, $F_y<0$, so $G_x{F}_y-F_x{G}_y<0$ holds.

According to the equation $F_x+G_y=0$, we obtain that

$$\begin{array}{c}(1-\frac{2x_{*}}{K}-\frac{y_{*}}{2+h{y}_{*}})q\left(y_{*}\right)+\epsilon (\frac{2{\alpha }^2}{x_{*}}-\alpha )=0.\end{array}$$

(13)

Replacing $y_{*}=\frac{x_{*}-2\alpha }{\alpha h}$ in Equation (13), we have

$$\begin{array}{c}(1-\frac{2x_{*}}{K}-\frac{\frac{x_{*}-2\alpha }{\alpha h}}{2+h\frac{x_{*}-2\alpha }{\alpha h}})q\left(\frac{x_{*}-2\alpha }{\alpha h}\right)+\epsilon (\frac{2{\alpha }^2}{x_{*}}-\alpha )=0,\end{array}$$

(14)

i.e.,

$$\begin{array}{c}\tau =\frac{h{x}_{*}}{x_{*}-2\alpha }-\frac{h{x}_{*}^2}{\epsilon \alpha {(x_{*}-2\alpha )}^2}(1-\frac{2x_{*}}{K}-\frac{x_{*}-2\alpha }{h{x}_{*}}).\end{array}$$

(15)

The proof is completed. □

Let $\lambda \left(\tau \right)=p\left(\tau \right)\pm iq\left(\tau \right)$ be the pair of pure imaginary roots of Equation (10) that satisfy $p\left(\tau \right)=0$, $q\left(\tau \right)\ne 0$. Then, the following transversality condition holds.

Lemma 3.

For Equation (10), we have ${\lambda }^{\prime }\left({\tau }^{*}\right)\ne 0$.

Proof. 

The condition for Equation (10) has a pair of pure imaginary roots:

$$\begin{array}{c}\hfill Re\left(\frac{d\lambda }{d{\tau }^{*}}\right)=\frac{d{F}_x}{2d{\tau }^{*}}=(1-\frac{2x_{*}}{K}-\frac{y_{*}}{2+h{y}_{*}})\frac{y_{*}(2+h{y}_{*})}{{(2+h{y}_{*}-{\tau }^{*}{y}_{*})}^2}\ne 0.\end{array}$$

Denote ${\tau }^{*}$$=\frac{h{x}_{*}}{x_{*}-2\alpha }-\frac{h{x}_{*}^2}{\epsilon \alpha {(x_{*}-2\alpha )}^2}(1-\frac{2x_{*}}{K}-\frac{x_{*}-2\alpha }{h{x}_{*}})$. Combining Lemmas 2 and 3, we conclude what follows.

Theorem 1.

For System (7), if $\tau ={\tau }^{*}$, then the equilibrium $E_{*}(x_{*},y_{*})$ is a Hopf bifurcation point and a limit cycle occurs.

3.3. Dynamics of Limit Systems

Consider the limit system of System (7). Letting $\epsilon \to 0$ in (7), we reach the fast subsystem

$$\begin{array}{cc}\hfill & \frac{dx}{dt}=x^{\prime }=(x(1-\frac{x}{K})-p(x,y)y)q\left(y\right),\hfill \\ \hfill & \frac{dy}{dt}=y^{\prime }=0.\hfill \end{array}$$

(16)

Considering the slow time scale $\tau =\epsilon t$ and letting $\epsilon \to 0$, we obtain a slow subsystem

$$\begin{array}{cc}\hfill & 0=(x(1-\frac{x}{K})-p(x,y)y)q\left(y\right),\hfill \\ \hfill & \frac{dy}{dt}=\dot{y}=y(p(x,y)-\alpha ),\hfill \end{array}$$

(17)

which is a differential-algebraic equation on the critical set $C_0=C_1\cup C_2\cup C_3$, with

$$\begin{array}{cc}\hfill & C_1=\left\{(x,y)|x=0,y\ge 0\right\},\hfill \\ \hfill & C_2=\left\{(x,y)|0<x<2,y=\frac{2(K-x)}{K-hK+hx}\right\},\hfill \\ \hfill & C_3=\left\{(x,y)|x\ge 2,y=\frac{2(K-x)x}{h{x}^2+K(2-hx)}\right\}.\hfill \end{array}$$

Assume that $0<h<1$ and $K>4$. Within this constraint, we find that the branch $C_3=\left\{(x,y)|x\ge 2,y=\frac{2(K-x)x}{h{x}^2+K(2-hx)}\right\}$ of $C_0$ has a unique generic fold point $D=(x_M,y_M)=(\frac{K}{2},\frac{2K}{8-hK})$. The other two branches $C_1$ and $C_2$ intersect at point $E_{*}(x,y)=(0,\frac{2}{1-h})$, and the two branches $C_2$ and $C_3$ intersect at point $F(x,y)=(2,\frac{2(K-2)}{K-hK+2h})$ (See Figure 2). Then, the critical set $C_0$ is divided into five parts by points D, E, and F:

$$\begin{array}{cc}\hfill & C_1^a=\left\{(x,y)|x=0,y>\frac{2}{1-h}\right\},C_3^a=\left\{(x,y)|x>\frac{K}{2},y=\frac{2(K-x)x}{h{x}^2+K(2-hx)}\right\},\hfill \\ \hfill & C_1^r=\left\{(x,y)|x=0,y<\frac{2}{1-h}\right\},C_3^r=\left\{(x,y)|2<x<\frac{K}{2},y=\frac{2(K-x)x}{h{x}^2+K(2-hx)}\right\},\hfill \\ \hfill & C_2=\left\{(x,y)|0<x<2,y=\frac{2(K-x)}{K-hK+hx}\right\},\hfill \end{array}$$

where $C_1^r$ and $C_3^r$ are normally hyperbolic repelling, $C_1^a$, $C_3^a$ and $C_2$ are normally hyperbolic attracting.

4. Relaxation Oscillation

We denote

$$\begin{array}{c}I\left(\overline{y}\right)={\int }_{\overline{y}}^{p_0\left(y\right)}\frac{f(0,y,0)}{g(0,y,0)}dy={\int }_{\overline{y}}^{p_0\left(y\right)}\frac{y-2-hy}{\alpha y(2+hy)}·\frac{2+hy}{2+hy-\tau y}dy.\end{array}$$

(18)

Lemma 4.

For System (7), there exists a unique ${\overline{y}}^{*}$($0<{\overline{y}}^{*}<\frac{2}{1-h}$) such that

$$\begin{array}{c}{\int }_{{\overline{y}}^{*}}^{y_M}\frac{f(0,y,0)}{g(0,y,0)}dy=0.\end{array}$$

(19)

Proof. 

Consider System (7), for $\overline{y}\in (0,\frac{2}{1-h})$,

$$\begin{array}{c}I\left(\overline{y}\right)={\int }_{\overline{y}}^{y_M}\frac{y-2-hy}{\alpha y(2+hy)}·\frac{2+hy}{2+hy-\tau y}dy.\end{array}$$

(20)

Therefore,

$$\begin{array}{c}I\left(\overline{y}\right)\to -\infty as\overline{y}\to 0.\end{array}$$

(21)

Furthermore,

$$\begin{array}{c}I\left(\frac{2}{1-h}\right)={\int }_{\frac{2}{1-h}}^{y_M}\frac{f(0,y,0)}{g(0,y,0)}dy={\int }_{\frac{2}{1-h}}^{y_M}\frac{y-2-hy}{\alpha y(2+hy)}·\frac{2+hy}{2+hy-\tau y}dy>0,\end{array}$$

(22)

and

$$\begin{array}{c}{I}^{\prime }\left(\overline{y}\right)=-\frac{y-2-hy}{\alpha y(2+hy)}·\frac{2+hy}{2+hy-\tau y}dy>0.\end{array}$$

(23)

Then, we combine (22) and (23) and conclude that there exists a unique ${\overline{y}}^{*}$($0<{\overline{y}}^{*}<\frac{2}{1-h}$) such that

$$\begin{array}{c}{\int }_{{\overline{y}}^{*}}^{y_M}\frac{f(0,y,0)}{g(0,y,0)}dy=0.\end{array}$$

(24)

Let us define $x_r$ to be the x-coordinate of the intersection point of $y={\overline{y}}^{*}$ and $C_3^a$. Define a singular slow–fast cycle ${\gamma }_0$, which include two slow segments on $C_3^a$ from $(x_r,{\overline{y}}^{*})$ to $(x_M,y_M)$ and on the positive y-axis from $(0,y_M)$ to $(0,{\overline{y}}^{*})$. Additionally, ${\gamma }_0$ contains two fast connections from $(x_M,y_M)$ to $(0,y_M)$ and $(0,{\overline{y}}^{*})$ to $(x_r,{\overline{y}}^{*})$, respectively. The next theorem shows the existence and uniqueness of stable relaxation oscillations. □

Theorem 2.

For System (7), assume $0<\epsilon \ll 1$, $0<\alpha <1$, $K>2$, $0<h<1$. Let V be a tubular neighborhood of ${\gamma }_0$. If $h<\frac{8(1-\alpha )}{K}$; then, for each fixed $\epsilon >0$ sufficient small, System (7) has a unique stable limit cycle ${\gamma }_{\epsilon }\subset V$, which is strongly attracting. Moreover, the cycle ${\gamma }_{\epsilon }$ is the unique stable limit cycle, which converges to ${\gamma }_0$.

Proof. 

According to Fenichel’s theory, the critical submanifold $C_3^a$ perturbs to the nearby slow manifold $C_{\epsilon }^a$, which is $O\left(\epsilon \right)$ near $C_3^a$. By Theorem 2.1 of [9] on the analysis of a jump point, the slow manifold $C_{\epsilon }^a$ jumps to another attracting branch $C_0^{+}$ when it reaches near the fold point $D(x_M,y_M)$.

Define two vertical sections ${\Delta }^{in}$ and ${\Delta }^{out}$ as shown in Figure 2:

$$\begin{array}{c}{\Delta }^{in}:=\left\{(2,y)\in V|y\in I_{in}\right\}and{\Delta }^{out}:=\left\{(2,y)\in V|y\in I_{out}\right\},\end{array}$$

(25)

where $I_{in}$ and $I_{out}$ are closed intervals centered at $y_M$ and $y_{out}=p_0\left(y_M\right)$, respectively. Through the flow of System (2), we define the transition map $\Pi :{\Delta }^{in}\to {\Delta }^{in}$, which is a composition of the next two maps

$$\begin{array}{c}{\Pi }_1:{\Delta }^{in}\to {\Delta }^{out}and{\Pi }_2:{\Delta }^{out}\to {\Delta }^{in}.\end{array}$$

(26)

Then, $\Pi :{\Delta }^{in}\to {\Delta }^{in}$ is given by the composition $\Pi :={\Pi }_2\circ {\Pi }_1$.

Now, we analyze the properties of these two maps ${\Pi }_1$ and ${\Pi }_2$.

(a)

Analysis of ${\Pi }_1$. Using the same proof as in Lemma 4, for each $y_0\in I_{in}$, we can define $p_0\left(y_0\right)$, with $0<p_0\left(y_0\right)<\frac{2}{1-h}$ through the formula

$$\begin{array}{c}{\int }_{y_0}^{p_0\left(y_0\right)}\frac{f(0,y,0)}{g(0,y,0)}dy=0.\end{array}$$

(27)

(b)

Analysis of ${\Pi }_2$. System (7) is continuous at line $x=2$. Consider two orbits ${\gamma }_{\epsilon }^1$ and ${\gamma }_{\epsilon }^2$ that start on ${\Delta }^{out}$. By Fenichel’s theory, ${\gamma }_{\epsilon }^1$ and ${\gamma }_{\epsilon }^2$ will be attracted to $C_{\epsilon }^a$ at the exponential rate $O\left(e^{\frac{-1}{\epsilon }}\right)$. Based on Theorem 2.1 of [9], ${\gamma }_{\epsilon }^1$ and ${\gamma }_{\epsilon }^2$ pass by the generic fold point $D(x_M,y_M)$ and then exponentially are contracted toward each other. After that, they fly to ${\Delta }^{in}$.

Therefore, according to the results of (a) and (b), we see that the transition map $\Pi :{\Delta }^{in}\to {\Delta }^{in}$ is a contraction at the exponential rate $O\left(e^{\frac{-1}{\epsilon }}\right)$. Furthermore, we obtain that $\Pi$ has a unique fixed point in ${\Delta }^{in}$ by the contraction mapping theorem, which must be stable. This fixed point gives a unique stable relaxation oscillation cycle ${\gamma }_{\epsilon }\subset V$ of System (7) passing through ${\Delta }^{in}$ for each $0<\epsilon \ll 1$.

By Fenichel’s theory and theorem 2.1 of [9], we obtain that the relaxation oscillation cycle ${\gamma }_{\epsilon }$ converges to the slow–fast cycle ${\gamma }_0$ as $\epsilon \to 0$.

Therefore, the relaxation oscillation cycle ${\gamma }_{\epsilon }$ is a unique stable limit cycle of System (7) located in V for each $0<\epsilon \ll 1$.

This completes the proof. □

5. Numerical Simulation

In this section, we present an example to illustrate Theorem 2.

Set the parameter values $(h,\alpha ,\epsilon ,K,\tau )=(0.3,0.55,0.01,10,0.005)$, and it is easy to show that System (7) has a unique positive equilibrium point $E_{*}(1.48,2.29)$. Numerical simulation shows us that it has a unique limit cycle. Figure 3 and Figure 4 show the phase portrait and the time series of the relaxation oscillation cycle, respectively.

In addition, we compare the phase portraits and time series of System (7) with System (5). Figure 5 shows the phase portraits of System (7) and System (5). It can be seen that the two curves are very close to each other. Figure 6 and Figure 7 are the time series of System (7) and System (5), respectively.

6. Conclusions

In this paper, we found that the geometric singular perturbation theory is effective for the study of slow–fast predator–prey models with piecewise smooth response functions. First, the condition that Equation (10) has a pair of pure imaginary roots is given using linearization analysis method, and the transversal condition of Hopf bifurcation condition is verified. Lemmas 1–3 and Theorem 1 are then derived. The existence and uniqueness of relaxation oscillations are proven using geometric singular perturbation theory and the entry–exit function, as shown by Theorem 2. Within the time delay, the dynamic characteristics of the system may change, so the influence of time delay cannot be ignored in this model. Finally, numerical simulation is carried out with fixed parameter values, and the results are consistent with the theoretical analysis.