Bifurcation analysis of a predator – prey system with self-and cross-diffusion and constant harvesting rate

In this paper, we focus on a ratio dependent predator–prey system with selfand cross-diffusion and constant harvesting rate. By making use of a brief stability and bifurcation analysis, we derive the symbolic conditions for Hopf, Turing and wave bifurcations of the system in a spatial domain. Additionally, we illustrate spatial pattern formations caused by these bifurcations via numerical examples. A series of numerical examples reveal that one can observe several typical spatiotemporal patterns such as spotted, spot-stripelike mixtures due to Turing bifurcation and an oscillatory wave pattern due to the wave bifurcation. Thus the obtained results disclose that the spatially extended system with self-and cross-diffusion and constant harvesting rate plays an important role in the spatiotemporal pattern formations in the two dimensional space.


Introduction
In population dynamics, many ecologists and mathematicians are interested in Michaelis-Menten-type predator-prey model, so-called a ratio-dependent predator-prey system [2,6,10,14,20] as follows: where U and V stand for prey and predator, respectively.All parameters are positive constants.Especially, the parameters r, K, c and m stand for the prey intrinsic growth rate, the carrying capacity of prey, the capturing rate, and the half-saturation constant, respectively.The predator grows logistically with the growth rate D and the conversion efficiency f .Following [10,25,26], with the scaling

H. Baek
we can arrive at the following equations containing dimensionless quantities: From the point of view of human needs, the exploitation of biological resources and the harvest of populations are commonly practiced in fishery, forestry and wildlife management.Concerning the conservation for the long-term benefits of humanity, there is a wide range of interest in the use of bioeconomic modeling to gain insight in the scientific management of renewable resources like fisheries and forestries [24,25,26,27,30].Thus, the authors in [25,26,27] considered the following systems and studied dynamics and bifurcation phenomena of the systems: Generally speaking, in nature, the tendency of the prey would be to keep away from predators, and hence the escape velocity of the prey may be taken as proportional to the dispersive velocity of the predators.Also, the tendency of predators would be to get closer to the prey, and hence the chase velocity of predators may be considered to be proportional to the dispersive velocity of the prey.In this context, there has been considerable interest in investigating the stability behavior of systems of interacting population by taking into account the effect of self as well as cross-diffusion [5,11,17].Thus, in this paper, we will consider the following system with diffusion effects and constant harvesting rate: where d 11 , d 22 represent the positive self-diffusion coefficients and d 12 , d 21 the cross-diffusion coefficients of prey and predator, respectively and ∇ 2 = ∂ ∂x + ∂ ∂y is the usual Laplacian operator in the two-dimensional space Ω.The boundary conditions are taken as ∂u ∂n = ∂v ∂n = 0 on ∂Ω, which means that no external input is imposed from outside.Here, n is the outward unit normal vector of the boundary ∂Ω.
In [31], the authors have investigated spatiotemporal dynamics of systems (1.5) and (1.6) when Bifurcation analysis of a predator-prey system 3 The values of d 12 and d 21 imply that the prey species approaches towards the lower concentration of the predator species and the predator species tends to diffuse in the direction of higher concentration of the prey species.In many cases, the predator prefers to avoid group defence by a huge number of prey and chooses to catch its prey from a smaller concentration group unable to sufficiently resist [5,8].For this reason, it is reasonable to assume that d 12 could be any number while d 21 is positive.In addition, we assume that which indicates that self-diffusion is stronger than cross-diffusion.In other words, the flow of the respective densities in the spatial domain depends more strongly on their own density than on the others.Studies of reaction-diffusion systems have led to the characterization of three basic types of symmetry-breaking bifurcation responsible for the emergence of spatiotemporal patterns.The classification of these bifurcations is based on linear stability analysis of a homogeneous state.The space-independent Hopf bifurcation breaks the temporal symmetry of a system and gives rise to oscillations that are uniform in space and periodic in time.The (stationary) Turing bifurcation breaks spatial symmetry, leading to the formation of patterns that are stationary in time and oscillatory in space.The wave (oscillatory Turing or finite-wavelength Hopf) bifurcation breaks both spatial and temporal symmetries, generating patterns that are oscillatory in space and time [1,4,9,11,12,13,15,17,21,22].Thus, in this paper, we will focus on studying bifurcation phenomena of systems (1.5) and (1.6).

