Turing–Hopf Bifurcation Analysis in a Diffusive Ratio-Dependent Predator–Prey Model with Allee Effect and Predator Harvesting

This paper investigates the complex dynamics of a ratio-dependent predator–prey model incorporating the Allee effect in prey and predator harvesting. To explore the joint effect of the harvesting effort and diffusion on the dynamics of the system, we perform the following analyses: (a) The stability of non-negative constant steady states; (b) The sufficient conditions for the occurrence of a Hopf bifurcation, Turing bifurcation, and Turing–Hopf bifurcation; (c) The derivation of the normal form near the Turing–Hopf singularity. Moreover, we provide numerical simulations to illustrate the theoretical results. The results demonstrate that the small change in harvesting effort and the ratio of the diffusion coefficients will destabilize the constant steady states and lead to the complex spatiotemporal behaviors, including homogeneous and inhomogeneous periodic solutions and nonconstant steady states. Moreover, the numerical simulations coincide with our theoretical results.


Introduction
The predator-prey interaction is a significant topic in the studies of populations, communities, and ecosystems and has attracted much attention from scholars.Since the introduction of the classical Lotka-Volterra model, predator-prey models have been continuously improved and developed, but there are still many ecological problems that need attention and solving [1][2][3][4].We can use the following form of an ordinary differential equation to represent the predator-prey system where u and v are the densities of prey and predators, respectively.Here, f 1 (u) is the natural growth rate of the prey population in the absence of predators, f 2 (v) is the predator growth rate without prey, c is the conversion rate from predation, and p(u, v) is the functional response.Generally, functional responses can be divided into two categories: prey-dependent functional response and predator-dependent functional response.Preydependent functional response means that p(u, v) is the function of prey density only; see [4][5][6][7], for example.However, in predator-dependent predator-prey models, the function p(u, v) depends on the densities of both predator and prey.In [8], the authors proposed that a functional response should depend on the ratio of prey-to-predator density, which is supported by several laboratory observations and has been widely used in predator-prey systems.Since then, the ratio-dependent predator-prey models have gained much attention over many decades.Compared with the traditional prey-dependent predator-prey models, the ratio-dependent models can exhibit richer and more reasonable dynamical behaviors [9][10][11][12][13].
In fact, there are still other factors, as well as the functional response, that should be considered in the study of the predator-prey systems.Stephens et al. [14] suggested that the Allee effect, proposed by Allee [15], can be defined as a positive relationship between any component of individual fitness and the number or density of conspecifics.In the 2000s, Kramer et al. [16] suggested that the Allee effect can be caused by mate limitation, predator satiation, cooperative feeding or defense, habitat alternation, dispersal, etc.Their study showed that the Allee effect plays a key role in numerous systems.Later, Merdan [17] studied the effect of the Allee effect on the stability of a Lotka-Voterra model.The study demonstrated that the presence of the Allee effect makes the system take a longer time to reach the stable equilibrium and reduces the densities of both prey and predators at the stable equilibrium.Due to the significant effect on population dynamics, the Allee effect has gained increasing attention; see [18][19][20][21][22][23][24][25][26][27], for example.In view of human needs, the harvesting of biological resources should also be taken into account in predator-prey models.In [28], Xiao et al. considered the ratio-dependent predator-prey model with constant harvesting of predators as follows: where h represents the harvesting rate; for more background about (1), we refer readers to [28] and the references therein.For system (1), the authors obtained the occurrence of numerous types of bifurcations, including the Bogdanov-Takens bifurcation, the saddle-node bifurcation, and the supercritical and subcritical Hopf bifurcations.From the perspective of biology, it is more reasonable that the harvesting rate function should be proportional to the harvested population.In [29], Chakraborty et al. proposed the modified ratio-dependent predator-prey system with nonconstant predator harvesting as follows: where h > 0 is the harvesting rate; for more details concerning system (2), we refer readers to [29] and the references therein.This study revealed that when the harvesting rate is very high or low, the predator will eventually be extinct.In [30], the authors discussed the local stability of system (2).They gave the conditions under which the interior equilibrium is stable or unstable, a focus or a center.Their study showed that the predator harvesting rate plays a key role in the stability of the interior equilibrium of system (2), and the presence of predator harvesting makes system (2) exhibit much richer dynamical behaviors.
As is well known, the predators and prey distribute inhomogeneously in different locations; therefore, diffusion should be taken into account in ecological and biological models.Based on the fact that diffusion may destabilize the steady state and induce the occurrence of Turing instability, many scholars have investigated the diffusive predatorprey systems; see [13,[31][32][33][34], for example.Although the ratio-dependent predator-prey system has been extensively investigated, a study concerning the system incorporating the Allee effect, predator harvesting, and diffusion has not been seen yet.Based on this, we consider a diffusive system with the Allee effect in prey and predator harvesting as follows: where K is the carrying capacity for the prey, γ is the maximum per capita growth rate of the prey, G is the capturing rate, µ is the conversion rate, L is the predator death rate, and H is the harvesting rate.Further, Q 0 (Q 0 < K) is the Allee threshold, and d 1 and d 2 are the diffusion coefficients.Let and, by dropping the hats of the notations, we can obtain the corresponding diffusive system with a homogeneous Neumann boundary and initial conditions as follows: For system (4) without the Allee effect, Gao et al. [34] studied the existence and properties of a Hopf bifurcation, provided the conditions for the occurrence of Turing instability induced by diffusion, and proved the existence and non-existence of the non-constant steady states.When the harvesting term in ( 4) is absent, i.e., h = 0, Rao and Kang [13] investigated the effect of diffusion and the Allee effect on the dynamical complexity of the system.Their results reveal that the strength of the Allee effects plays a key role in the formation of distinct spatial patterns.
In this paper, we aim to explore the joint effect of diffusion and harvesting effort on the dynamics in system (4).Notice the fact that the term uv u+v has no definition at (0, 0); we assume that uv u+v (0,0) = 0 as in [9,10].The rest of this article is organized as follows.In Section 2, we first discuss the existence of positive equilibria.Then, we investigate the dynamics of the ODE system corresponding to system (4).In particular, using the harvesting rate as the bifurcation parameter, we study the stability of the positive equilibria, verify the existence of a Hopf bifurcation, and derive the explicit formulas for determining the properties of the bifurcating periodic solutions by applying the center manifold theory and normal form method.In Section 3, we give the sufficient conditions for the occurrence of the Turing-Hopf bifurcation.In Section 4, to illustrate the complex dynamics of system (4), we calculate the normal form near the Turing-Hopf bifurcation point.In Section 5, we give some numerical simulations to illustrate our theoretical results.

Dynamics of the ODE Model
When diffusion is absent, system (4) becomes It is easy to see that system (5) always has three boundary equilibria: E 0 = (0, 0), E 1 = (b, 0), and E 2 = (1, 0).Obviously, the interior equilibria should satisfy the following equation which yields to Therefore, when the following condition , and Then, we analyze the stability of all equilibria of system (5).Note that the Jacobian matrix cannot be evaluated at E 0 since uv u+v is not differentiable at (0, 0).We can obtain, from the first equation of (5), that Similarly, from the second equation of system (5), we have that which means u(t) > 0 and v(t) > 0 when u(0) > 0, v(0) > 0. From the equation for the prey density, we have for 0 It can be proven that u(t) < u 0 e −aAt , thus lim so we mainly consider two cases: Case I: If v(t) has infinite extremum for t > 0, then the maximal values are determined by dv/dt = 0, which can be given by for a sufficiently large t.It is easy to prove that v(t) < v 0 e −(1−ε)(d+h)t , which implies that v(t) → 0 when t → +∞.It contradicts the assumption B > 0. Consequently, lim t→+∞ v(t) = 0. Based on the above discussion, we obtain that E 0 is locally asymptotically stable.For E i (i = 1, 2), the corresponding Jacobian matrix can be given by respectively, so we can obtain the local stability of boundary equilibria as follows.
Lemma 1.For system (5): ) can be written as and TrJ DetJ According to Equations ( 7) and (10), we can conclude that TrJ When (H1) holds, we can obtain that lim So, under the assumptions (H1) and Summarizing the previous discussion, we conclude that: Theorem 1. Assume that (H1) is satisfied.Then, system (5) has two positive equilibria where h H is the positive zero of (11).
In fact, we can also obtain the conclusion as follows.

