Dynamic analysis of a Leslie–Gower-type predator–prey system with the fear effect and ratio-dependent Holling III functional response

. In this paper, we extend a Leslie–Gower-type predator–prey system with ratio-dependent Holling III functional response considering the cost of antipredator defence due to fear. We study the impact of the fear effect on the model, and we ﬁnd that many interesting dynamical properties of the model can occur when the fear effect is present. Firstly, the relationship between the fear coefﬁcient K and the positive equilibrium point is introduced. Meanwhile, the existence of the Turing instability, the Hopf bifurcation, and the Turing–Hopf bifurcation are analyzed by some key bifurcation parameters. Next, a normal form for the Turing–Hopf bifurcation is calculated. Finally, numerical simulations are carried out to corroborate our theoretical results.


Introduction
As one of the most common mutual relationships between two populations in nature, predator-prey relationship plays a significant role in ecology and mathematical biology [2]. Since the dynamic behaviours of predator-prey models were formulated by Lotka and Volterra, many experts and scholars have studied various types of predator-prey models over the past decades, thus establishing the theoretical basis for interspecific interactions [1,4,12,14]. In 1948, Leslie and Gower introduced a functional response called Leslie-Gower, which was based on the Logistic model proposed by Verhulst. This functional response can be adapted to describe the correlation between a reduction in the number of predators and the per capita availability of the predator's preferred food [2,17]. The Leslie-Gower predator-prey system typically takes the following form: x(0) > 0, y(0) > 0, r, s, K, p > 0, where p is the prey-to-predator conversion factor, and the term y/(px) is known as the Leslie-Gower term, which means the scarcity of prey leads to a loss of predator density. The functional response is an important component of the predator-prey model describing the way in which predator-prey interactions occur. In proposing this function, an essential role is played by the behavior of prey and predator [17,18]. Some common types of functional responses are listed here. For example, Holling I-III functional response [20], Beddington-DeAngelis functional response [10], Leslie-Gower functional response [11], Crowley-Martin functional response [9], and ratio-dependent functional response [6].
The ratio-dependent functional response, a predator-dependent functional response, is an important form in characterizing the functional responses of predator and has a better modeling for a predator-prey system [25]. Based on the existed studies and the above considerations, many experts and scholars have replaced f (x) in (1) with a ratio-dependent functional response and considered its spatial distribution patterns and dispersal mechanisms. Making a suitable nondimensional scaling, system (1) reduces to ∂u(x, t) ∂t where u(x, t) and v(x, t) denote, respectively, the densities of the prey and the predator species at the space location x and time t. ∆ is the Laplacian operator. The positive constants d 1 and d 2 represent separately the diffusion coefficients of the prey and the predator. Parameters r, m, and β are all positive constants. The nonlinear term βu 2 /(u 2 + mv 2 ) in (2) is called ratio-dependent Holling type III functional response [5,21], and it indicates the growth rate of predator per capita, which is a function of the ratio of the number of prey to predator. The term v/u is known as the Leslie-Gower term [8], and it means the scarcity of prey leads to a loss of predator density [16]. There have been several reports on the dynamic behavior of this ratio-dependent predator-prey system (2). The Turing and Turing-Hopf bifurcations of system (2) were investigated and the spatiotemporal patterns of system (2) in two-dimensional space were discovered by Chen and Wu [3]. In [16], under homogeneous Neumann boundary conditions, the existence conditions of Hopf bifurcation, Turing instability of spatial uniformity and Turing-Hopf bifurcation in onedimensional space were shown. In Chang and Zhang [2], they used the Leray-Schauder degree theory and the implicit function theorem to report the nature of the spatially homogeneous Hopf bifurcation and the nonexistence and existence of nonconstant steady states.
In proposing model (2), the fear effect of prey for predator has been neglected. All prey respond to the risk of predation and then exhibit a variety of antipredator responses. The common ones are habitat change, foraging, vigilance, and other different physiological changes [15,19,[22][23][24]. Some previous experimental studies have confirmed this result. In addition, the fear of predation has immense impact on prey species. The mortality of prey is directly increased by predation. Moreover, the fear of prey may affect the physiological condition of juvenile prey, which is harmful to its survival in adulthood, for example, reducing the reproduction of prey [19,22]. Furthermore, strong behavioral responses can be elicited by the fear of predation enough to affect the population and life history of entire prey populations [13]. According to their ideas [22][23][24], we extend model (2) by multiplying the production term by a factor f (K, v) = 1/(1 + Kv) taking in account the cost of antipredator defence due to fear [22], then model (2) becomes where parameter K > 0 refers to the level of fear. The organizational structure of this paper is as follows. In Section 2, both the existence of positive equilibrium state and the linear stability of system (3) are analysed. Furthermore, the bifurcation analysis is investigated in Section 3, where the existence of the Turing bifurcation, the Hopf bifurcation, and the Turing-Hopf bifurcation are shown by choosing the d 2 , r, and K as bifurcation parameters, respectively. Moreover, the normal form of the Turing-Hopf bifurcation with r and d 2 as parameters for system (3) near the unique positive constant equilibrium is obtained in Section 4. Finally, in Section 5, numerical simulations are carried out to verify the obtained theoretical conclusions.

Existence and stability of equilibria
Firstly, the existence of coexisting equilibrium is analyzed by considering the following equation: Obviously, (1, 0) is a boundary equilibrium point of system (3). In this article, we mainly study the relevant properties of the coexistence equilibrium of system (3). By using the second equation of (4), we are able to obtain u = v. Substituting this into the first equation of (4), we obtain the expression h(u) = −K(1 + m)u 2 − (1 + m + Kβ)u + (1 + m − β). Theorem 1. The following statements are correct: (i) Assume that 1 + m − β > 0 holds, then h(u) = 0 must have a positive real solution. (ii) The equilibrium of prey u * is a strictly decreasing function with respect to K when 1 + m − β > 0.
(ii) In the case with the fear factor, i.e., K > 0, for system (3), when 1 + m − β > 0, we can obtain the equilibrium of prey Simple computation shows that which indicates that u * is a strictly decreasing function with respect to K. That is, increasing the level of K can decrease the value of the prey equilibrium u * .
Theorem 1 means that system (3) has a positive equilibrium point. We suppose that (u * , v * ) is the positive equilibrium point of system (3) for the remainder of the article. Furthermore, it also means that if the model has a unique positive equilibrium, the fear coefficient K has an effect on the positive equilibrium point: as the fear coefficient K increases, the positive equilibrium point of prey u * gradually decreases (see Fig. 1). Next, we conduct a linear stability analysis of system (3). For the positive equilibrium (u * , v * ), the linearized system of (3) is where and It is obvious to know that the characteristic equation of the linearized system (5) is where µ n = n 2 l 2 , TR n = −(d 1 + d 2 )µ n − r + a 1 , and the eigenvalues of system (3) are given by Then we make the following hypotheses: (A1) a 1 − r < 0; (A2) a 1 + a 2 < 0.
If Assumptions (A1) and (A2) are both valid, then when n = 0, there is TR 0 < 0 and DET 0 > 0, that is to say, the real parts of the eigenvalues of system (3) are all less than zero. Therefore, the following theorem is obtained.
Theorem 2. Suppose that (A1) and (A2) hold. Then the ordinary differential equation system corresponding to system (3) is locally asymptotically stable at the positive equilibrium point (u * , v * ).

Bifurcation analysis 3.1 Turing instability
In this section, the existence conditions of the Turing instability are analyzed. Under the assumptions that (A1) and (A2) are established, it is known from Theorem 2 that there is TR n < TR 0 < 0 for n ∈ N 0 . Then d 2 (see (7)) is selected as the Turing bifurcation line parameter. Let a 1 > 0 and discuss in the following three situations: After the analysis and discussion of the above three situations, we can get the following theorem.
Theorem 3. Suppose (A1) and (A2) hold. Then the following statements are correct: (i) In Case 1 or Case 2, the positive equilibrium point (u * , v * ) is locally asymptotically stable. (ii) In Case 3, if there does not exist a µ n (n ∈ N 0 ) such that DET n < 0, then the positive equilibrium point (u * , v * ) of system (3) is locally asymptotically stable; conversely, if there exists at least a µ n (n ∈ N 0 ) such that DET n < 0, then the positive equilibrium point (u * , v * ) of system (3) is Turing unstable.
Proof. Suppose that (A1) and (A2) hold. Under the condition that the parameters of Case 1 or Case 2 are satisfied, there is TR n < TR 0 < 0 for n ∈ N 0 , then if DET n > 0 (n ∈ N 0 ), system (3) has the eigenvalue of the negative real part. Then statement (i) is proved. When the parameter relationship belongs to Case 3 and there is not n ∈ N 0 such that DET n < 0, then a similar method can be used to prove the conclusion. When the parameter relationship belongs to Case 3 and there is a n 1 ∈ N 0 such that DET n 1 < 0, then the real part of the eigenvalue λ (n 1 ) 1 = (TR n 1 + TR 2 n 1 −4 DET n 1 )/2 of system (3) will be positive, which means that the positive equilibrium point (u * , v * ) of system (3) becomes no longer stable. Then the statement (ii) is proved.
Next, we further analyze the necessary conditions and the boundary for the occurrence of Turing instability. Assume that We can find that the curve ε 2 (d 1 )d 1 increases monotonically with d 1 until the first inflection point d A , after which the curve ε 2 (d 1 )d 1 decreases linearly until the second inflection point d B , and the value of ε 2 (d 1 )d 1 is less than zero when d 1 is greater than the second inflection point d B (see Fig. 2(a)). To simplify the study, we only study the part of the curve ε 2 (d 1 )d 1 before the first inflection point d A . In this part, ε 2 (d 1 )d 1 is monotonically increasing in d 1 and intersects ε 1 d 1 at the point Hence, we have the following lemma.
Since x > 0, then a 1 ε−r > 0. If we want DET n min < 0 with the condition a 1 ε−r > 0, ε must satisfy ε > (−a 1 r − 2a 2 r)/a 2 1 + 2 (a 1 a 2 r 2 + a 2 2 r 2 )/(a 4 1 ) = ε 1 . Then we can obtain When x takes a value on the axis of symmetry, n can be taken to the minimum, and we can get where d n = (a 1 l 2 + a 2 l 2 )/n 2 , then DET n = 0 when d 2 = d n 2 (d 1 ). We can find that the curve d n 2 (d 1 ) increases monotonically with d 1 until the first inflection point d n A , after which the curve d n 2 (d 1 ) decreases linearly until the second inflection point d n B , and the value of d n 2 (d 1 ) is less than zero when d 1 is greater than the second inflection point d n B (see Fig. 2(b)). For simplicity, only the part of the curve d n 2 (d 1 ) before the first inflection point d n A is studied in this section. Since both d 1 and d 2 are positive parameters, this qualification makes sense.
, has the following properties: (i) For any given In Fig. 3(a), we characterize a graph of functions , n = 1, 2, 3, which will help us understand the results of Lemmas 1 and 2. In Fig. 3(b), we present a graph of the Turing bifurcation line

Hopf bifurcation and Turing-Hopf bifurcation with r and d 2 as parameters
In this section, the existence conditions of the Hopf bifurcation with r as the parameter are first analyzed. Denote Obviously, DET 0 (r 0 ) = −a 1 (a 1 + a 2 ) > 0 under hypothesis (A2). Denote Λ = {n ∈ N 0 | DET n > 0 and r n > 0}.
After analysis, we can get the following theorem. https://www.journals.vu.lt/nonlinear-analysis Theorem 5. Suppose (A2) holds. System (3) undergoes a Hopf bifurcation at (u * , v * ) when r = r n for n ∈ Λ. Moreover, the bifurcating periodic solution is spatially homogeneous when r = r 0 and spatially nonhomogeneous when r = r n for n ∈ Λ and n = 0.
(i) When r = r n , we can get TR n (r) = 0 and DET n (r) > 0 for n ∈ Λ, then Eq. (6) has a pair of pure imaginary roots λ n (r n ) = ±i DET n (r n ).
Then we can obtain dα n (r)/dr = −1/2 < 0. That is to say, for each r n , n ∈ Λ, the transversal condition holds. This completes the proof.
Next, the existence conditions of the Turing-Hopf bifurcation with r and d 2 as parameters are analyzed. System (3) will undergo a Turing-Hopf bifurcation when the following conditions are satisfied: (i) When n = 0, Eq. (6) has a pair of pure imaginary roots ±iω. This phenomenon can be produced when r = r 0 = a 1 . (ii) When n > 0, Eq. (6) has a single zero root. This phenomenon can be produced when DET n = 0.
In this section, we assume (A2) always holds. Denote such that From the Theorem 5 we can know that system (3) will have a Hopf bifurcation at the positive equilibrium point (u * , v * ) when r = r 0 = a 1 . Therefore, r * = a 1 when n = n * . d n 2 are the Turing bifurcation lines, and r * is the Hopf bifurcation line. When n = n * , we hope to find the first intersection of these two types of bifurcation lines in the first quadrant as the Turing-Hopf bifurcation point. After the above analysis, we can get the following theorem. Q := (r, d 2 ) r > r * , 0 < d 2 < r d 1 µ n * − (a 1 + a 2 ) µ n * (a 1 − d 1 µ n * ) .
Proof. In r, d 2 -plane, the Turing bifurcation curves are The Hopf bifurcation curve is H 0 : r = r * . (i) If S = ∅, then there is no intersection point between Turing bifurcation curves L n and the Hopf bifurcation curve H 0 in the first quadrant. This indicates that system (3) does not undergo Turing-Hopf bifurcation.

Hopf bifurcation and Turing-Hopf bifurcation with K and d 2 as parameters
In this section, the existence conditions of the Hopf bifurcation with K as the parameter are first analyzed. Denote Obviously, DET 0 (K 0 ) = −r(a 1 + a 2 ) > 0 under hypothesis (A2). Denote Λ = {n ∈ N 0 | DET n > 0 and K n > 0}. Then we make the following hypotheses: After analysis, we can get the following theorem.
Next, the existence conditions of the Turing-Hopf bifurcation with K and d 2 as parameters are analyzed. In this section, we assume (A2) always holds. Denote such that . From Theorem 7 we can know that system (3) will have a Hopf bifurcation at (u * , v * ) when K = K 0 . Therefore, K = K 0 when n = n * . d n 2 are the Turing bifurcation lines, and K 0 is the Hopf bifurcation line. When n = n * , we hope to find the first intersection of these two types of bifurcation lines in the first quadrant as the Turing-Hopf bifurcation point. After the above analysis, we can get the following theorem. The proofs of Theorems 7 and 8 are similar to those of Theorems 5 and 6, respectively, and will not be repeated here.
By [7], we know: A 0 is the coexistence equilibrium; A ± 1 are spatially inhomogeneous steady states; A 2 is spatially homogeneous periodic solution; A ± 3 are spatially inhomogeneous periodic solutions.

Conclusions
In this paper, we have formulated a Leslie-Gower-type predator-prey system with the fear effect and ratio-dependent Holling III functional response to consider the cost of fear on prey demography. The main focus of this study is to investigate the influence of antipredator behaviour due to fear of predators analytically and numerically. We obtain many interesting results. Mathematically, we analyze the existence and stability of the positive equilibria of the model and give conditions for the existence of the Turing instability, the Hopf bifurcation, and the Turing-Hopf bifurcation by selecting d 2 , r, and K as the bifurcating parameters, respectively. In addition, we calculate the normal form of system (3) and divide the plane of parameters σ 1 , σ 2 into six regions D 1 -D 6 . In each region, we demonstrate that the predator-prey model exhibits complex spatiotemporal dynamics including spatially homogeneous periodic solutions, spatially inhomogeneous periodic solutions, and spatially inhomogeneous steady-state solutions. Ecologically, we show that if the model has a unique positive equilibrium, the fear effect can reduce the density of predator and prey: as the level of fear K increases, both the predator and prey density gradually decrease. Moreover, depending on the predator and prey's different states in each region, we can give solutions accordingly.