Bifurcation analysis
2.1 Bifurcations of system (1.5)In this subsection, first, we will take into account system (1.5) to investigate its dynamics and bifurcation phenomena.
In order to accomplish the purposes, we need to consider the nonspatial system (1.3) of system (1.5).
It is easy to see that simple calculation yields that there are at most four equilibria for system (1.3) as follows; where Hence for the persistence of the ecosystem, the equilibrium of the greatest interest would be an equilibrium interior to the first quadrant.In fact, if the conditions H. Baek hold, there exist exactly the four equilibria mentioned above.
From [25], the existence and the stability of the equilibria for system (1.3) is given by the following lemma.
It follows from [25] that the stability of the equilibrium point (u * 2 , v * 2 ) depends on the value of h when the condition 0 < b − d < b a < 1 hold.For this reason, throughout the paper, we assume that the conditions are satisfied.Thus from the biological point of view, we focus our concern on studying the stability behavior of the interior equilibrium point (u * 2 , v * 2 ).To simplify notation, we let (u * , v * ) stand for the equilibrium (u * 2 , v * 2 ).Diffusion is often considered as a stabilizing process, yet it is the diffusion-induced instability of a homogeneous stable steady state that results in a reaction-diffusion system's spatial pattern formation [18].It is easy to show that the equilibrium point (u * , v * ) for homogeneous system (1.3) is still a steady state for system (1.5).Now we investigate the local dynamical behavior of the spatial system (1.5) in a two-dimensional space by virtue of linear stability analysis.In order to accomplish this process, we test how perturbation of a homogeneous steady-state solution behave in the long time limit.For this we choose two-dimensional Fourier modes ū = exp((k where x = (x, y) and k = (k x , k y ) and λ is the frequency.It is worth pointing out that the terms mk x x, nk y y with m, n > 1 are not needed to be considered in the Fourier series since we take into account linear stability analysis of a spatiotemporal perturbation.Replacing u and v in system (1.5) by u * + ū and v * + v, respectively, to linearize the diffusion terms with help of Taylor expansion around the positive equilibrium point (u * , v * ).Thus we can obtain the characteristic equation where D is the diffusion matrix given by and the matrix A is given by k stands for the wave number and I represents the 2 × 2 identity matrix.
Bifurcation analysis of a predator-prey system 5 Equation (2.5) can be rewritten by an elementary calculation as (2.9) Therefore, the solutions of equation (2.8) yield the dispersion relation Now, we will discuss bifurcation phenomena of system (1.5) by considering the three bifurcation critical lines as mentioned in the Introduction.
It is well known that the onset of Hopf bifurcation corresponds to the case when a pair of imaginary eigenvalues cross the real axis from the negative to the positive side and that this bifurcation occurs only when the diffusion vanishes [16,17,21,22].Thus we can get the next theorem.
Remark 2.3.From Theorem 2.2, we can get the critical value of the Hopf bifurcation parameter h which equals to (2.12) Moreover, it follows from [25] that system (1.5) has at least one stable limit cycle if h > h H .In fact, Hopf bifurcation breaks the temporal symmetry of system (1.5), which gives rise to oscillations that are uniform in space and periodic in time with the frequency and the corresponding wavelength is where Im(z) represents the imaginary part of the number z.
In order to illustrate some numerical examples to substantiate theoretical results of this paper, we solve the partial differential equation (1.5) by discretizing the space and time and thus we regard the continuous domain Ω in system (1.5) as 200 × 200 lattice sites and we set the spacing between the lattice points to be ∆x = ∆y = 0.25.The zero-flux boundary condition is employed for the numerical simulations.We adopt a finite difference scheme for the spatial derivatives and an explicit Euler method for the time integration with the time step ∆t = 0.01.It is well known that the spatiotemporal dynamics of a diffusion-reaction system H. Baek depends on the choice of initial conditions [7,13,19,23].In the paper, an initial conditions is taken as a small amplitude random perturbation around the steady state (u * , v * ) since it is very natural from the biological point of view.We stop the simulation when the numerical solutions either reach a stationary state or show oscillatory behavior.2.1 shows snapshots of contour pictures of the time evolution of prey population in system (1.5) when d 21 = 0.004 and h = 0.118 > h H .However, it is a little bit hard to observe an oscillatory phenomenon arising from Hopf bifurcation from these snapshots.Thus, we exhibit the local phase portrait of system (1.5) for a fixed point (0.4171, 0.0338) in Figure 2.2 to observe the existence of a stable limit cycle.In addition, we can calculate numerically that the frequency of the periodic oscillations in time ω ≈ 0.1396 and the corresponding wave length λ ≈ 45.1.Also, it follows from (2.13) and (2.14) that the theoretical frequency of the periodic oscillations in time is ω H = 0.1364 and the corresponding wave length λ H = 46.0517.On the other hand, the Turing condition is one in which the uniform steady state of a reaction-diffusion partial differential equation (PDE) is stable for the ordinary differential equation, but it is unstable for the corresponding PDE with diffusion terms [1,3,8,13,31].In this context, we can get conditions that Turing bifurcation takes place in the following theorem.At the Turing threshold h T , the spatial symmetry of the system is broken and the patterns are stationary in time and oscillatory in space with the wavelength  In reaction-diffusion systems, most studies have been devoted to the study of Hopf bifurcation and Turing structures arising from Turing bifurcation [1,3,8,11,13,15,16,19,22,25,26].In the present decade, attention has turned toward patterns arising from the wave bifurcation(instability).In fact, the wave bifurcation caused by the wave instability plays an important part in pattern formations in many areas [11,12,21,28,29].Similar to Hopf bifurcation, the wave bifurcation take places when a pair of imaginary eigenvalues cross the real axis from the negative to the positive side for k = k w = 0 in equation (2.8).Thus we get the following theorem for the wave bifurcation.
Theorem 2.8.The wave bifurcation occurs if h > h w , where (2.23) Proof.For the occurrence of a wave bifurcation, the equation (2.8) must have purely imaginary roots when k is not zero.Thus the conditions α(k 2 ) = 0 and β(k 2 ) > 0 must hold for some k = 0 for the existence of the wave bifurcation.From the condition α(k 2 ) = 0, we can get Bifurcation analysis of a predator-prey system 9 Moreover, it follows from an elementary calculation that if h > h w then β(k 2 w ) > 0.  Remark 2.9.It is well known that, at the wave threshold h w , both spatial and temporal symmetries are broken and the patterns are oscillatory in space and time with the wave length where