Proof.
Denote where Let λ(h) = κ(h) ± iω(h) be a pair of complex roots of Since For the sake of convenience, we still denote ũ and ṽ by u and v. Thus, system (5) can be rewritten as Rewrite system (14) as where with Let Then, we obtain that Using (u, v) T = P(p, q) T , system (15) becomes where The polar coordinate form of ( 16) can be written as: By Taylor expanding (17 To determine the stability of the bifurcating periodic solutions, we need to discuss the sign of α(h H ), which is given by where all partial derivatives are evaluated at the bifurcation point (p, q, b) = (0, 0, h H ), and Thus, we can obtain the sign of the coefficient α(h H ) in (19).Note that κ ′ (h H ) < 0; we draw the following conclusions.
Theorem 2. Suppose that (H1) and (H2) hold.Then, system (5) undergoes a Hopf bifurcation at E * 2 when h = h H . Furthermore: (i) If α(h H ) < 0, the bifurcating periodic solutions are orbitally asymptotically stable, and periodic solutions occur as h decreases and passes h H . (ii) If α(h H ) > 0, the bifurcating periodic solutions are unstable, and periodic solutions occur as h decreases and passes h H .

Turing Instability Induced by Diffusion and Turing-Hopf Bifurcation
In this section, we discuss the effect of diffusion on the stability of E * 2 and give the sufficient conditions for the occurrence of Turing instability induced by diffusion.By computation, we have that the characteristic equations corresponding to E * 2 can be given by which can be equivalent to where with a ij (i, j = 1, 2) defined as in (13).In the following analysis, we always assume that (H1) and (H2) are satisfied.
(ii) The Turing-Hopf bifurcation occurs at E * 2 when (h where Proof.From Theorem 1, we know that Turing instability occurs as ∆ n (λ) has roots with a positive real part for some n ∈ N when h > h H . From (20), T n < 0 can be satisfied provided h > h H . Thus, we only need to seek the condition for J n = 0 to ensure the occurrence of Turing instability.In fact, J n = 0 is equivalent to .
By calculation, we have where u * 2 (h) is a function of h defined as in (7).Since Thus, lim > 0 holds true.Denote To guarantee the positiveness of D 0 on the curve D 0 = S n (h), we should ensure that the condition h < h (n) holds.Denote the Turing bifurcation curve as L n , that is When L n intersects the critical Hopf bifurcation curve h = h H at (h H , D * 0(n) ), system (3) undergoes the Turing-Hopf bifurcation at E * 2 as (h, D 0 ) = (h H , D * 0(n) ).
For n ≥ k * , to seek the maximum of D * 0(n) , we take the derivative of D * 0(n) with respect to n 2 .We can have that In fact, dD * 0(n) dn 2 has the same sign as . We can see that So, there must exist a x * > 0 satisfying Φ(x * ) = 0 and Φ(x we can obtain that D * 0(k * 0 ) = max n∈N D * 0(n) .Then, we complete the proof.
By calculation, we obtain that , , Following the techniques in [35], by a recursive transformation, we can obtain that the normal form for the Turing-Hopf bifurcation can be given by where and Next, we need to calculate C ijk , D ijk , and E ijk .
Then, we can figure out by calculation that Thus, we can obtain So, the normal form truncated to the third-order terms for the Turing-Hopf bifurcation can be written as where , κ 22 = sign(B 003 ).
Then, we obtain the bifurcation diagram in the µ 1 − µ 2 plane and the corresponding phase portraits of system (28) in the ρ − ς plane, as shown in Figure 3. Clearly, the above curves divide the µ 1 − µ 2 plane into six regions, denoted as R i , i = 1, 2 • • • , 6 (see Figure 3).The existence and stability properties of the steady states in the six regions are listed in Table 1.Obviously, the equilibria Ã0 , Ã1 , Ã± 2 , and Ã± 3 of system (28) correspond to the constant equilibrium, the spatially homogeneous periodic solution, the nonconstant steady state, and the spatially inhomogeneous periodic solution of system (4).Thus, the dynamical behaviors of system (4) nearby the Turing-Hopf singularity in the h − D 0 plane can be determined by the dynamical behaviors of system (28).

Conclusions and Discussion
In this paper, we investigate the ratio-dependent predator-prey system with the Allee effect in prey and predator harvesting.We give a detailed analysis of the joint effect of harvesting effort and diffusion on the spatiotemporal behaviors of system (4), and our results reveal that the presence of a harvesting term makes the system exhibit more interesting dynamical behaviors.
A ratio-dependent predator-prey system with a harvesting term is a relatively new issue that has been investigated by several researchers and has yielded many interesting results.Recently, Gao et al. [34] analyzed the existence of a Hopf bifurcation induced by harvesting rate and a Turing bifurcation induced by diffusion, respectively, for systems (4) without the Allee effect.However, they did not consider the dynamics of multi-parameter synergism, which is also a major difficulty in our research.
The main contribution of our paper is a detailed analysis of bifurcation near the posi-

Conclusions and Discussion
In this paper, we investigate the ratio-dependent predator-prey system with the Allee effect in prey and predator harvesting.We give a detailed analysis of the joint effect of harvesting effort and diffusion on the spatiotemporal behaviors of system (4), and our results reveal that the presence of a harvesting term makes the system exhibit more interesting dynamical behaviors.
A ratio-dependent predator-prey system with a harvesting term is a relatively new issue that has been investigated by several researchers and has yielded many interesting results.Recently, Gao et al. [34] analyzed the existence of a Hopf bifurcation induced by harvesting rate and a Turing bifurcation induced by diffusion, respectively, for systems (4) without the Allee effect.However, they did not consider the dynamics of multi-parameter synergism, which is also a major difficulty in our research.
The main contribution of our paper is a detailed analysis of bifurcation near the positive constant steady state of system (4) in the one-dimensional spatial domain (0, lπ).For the spatially homogeneous model, using the harvesting rate h as the bifurcation parameter, we analyze the stability of interior equilibria.By applying the center manifold theory and normal form method, we derive the formula determining the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions.For the reaction-diffusion model, we firstly verify the existence of Turing instability induced by diffusion, which reveals the existence of spatially inhomogeneous patterns, including the spatially inhomogeneous periodic solutions and non-constant steady-state solutions.Then, the normal form near the Turing-Hopf bifurcation point is derived.Our study demonstrates that as the harvesting rate h decreases and passes the critical value h H , the coexistence equilibrium E * 2 will lose its stability, and the Hopf bifurcation occurs.Finally, numerical simulations, which are consistent with the theoretical results, are performed to illustrate the theoretical analysis.
In ecosystems, humans, as higher animals, have the ability to harvest biological resources.Our research shows that if humans overharvest biological resources, it will lead to the unbalance of the ecosystem.Our study provides the critical harvesting rate h H without destroying the ecosystem, which not only ensures the health of the ecosystem but also maximizes the biological resources available to humans.At the same time, our results also show that humans can obtain better harvest times and easier harvest locations by controlling the ratio of the diffusion coefficients of prey and predators and the harvesting rate near the Turing-Hopf singularity (h H , D * 0(n) ).This can greatly reduce the difficulty for humans in obtaining biological resources.
In [36], the authors proposed the non-continuous harvesting function as follows H(x) = 0, x < T, qx, x > T.
They assumed that the harvesting may start as the population reaches some critical value T. To further the study, we will consider a non-continuous harvesting function in our system, compare the theoretical results to the results in the paper, and reveal the effect of the harvesting rate on the predator population.Moreover, we may consider the effect of spatial heterogeneity on the dynamics of (4).
then E 1 is a source, and E 2 is a saddle; (iii) If c < d + h, then E 1 is a saddle, and E 2 is a stable node.

Table 1 .
(28)ility of the steady states in different regions of system(28).