Relaxation oscillations of a piecewise-smooth slow-fast Bazykin’s model with Holling type I functional response

: In this paper, we consider the dynamics of a slow-fast Bazykin’s model with piecewise-smooth Holling type I functional response. We show that the model has Saddle-node bifurcation and Boundary equilibrium bifurcation. Furthermore, it is also proven that the model has a homoclinic cycle, a heteroclinic cycle or two relaxation oscillation cycles for di ff erent parameters conditions. These results imply the dynamical behavior of the model is sensitive to the predator competition rate and the initial densities of prey and predators. In order to support the theoretical analysis, we present some phase portraits corresponding to di ff erent values of parameters by numerical simulation. These phase portraits include two relaxation oscillation cycles, an unstable relaxation oscillation cycle surrounded by a stable homoclinic cycle; the coexistence of a heteroclinic cycle and an unstable relaxation oscillation cycle. These results reveal far richer and much more complex dynamics compared to the model without di ff erent time scale or with smooth Holling type I functional response.


Introduction
The study of the long term symbiosis and complex dynamics among different interacting species is a hot topic and has attracted more and more attention of mathematicians and biologists over decades.Under the observations of interactions between various species, researchers purposed many appropriated mathematical models.Among these models, the Lotka-Volterra predator-prey model proposed by Lotka [1] and Volterra [2] is a widely known model due to its universal application and significance.Now, the various Lotka-Volterra models are developed to study the complex population dynamics in different ecological situations.Hence, taking into consideration prey competition and different functional responses, the original Lotka-Volterra model is extended to the generalized Gause type case [3,4] du dt = au − bu 2 − vp(u), where u and v separately represent the amount of prey and predators; a and b separately represent the intrinsic growth rate and compete rate of prey; c and d separately represent the death rate and maximal growth rate of predators; p(u) is the functional response describes the variation of the amount of prey affected by the attacks of predators.However, it is worth noticing that not only is there competition among prey, but predators also compete each other for the limited resource such as the size of the habitat to live and reproduce [5].Hence, by introducing predator competition, system (1.1) can be rewritten as where h is the competition rate of predators.This is the well-known Bazykin's predator-prey model proposed in [6].The functional response p(u) is usually classified into four types firstly purposed by the biologist Holling [7,8], for example the Holling Type I functional response where m is a maximal per capita consumption rate and α is the half-saturation constant-the prey's number at which the per capita consumption rate is half of its maximum m, and other three Holling Type functional responses, see Table 1.Note that these functional response are all bounded function which is suitable for actual field data and the Holling type I functional response (1.3) is continuous but not smooth, which is also called the piecewise-smooth Holling type I functional response.
Table 1.Holling types II, III and IV functional responses and their generalizations.