Bifurcations of system (1.6)
As shown in [26], the expressions, which depend on all parameters in system (1.4), of the equilibria of the nonspatial system (1.4) of system (1.6) are too complicated to analyze their stabilities or bifurcation phenomena.Thus, in this subsection, we will restrict our attention to the case a = 1.Although the case a = 1 is a special case for system (1.6), mathematically and biologically the procedure of investigating for bifurcations of system (1.6) with a = 1 is generic since it can be applied to explore bifurcations of system (1.6) in other cases of parameters [26].It is from [26] that the nonspatial system (1.4) has no equilibrium points or a unique equilibrium or two equilibria according to the conditions of parameters.Among them, for the persistence of the ecosystem, we will pay attention to the case that the nonspatial system (1.4) has two equilibria, (x 1 , y 1 ) and (x 2 , y 2 ) in the first quadrant, where (2.27) Thus we focus our concern on studying the stability behavior of the interior equilibrium points (x 1 , y 1 ) and (x 2 , y 2 ).Since the necessary conditions for the existence of the two equilibria (x 1 , y 1 ) and (x 2 , y 2 ) are 0 < h < b (1+ √ b) 2 and d < b.From now on, we assume that the conditions (2.28) hold.In fact, (x 2 , y 2 ) is a hyperbolic saddle and (x 1 , y 1 ) is a focus (or a center or a node) [26].
For this reason, in this section, we consider the positive equilibrium point (x 1 , y 1 ).We denote it briefly by (x * , y * ), where Proof.See [26].
Theorem 2.12.Turing bifurcation occurs if h > h T .Here, the critical value h T satisfies the equation

Conclusion and discussion
In this study, we investigate bifurcation phenomena of a ratio dependent predator-prey system with self-and cross-diffusion and constant harvesting rate.It has been shown from the bifurcation analysis that the system has Turing, Hopf and wave bifurcations.Furthermore, we demonstrate by numerical simulations that typical Turing patterns such as spotted and spotstripelike mixtures patterns, a stable limit cycle due to Hopf bifurcation and an oscillatory wave phenomenon can be emerged as the harvesting rate h varies.In fact, in many researches [1,3,13,15,16,19,21,22,25,26], investigation of the spatial pattern of a predator-prey system with self-diffusion have been done.It revealed that regular Turing patterns can be induced by self-diffusion.Moreover, the authors in [5,8,17] have shown that both self-and cross-diffusion can lead the population to be in regular distribution.

H. Baek
Especially, they looked into that Turing patterns can be induced by different mechanisms, not only self-diffusion.In the present paper, the great difference between the mathematical systems in [5,8,17] and the system dealt with in the paper is consideration of the prey or predator harvesting.Our bifurcation analysis and numerical computations unveil that Turing patterns including spotted and spot-stripelike mixtures patterns can also be observed, which are similar phenomena to a predator-prey system without any harvesting rate, though the harvesting effect exists.Based on these facts, we reckon that harvesting plays a significant role in pattern formation in predator-prey systems with self-and cross-diffusion.In addition, we studied the wave bifurcation, which is also called oscillatory Turing or finite-wavelength Hopf bifurcation.More recently, attention has turned toward patterns arising from the wave bifurcation [12,28].From the wave bifurcation, oscillatory wave patterns which are different from spatial patterns caused by Turing bifurcation can be obtained.Thus, this research will be useful for understanding the dynamic complexity of ecosystems or physical system when there is harvesting effect on prey or predator population.
and 2.4, typical Turing patterns, called stripelike and spot-stripelike patterns, of prey in system (1.5) can be observed for h = 0.1174 and h = 0.118, respectively.

Example 2 . 10 .
Now, consider the same parameter values as in Example 2.7 except for the values d 21 = 0.001 and h = 0.117.It follows from Theorem 2.8 that the wave bifurcation occurs under h T = 0.1173 > h > h w = 0.0126.In fact, if one takes the harvesting rate h = 0.117, one can observe wave phenomena numerically as shown in Figure 2.5.Figure 2.5 (d) is a space-time plot obtained by piling up the prey population for the lines y = 100 in each snapshots as time progresses.Also, numerical calculations yield that the frequency of the periodic oscillations in time k w ≈ 0.3740 and the corresponding wavelength λ w ≈ 45.1.From (2.24) and (2.25), the values of the frequency and the corresponding wavelength can be obtained as k w = 0.4226 and λ w = 35.1791,respectively.