Holling type Definition
Recent researches related to the Bazykin's model with different functional responses obtain some interesting complex dynamics and bifurcations.For example, Bzaykin [9] showed that there is a threshold value c 1 such that system (1.2) with the functional response p(u) = mu has the globally asymptotically stable boundary equilibrium ( a b , 0) if c > c 1 and the unique globally asymptotically stable positive equilibrium if c < c 1 ; many researchers [5,6,[9][10][11][12][13][14][15][16][17] including Bazykin et al. [5,6,9,10], Hainzl et al. [11,12] and Lu and Huang [13] investigated system (1.2) with Holling type II functional response by the theoretical analysis or numerical methods and got complex dynamics and rich bifurcation phenomenons such as the location and stability of equilibriums, the existence of limit cycles, the degenerate Bogdanov-Takens bifurcation, the Hopf bifurcation, etc.
In this paper, we investigate Bazykin's predator-prey model (1.2) with piecewise smooth Holling type I functional response (1.3).Before going into details, we apply the following rescaling transformation and removing the bar notation, then the system (1.2) can be rewrote as the following non-dimensional system where u and v are the non-dimensional variables; b = bd a , ϵ = dm a , h = haα dm 2 and c = c dm are nondimensional parameters.Furthermore, we assume the growth rate of predators is much smaller than that of prey which means that system (1.4) is a piecewise-smooth slow-fast system with the small parameter 0 < ϵ ≪ 1.Note that our assumption is reasonable because species at different trophic levels have different growth rates, which may vary several orders of magnitude, and the growth time of individuals increases gradually from the bottom of the food chain to the top [18].Moreover, many observations of interactions between prey and predators such as hares and lynx [19], phytoplankton and zooplankton [20], insects and birds [21], etc. indicate the prey grow much faster than predators.
So far, there are few studies on the slow-fast Bazykin's model and these studies mainly focus on the dynamics of the model (1.2) with smooth Holling II functional response [15,16] such as relaxation oscillation cycles, canard phenomenon and so on.Hence, Our main aim in this paper is to study the dynamics of the slow-fast Bazykin's model with a piecewise-smooth functional response (1.4) by using the geometrical singular perturbation theory [22][23][24][25][26][27][28][29] and piecewise-smooth dynamical system theory [30][31][32][33][34][35].For the various values of parameters, we will show that there are the coexistence of two relaxation oscillation cycles with different stability, an unstable relaxation oscillation cycle surrounded by a stable homoclinic cycle and a heteroclinic cycle enclosing an unstable relaxation oscillation cycle.Furthermore, the system (1.4) undergoes a series of bifurcations such as Saddle-node bifurcation and Discontinuous saddle-node bifurcation of codimension 1.Moreover, we also give some conditions of the global stability of the unique positive equilibrium.Numerical simulations with the help of the "PPlane8" tool of Matlab [36] are presented to illustrate the theoretical results.
The rest of this paper is organized as follows.In Section 2, we study the critical manifold and the existence and types of equilibriums of system (1.4).In Section 3, we show that the dynamics and bifurcations of system (1.4) such as relaxation oscillation cycles, homoclinic cycle, heteroclinic cycle, etc.A brief discussion is given in the last section.

The critical manifold and equilibriums
Based on the ecological viewpoint, we are keen on the dynamics of the system (1.4) in the first quadrant R 2 + = {(u, v) | u ≥ 0, v ≥ 0} and the switching boundary Σ = {(u, v) | u = 1} splits the first quadrant into two regions denoted by The system (1.4) is smooth in Σ (∓) and determined by the following systems and Since the vector field is locally Lipschitz, the fundamental existence and uniqueness theory is true for the system (1.4) [31].Note that the trajectories of the system (1.4) transversally pass through the switching boundary Σ.Moreover, in order to guarantee the density of preys can support the growth of predators, we assume throughout the paper that 0 < b < 1 4 and 0 < c < 1.These are also the necessary condition for the existence of relaxation oscillation cycles.Hence, we will study the system (1.4) in the parametric region Before going into the detail dynamics, we firstly show some basic properties of the system (1.4).
Lemma 2.1.The system (1.4) has the invariant set and it also has no limit cycle which is entirely lied in the region Σ (−) .
Proof.It is clear that u = 0 and v = 0 are two invariant straight lines of system (1.4) and the boundary of set Ω is b , which means that the vector field of the system (1.4) in the line ∂Ω never points outside.Hence, any trajectories of the system (1.4) are confined in the set Ω as they enter Ω in finite time.
For (u, v) ∈ Σ (−) , we construct the Dulac function φ(u, v) = 1 uv and have which implies the system (1.4) has no limit cycles in the region Σ (−) because of the Dulac's criteria.□ Hence, if the system (1.4) has a limit cycle, it will either entirely locate in the region Σ (+) or consist of trajectories in the regions Σ (+) and Σ (−) and transversally pass through the switching boundary Σ.

Critical manifold
Under the time scale transformation τ = ϵt, we can get the equivalent system The systems (1.4) and (2.1) are individually known as the fast system and the slow system because of their different time scales.By setting ϵ = 0 in systems (1.4) and (2.1), we get the degenerate system and the layer system Hence, the critical manifold is , which is composed of the singular points of layer system (2.3) and can be divided into three parts as follows, see Figure 1, Note that the sub-manifold S 2 0 is a normally hyperbolic attracting sub-manifold and the non-normally hyperbolic points A(0, 1) and M( 12b , 1 4b ) separately split the sub-manifolds S 1 0 and S 3 0 into the normally hyperbolic attracting sub-manifolds and normally hyperbolic repelling sub-manifolds For 0 < ϵ ≪ 1, the Fenichel theory [22,23] indicates the normally hyperbolic attracting submanifolds S 2 0 , S 1a 0 and S 3a 0 can be perturbed to the attracting slow manifolds S 2 ϵ , S 1a ϵ and S 3a ϵ and the normally hyperbolic repelling sub-manifolds S 1r 0 and S 3r 0 can be perturbed to the repelling slow manifolds S 1r ϵ and S 3r ϵ .Hence, the trajectory starting near S 1a 0 is attracted to the v-axis, then it moves down with O(ϵ) speed.It is clear that the trajectory passes the non-hyperbolic point A(0, 1) and leaves the O(ϵ) neighborhood of S 1r 0 at the point (0, p ϵ (v)) satisfying lim ϵ→0 p ϵ (v) = p 0 (v).The function p 0 (v) is called the entry-exit function [28,29] and satisfies the following lemma.Lemma 2.2.For system (1.4) and v 0 ∈ (1, +∞), there is a unique p 0 (v 0 ) ∈ (0, 1) satisfying p 0 (v 0 ) Proof.Let it is easy to verify that I(1) < 0 and lim ṽ→0 I(ṽ) = +∞.Since I ′ (ṽ) = ṽ − 1 ṽ(hṽ + c) < 0, ṽ ∈ (0, 1), there is a unique ṽ * ∈ (0, 1) such that I(v * ) = 0 which implies p 0 (v 0 ) = v * .□ Note that the above proof indicates the entry-exit function p 0 (v 0 ) is monotone decreasing function in (1, +∞).

Equilibriums and their types
In this subsection, we mainly study the existence and stability of equilibriums of system (1.4).We can get the equilibriums by a straight calculation of the following equations 1.The critical manifold and equilibriums of system (1.4) with different parameters in region D. and the dynamics of system (1.4) near these equilibriums are determined by the Jacobian matrix of system (1.4) at these equilibriums, which is given by and here (u * , v * ) is the coordinate of these equilibriums.Now, set we give the following lemmas about the existence and stability of equilibriums.
Proof.The number and locations of equilibriums of system (1.4) can be obtained by the straight calculation of Eq (2.6).Next, we determine the types of equilibriums O(0, 0), B( 1 b , 0) and E 1 (u 1 , v 1 ).For equilibrium O(0, 0), the eigenvalues of the Jacobian matrix J (−) are λ 1 = 1 and λ 2 = −ϵc < 0 which implies the equilibrium O(0, 0) is a saddle point.Similarly, the eigenvalues of the Jacobian matrix 0 , the similar calculation gives the eigenvalues of the Jacobian matrix J (+) as 0 is a stable node.For equilibrium E 1 (u 1 , v 1 ) ∈ S 2 0 , we calculate the determinant and trace of Jacobian matrix J (−) as which implies E 1 (u 1 , v 1 ) ∈ S 2 0 is stable.Next, we will prove the global stability of equilibrium E 1 when h > h 2 .We construct the trapping region and assert that the vector field of system (1.4) point inside at the boundary ∂Ω 1 .So the trajectories cannot leave after entering the region Ω 1 .Furthermore, we construct the same Dulac function in Lemma 2.1 as φ(u, v) = 1 uv and have Hence, based on the Dulac's criteria, the system (1.4) has no limit cycles in the region Ω 1 , which indicate the equilibrium E 1 is globally stable in Ω 1 .Moreover, the trajectories starting in the region Σ\Ω 1 enter into the region Ω in finite time because of the Fenichel theory and Lemma 2.2.Then we can assert that system (1.4) has a globally stable node E 1 when h > h 2 .□ Lemma 2.4.If h = h 1 or h = h 2 , then system (1.4) has two positive equilibriums.More precisely, 1) If h = h 1 , then system (1.4) has a stable equilibrium E 1 (u 1 , v 1 ) located in S 2 0 and an attracting saddle-node E M ( 1 2b , 1 4b ) which is the vertex point in S 3 0 , see Figure 1(b).2) If h = h 2 , then system (1.4) has a stable equilibrium E 1 (u 1 , v 1 ) located in S 2 0 and a boundary equilibrium E B (1, 1 − b) located in the switching boundary Σ, see Figures 1(d) and 2(b).
Proof.The existence and location of two positive admissible equilibriums follows from the straight calculation of Eq (2.6).By some simple calculations, we get and the eigenvalues of equilibrium E M ( 1 2b , 1 4b ) are λ 1 = 0 and λ 2 = −ϵ(1−c) < 0. Hence, the equilibrium E 1 is stable and the equilibrium E M is an attracting saddle-node which can be proven based on the conditions in [37].
) is a saddle in S 3r 0 and E 3 (u 3 , v 3 ) is a stable node in S 3a 0 , see Figure 1(c).The proof of Lemma 2.5 is similar to that given in Lemmas 2.3 and 2.4 and so is omitted.

Dynamics of the slow-fast system
In this section, we are keen on various possible dynamics of the slow-fast system (1.4) including relaxation oscillation, saddle-node bifurcation, boundary equilibrium bifurcation and so on.

Saddle-node bifurcation
According to Lemmas 2.3-2.5, we know that is a saddle-node bifurcation surface.When parameters change from one side of the surface S N to another one, the number of positive equilibriums of system (1.4) in Σ (+) changes from zero to two.So there is a critical competition rate h 1 such that the two species coexist in the form of a positive equilibrium for suitable choices of initial values when In what follows, we show that system (1.4) has homoclinic orbit, heteroclinic orbit and relaxation oscillation cycles in both sides of saddle-node bifurcation surface S N.
First, we investigate the relaxation oscillation cycles of system (1.4) when h < h 1 .Now, we denote the function Î(ṽ) as the function (2.5) with v 0 = 1 4e and state the main theorem as follows Theorem 3.1.
then system (1.4) has a hyperbolically stable relaxation oscillation cycle Γ 1 ϵ and a hyperbolically unstable relaxation oscillation cycle γ ϵ which separately converges to in the Hausdorff distance as ϵ → 0, see Figures 3(a) and 4(a).
Proof.Our first goal is to prove the existence of the stable relaxation oscillation cycle Γ 1 ϵ .To begin with, it is easy to see that the vertex point M( 1 2b , 1 4b ) is a general fold and there exists v * ∈ (0, 1 − b) such that p 0 ( 1 4b ) = v * based on Î(1 − b) < 0 and the properties of function Î(ṽ).Hence, we construct the singular slow-fast cycle Γ 1 0 , see Figure 3(a), as which separately are vertical region of ⌢ MK 1 and ⌢ Q 1 D 1 and δ is a small positive parameter.We construct the Poincaré map Π as For the trajectory starting in ∆ 0 , it will arrive at the neighborhood of critical manifold S 1a 0 and move downwards until reach the neighborhood of Q 1 .The trajectory will move along the layers of system (2.3) and pass through the region ∆ 1 before reaching near the critical manifold S 3a 0 because system (1.4) is continuous on the switching boundary Σ.Then the trajectory will move along S 3a 0 and jump into ∆ 0 in the neighborhood of fold point M. Hence, the relaxation oscillation cycle Γ 1 ϵ is equivalent to the fixed point of Poincaré map Π which is ) by Lemma 2.2 and Theorem 2.1 in [26].According to the Contraction Mapping Theorem, there is an unique fixed point of Π corresponding to the relaxation oscillation cycle Γ 1 ϵ .Furthermore, it is clear that Γ 1 ϵ is a hyperbolically stable limit cycle and converges to the singular slow-fast cycle Γ 1 0 as ϵ → 0. Next we prove the existence of the unstable relaxation oscillation cycle γ ϵ .Since E 1 (u 1 , v 1 ) is stable and the relaxation oscillation cycle Γ 1 ϵ is hyperbolically stable, we can conclude that there exists at least an unstable limit cycle inside Γ ϵ based on Poincaré-Bendixson Theorem.According to the Geometric Singular Perturbation Theory [23], the limit cycles of slow-fast system (1.4) are usually perturbed by singular slow-fast cycles.Clearly, system (1.4) has an unique singular slow-fast cycle γ 0 inside the relaxation oscillation cycle Γ 1 ϵ for ϵ → 0, see Figure 3(a), which are constructed as ) and K 2 (0, v * ).Note that p 0 (v * ) = 1 − b and ū * is the smaller root of equation H 2 (u) = v * .Hence, if the limit cycles exist, they must be in the neighborhood of γ 0 .By the Corollary 4.3 in [38], the sign of the following integral determines the limit cycle is hyperbolically unstable if it exists.Thus, we can claim that there is an unique relaxation oscillation cycle γ ϵ near γ 0 because two adjacent limit cycles don't have the same stabilities.□ Combined with the above analysis and applying numerical simulations by the "PPlane8" tool in Matlab, we give the phase portrait of system (1.4) when b = 0.2, c = 0.5 and h = 0.3, see Figure 4(a).Note that system (1.4) has a hyperbolically unstable relaxation oscillation cycle γ ϵ surrounded by a hyperbolically stable relaxation oscillation cycle Γ 1 ϵ .Second, we investigate the existence of a homoclinic cycle of system (1.4) with h = h 1 and have the following theorem.Theorem 3.2.When 0 < b < 1 4 , 0 < c < 1 and h = h 1 , system (1.4) has a stable equilibrium E 1 (u 1 , v 1 ) and a attracting saddle-node E M ( 1 2b , 1 4b ).Moreover, if Î(1 − b) < 0, then system (1.4) has a stable homoclinic cycles Γ 2 ϵ and a unstable relaxation oscillation cycle γ ϵ which separately converge to in Hausdorff distance as ϵ → 0, see Figures 3(b) and 4(b).
Proof.For ϵ = 0, it is clear that there is a singular orbit Γ 2 0 in saddle-node E M , which can be constructed by the same way of Γ 1 0 in Theorem 3.1, see Figure 3(b).Furthermore, the layer of system (2.3) at E M will be perturbed to one dimensional unstable manifold W u ϵ (E M ) of system (1.4) at E M if 0 < ϵ ≪ 1.By Lemma 2.2, the unstable manifold W u ϵ (E M ) will be attracted to the neighborhood of S 1a 0 and move down until it reach near Q 1 .Then W u ϵ (E M ) move quickly from the neighborhood of S 1r 0 to the neighborhood of S 3a 0 .Next, we show that the unstable manifold W u ϵ (E M ) tends to saddle-node E M .We construct the trapping region Ω 1 which is same to (2.7) in Lemma 2.4 and the orbit cannot leave if it enters the region Ω 1 .Since the unstable manifold W u ϵ (E M ) enter this region Ω 1 and the orbits is monotone increasing in v, the unstable manifold W u ϵ (E M ) tends to saddle-node E M along its stable manifold W s ϵ (E M ) which is perturbed by W s (S 3a 0 ) based on the Fenichel theorem.Hence, when 0 < ϵ ≪ 1, there is a homoclinic cycle Γ 2 ϵ which converges to singular orbit Γ 2 0 as ϵ → 0. Based on the Theorem 5 in [37], the stability of the homoclinic cycle Γ 2 ϵ is determined by the sign of the following integral Hence, the homoclinic cycle Γ 2 ϵ is stable.The proof of existence and stability of the relaxation oscillation cycle γ ϵ is omitted because it is similar to that in Theorem 3.1.□ According to the analysis of Theorem 3.2 and using numerical simulation by the "PPlane8" tool in Matlab, the phase portrait of system (1.4) with b = 0.2, c = 0.5 and h = 0.4 is presented in Figure 4(b).Note that system (1.4) has a stable homoclinic cycle Γ 2 ϵ enclosing a hyperbolically unstable relaxation oscillation cycle γ ϵ .
Third, we study the existence of a heterocilinic cycle of system (1.4) when h > h 2 .Similarly, we denote Ĩ(ũ) as the function (2.5) with v 0 = 1−c h .By applying the similar analysis methods in Theorem 3.2 and Lemma 1 in [39], we can derive the following theorem.
then system (1.4) has a heteroclinic cycle Γ 3 ϵ and a unstable relaxation oscillation cycle γ ϵ which separately converge to in Hausdorff distance as ϵ → 0, see Figures 3(c) and 4(c).
According to the above analysis and some numerical simulations by the "PPlane8" tool in Matlab, we give the phase portrait of system (1.4) with b = 0.2, c = 0.5 and h = 0.45, see Figure 4(c).Note that the two unstable manifolds of saddle E 2 (u 2 , v 2 ) tend to the stable node E 3 (u 3 , v 3 ) and form the heteroclinic cycle.From Figure 4, we know that the number of equilibriums of system (1.4) changes from one to three when parameters pass through the saddle-node bifurcation surface.Furthermore, under some conditions, system (1.4) always has a small relaxation oscillation cycle γ ϵ which is separately surrounded by a big relaxation oscillation cycle, homoclinic cycle and heteroclinic cycle for different values of parameter h in the neighborhood of h 1 .

Boundary equilibrium bifurcation
In this subsection, we will discuss the boundary equilibrium bifurcation of system (1.4).From Lemma 2.4, we know that the boundary equilibrium E B (1, 1 − b) is an admissible equilibrium for both systems in Σ (−) and Σ (+) , which is a mark of boundary equilibrium bifurcation of codimension 1.Hence, incorporating Lemmas 2.3 and 2.5, we present the discontinuous saddle-node bifurcation surface as follows When the parameters vary from one side of the surface to another, system (1.4) has two positive equilibriums E 1 (u 1 , v 1 ) and E 2 (u 2 , v 2 ) which are separately located in regions Σ (−) and Σ (+) .Then the discontinuous saddle-node bifurcation gets two positive equilibriums.

Discussion
In this paper, we consider the slow-fast Bazykin's predator-prey model with piecewise smooth Holling I functional response.Our qualitative analysis on system (1.4) reveals that system has complex dynamics and bifurcations and the competition rate h plays a critical role which affects not only the number and type of equilibriums but also the type of bifurcations in this model.When the values of parameters vary, system (1.4) undergoes the saddle-node bifurcation and the discontinuous saddle-node bifurcation.For different values of parameters, it is clear that system (1.4) has a hyperbolically unstable relaxation oscillation cycle which is respectively surrounded by a hyperbolically stable relaxation oscillation cycle, a homoclinic cycle and a heteroclinic cycle.These complex dynamics cannot occur in the system (1.4) with single time scale and smooth Holling type I functional response.In fact, the system (1.4) with ϵ > 1 and smooth Holling type I functional response [9] has at least one positive equilibrium which is globally stable if c < c 1 , here c 1 is a threshold value.Hence, compared with our result in this paper, it is clear that the different time scale and piecewise-smooth functional response in system (1.4) can introduce more complex dynamical behaviour.
These complex dynamical phenomenons of system (1.4) show that the complexity of dynamical behaviors of the Bazykin's model.For the biological interpretations of these complex dynamical