# Spatiotemporal Dynamics of a Predator-Prey System with Linear Harvesting Rate

Authors: Hunki Baek
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We investigate a spatial ratio-dependent predator-prey system with linear harvesting rate. By
using linear stability and bifurcation analysis, we obtain the conditions for Hopf and Turing bifurcation in the
spatial domain. In addition, we classify spatial pattern formations of the system by making use of numerical
simulations. In fact, the numerical simulations unveil that the typical Turing patterns, such as spotted, spot-stripelike
mixtures and stripelike patterns, can be observed even if the system has the linear harvesting rate. In
order to analyze these patterns via the spatial frequency, the discrete Fourier transform is used. Moreover, we
discuss that the linear harvesting system is more realistic than a predator-prey system with constant harvesting.
Our results disclose that the spatially extended system with linear harvesting rate has more complex dynamic
patterns in the two-dimensional space. It may help to understand the effects of harvesting on species in the real
world.


---

## Body

## 1. Introduction

Since a classical Holling-Tanner type predator-prey system, so-called a ratio-dependent predator-prey system, has received great attention among mathematicians and mathematical biologists (see [1–4]), we will focus our concern on the following system of nonlinear coupled ordinary differential equations:
(1)dNdt=rN(1-NK)-cPNN+mP,dPdt=-δP+sPNN+mP,
where N and P 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 logically with the growth rate δ and conversion efficiency s. We adopt a ratio-dependent functional response in system (1) which roughly the per capita predator growth rate should be a function of the ratio of prey to predator abundant. One of the main reasons to consider such response function is that numerous laboratory experiments and observations showed that functional and numerical responses over typical ecological timescales ought to depend on the densities of both prey and predator, especially when predators must search for food (see [1, 5–8]).Following Kuang and Beretta in [3] and Xiao et al. in [7, 8], with the next scaling
(2)NK=N¯,mPK=P¯,t=Tr,a=cmr,b=sr,d=δr,
we can arrive at the following equations containing dimensionless quantities by omitting bar notations;
(3)dNdT=N(1-N)-aPNN+P,dPdT=-dP+bPNN+P.Although such kinds of predator-prey systems play an important role in studying about dynamics of prey or predator populations, we cannot infer any information about spatial distributions of prey or predator from these ordinary differential systems. Thus, it is needed to investigate spatial systems. In this context, recently, many authors in [9–15] have investigated spatiotemporal pattern formations and bifurcation analysis of predator-prey systems with reaction diffusion. Thus, firstly, we will take into account the following a ratio dependent predator-prey system with reaction diffusion:
(4)dNdT=D1∇2N+N(1-N)-aPNN+P,dPdT=D2∇2P-dP+bPNN+P,
where D1,D2 are the diffusion coefficients of prey and predator, respectively, and ∇2=∂/∂x+∂/∂y is the usual Laplacian operator in two-dimensional space Ω.In fact, Banerjee [10] and Wang et al. in [13] have studied the spatiotemporal complexity of a ratio-dependent predator-prey system similar to system (4).From the point of view of human needs, the exploitation of biological resources and the harvest of population 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 bioeconomics modeling to gain insight in the scientific management of renewable resources like fisheries and forestries (see [7, 8]). Thus, it is necessary to investigate the effect of the harvesting of species in some models. In fact, many researchers have studied complex bifurcation phenomena and dynamics of the models with nonzero constant rate of harvesting of either species or both species simultaneously (see [7, 8, 16–18]). However, in general, human needs are not immutable for a long time. Actually, human needs could be changed according to the amount of biological resources. For instant, let us think about an outbreak of algal bloom caused by phytoplankton or zooplankton (see [12, 19]). In this case, harvesting the populations of phytoplankton or zooplankton is one of simple ways to control the algal bloom. It is an undoubted fact that a nonconstant harvesting method considering prey or predator populations is more efficient and economic than the constant harvesting method. For the reason, it is needed to take into account nonconstant harvesting rate. Such harvesting activity has an effect on both species, prey and predator, directly. Thus, with the ideas discussed above, in the paper, we consider the following ratio-dependent predator-prey system with reaction-diffusion and linear harvesting rates h and u:
(5)dNdT=D1∇2N+N(1-N)-aPNN+P-hN,dPdT=D2∇2P-dP+bPNN+P-uP,
where h,u are nonnegative constants.However, since-dP-uP=-(d+u)P, the dynamics of system (5) are the same as the following system:
(6)dNdT=D1∇2N+N(1-N)-aPNN+P-hN,dPdT=D2∇2P-dP+bPNN+P,
where the boundary conditions are taken as ∂N/∂n=∂P/∂n=0 on ∂Ω, which means that no external input is imposed from outside. Here, n is the outward unit normal vector of the boundary ∂Ω.Thus, from now on, we will focus on studying spatial dynamics of system (6) instead of system (5).The main purpose of this paper is to investigate spatial dynamics of system (6). In Section 2, we investigate bifurcation phenomena and Hopf and Turing bifurcations by making use of linear stability theory and bifurcation analysis. In addition, we obtain the conditions of Hopf and Turing bifurcation in a spatial domain. In Section 3, using numerical simulations for system (6), we show that the spatial system (6) has various spatiotemporal patterns such as spotted and stripelike patterns. Also, we obtain the spatial frequency via discrete Fourier transforms of the spotted and stripelike patterns to analyze these patterns in terms of periodic components. Finally, we give conclusions and discussions in Section 4. Especially, we discuss that the linear harvesting rate is more realistic than the constant harvesting by comparing the above linear harvesting system with the constant harvesting system studied by Zhang et al. in [15].


## 2. Bifurcation Analysis

In the beginning, we will think about the spatially homogeneous system of the diffusive system (6) to look into spatial pattern formations of system (6). It is fundamental to mention that the nonspatial system of system (6) has at most three nonnegative equilibria, which correspond to spatially homogeneous equilibria of system (6), in the first quadrant ℝ+2,
(7)E0=(0,0),E1=(1-h,0),E*=(N*,P*),
where N*=1-a-h+ad/b and P*=((b-d)/d)N*. Moreover, the point E* lies in the interior of the first quadrant if 0<h<1-a+ad/b and b>d.It follows from [20] that the point E0 is a saddle point and the point E1 is either a saddle point when b-d>0 or a stable node when b-d≤0 for the nonspatial system of system (6).From the biological point of view, it is important to mainly take into account the dynamics of system (6) at the nontrivial stationary state E*=(N*,P*), which corresponds to the coexistence of prey and predator. Throughout the paper, we assume that parameters in system (6) satisfy the conditions
(8)0<h<1-a+adb,b>d.
Under condition (8), the equilibria E0 and E1 are both saddle points.In order to deal with the stability and bifurcation analysis of the spatiotemporal system (6), we linearize the dynamical system (6) around the spatially homogeneous fixed point E* for small space and time fluctuations. To accomplish this purpose, we assume that
(9)N(x→,t)~N*eλteik→·x→,P(x→,t)~P*eλteik→·x→,
where x→=(x,y) and k→·k→=k2, and k and λ are the wavenumber vector and frequency, respectively. Then we can obtain the corresponding characteristic equation |Ak-λI|=0, where Ak=A-k2D and D=diag⁡(D1,D2) is the diffusion matrix and the matrix A is given by
(10)A=(∂F∂N∂F∂P∂G∂N∂G∂P)(N*,P*),
whereF=N(1-N)-aPN/(N+P)-hN,G=-dP+bPN/(N+P). From an elementary calculation, we can get
(11)A=(1-2N*-aP*2(N*+P*)2-h-aN*2(N*+P*)2bP*2(N*+P*)2-d+bN*2(N*+P*)2)≡(FNFPGNGP).
Then the solutions of the characteristic equation yield the dispersion relation
(12)λk±=12(tr⁡(Ak)±tr⁡(Ak)2-4det⁡(Ak)),
where
(13)tr⁡(Ak)=FN+GP-k2(D1+D2)hhhhlh=-(D1+D2)k2-1+a-d+h+(b-a)(db)2,det⁡⁡(Ak)=FNGP-FPGNtr⁡(Ak)kklk-k2(FND2+GPD1)+k4D1D2tr⁡(Ak)h=D1D2k4+(ad2D2-bd2D1b2dD1+(1-a-h)D2hhhhhhhhhhhhhhhhhhh+ad2D2-bd2D1b2)k2hhhhlllhhh+d(b-d)(b-ab+ad-bh)b2.It is well known that reaction-diffusion systems have led to the characterization of two basic types of symmetry-breaking bifurcations, Hopf and Turing bifurcations, which are responsible for the emergence of spatiotemporal patterns (see [9–15, 21–27]).The space-independent Hopf bifurcation breaks temporal symmetry of a system and gives rise to oscillations that are uniform in space and periodic in time (see [9, 13, 15, 20]). This bifurcation takes place only when the diffusion vanishes. Thus, mathematically speaking, the Hopf bifurcation occurs if the conditions
(14)Im⁡(λk)≠0,Re(λk)=0atk=0
hold. From simple calculation, it is easy to know that the critical value of Hopf bifurcation parameter h is equal to
(15)h=1-a+d+(a-b)(db)2≡hH.
From now on, for the existence of the critical value of Hopf bifurcation, it is assumed that the parameters are satisfied with the following inequality:
(16)(1-a+d)b2+(a-b)d2>0.
Thus the nontrivial stationary state E*=(N*,P*) of system (6) is an unstable equilibrium if h>hH. Furthermore, as mentioned above, the bifurcation breaks the temporal symmetry of the system, which causes oscillations that give rise to uniform in space and periodic in time with the frequency
(17)ωH=Im⁡(λ0)=det⁡⁡(A)=d(b-d)(b-ab+ad-bh)b2
and the corresponding wavelength is
(18)λH=2πωH=2πb2d(b-d)(b-ab+ad-bh).On the other hand, Turing bifurcation is sometimes called Turing instability or diffusion-driven instability. It is well known that Turing bifurcation breaks spatial symmetry, leading to the formation of patterns that are stationary in time and oscillatory in space. Moreover, since this instability occurs only if the activator (N) diffuses more slowly than inhibitor (P), the condition D1≪D2 is assumed.In fact, Turing bifurcation occurs when the positive steady state(N*,P*) of system (6) is stable for the nonspatial ordinary differential system of system (6) and it is unstable for the diffusive system (6) (see [9, 10, 12, 13, 15, 21–23, 28]). The conditions for the steady state (N*,P*) to be stable are given by
(19)tr⁡(A0)=FN+GP<0,det⁡(A0)=FNGP-FPGN>0.
Since tr⁡(A0)<0, tr⁡(Ak) is always negative and hence the two roots λk± of the characteristic equation cannot be positive at the same time, therefore Turing  instability can happen only if det⁡(Ak)<0 because the stationary state is unstable to spatial perturbations for k≠0. To hold the condition det⁡(Ak)<0, the coefficient of k2 must be negative because of det⁡(A0)>0; that is,
(20)dD1+(1-a-h)D2+ad2D2-bd2D1b2<0.
However, condition (20) is necessary, not a sufficient condition. Thus we must require that the minimum value of det⁡(Ak) is below zero. Actually, det⁡(Ak) has the minimum value at
(21)kT2=12D1D2(bd2D1-ad2D2b2-dD1+(-1+a+h)D2hhhhhhhhll+bd2D1-ad2D2b2),
which is called the wavenumber. Therefore the value of det⁡(Ak) at k2=kT2 must be negative; that is,
(22)D1D2kT4+(dD1+(1-a-h)D2+ad2D2-bd2D1b2)kT2+d(b-d)(b-ab+ad-bh)b2<0.
In order to get the threshold of Turing bifurcation, let the inequality in (22) be the equality. Then we have the critical value of the bifurcation parameter h for the Turing instability as follows:
(23)h=1-a+ad2b2+(db-1)dD1D2±2d(b-d)abD1D2b2D2≡hT.At the Turing thresholdhT, the spatial symmetry of the system is broken and the patterns are stationary in time and oscillatory in space with the wavelength
(24)λT=2πkT.In order to discuss the bifurcations represented by the above formulae in the parameter space spanned by the parametersD1 and h, we will observe a bifurcation diagram of system (6) by letting a=1.3, b=0.8, d=0.4, and D2=0.2 in system (6). The bifurcation diagram shown in Figure 1 can be obtained using the linear stability analysis studied above. In this case, we consider the range D1≤0.2 to guarantee the existence of Turing instability. As seen in Figure 1, the bifurcation parametric space is separated into four distinct domains by the Hopf and the Turing bifurcation lines. In domain I located below two bifurcation lines, the steady state is the only stable solution of system (6). On the other hand, domains II and III are the regions of the pure Turing instability and the pure Hopf instability, respectively. Hopf and Turing instabilities take place in domain IV simultaneously.Figure 1
Bifurcation diagrams of system (6) with respect to D1 and h when a=1.3, b=0.8, d=0.4, and D2=0.2. Hopf bifurcation line: h=0.2250; Turing bifurcation line: h=(1/2)(26/5)D1-D1+1/40. Turing-Hopf bifurcation point: (0.0469,   0.2250).


## 3. Pattern Formation Analysis

In this section, we will investigate spatiotemporal pattern formations of the spatially extended system (6) via computational simulation in two-dimensional space. For this, let a=1.3, b=0.8, d=0.4, and D2=0.2. In order to solve the system (6) numerically, we employ a finite difference scheme for the spatial derivatives, an explicit Euler method for the time integration with the time step Δt=0.01, and the discrete domain of 200×200 lattice sites with the space step Δx=Δy=0.25. These numerical schemes for the diffusion equation (6) give stable solutions as long as the equalities D1Δt/(Δx)2,D2Δt/(Δy)2<1/2 are satisfied (see [29]).It is well known that spatiotemporal dynamics of a diffusion-reaction system depend on the choice of initial conditions (see [12, 13, 28]). However, from the biological point of view, it is natural to take an initial condition as a small amplitude random perturbation around the steady state (N*,P*). We run the simulations until the numerical solutions either reach a stationary state or show a behavior that does not seem to change its characteristic anymore.In the numerical simulation, we restrict our analysis of pattern formations to the distribution of prey since we have figured out that the distributions of predator and prey are always of the same type. Nevertheless, through the numerical simulation different types of spatiotemporal dynamics of the distribution of prey are observed.Now, fixing the parameter valueD1=0.01, the bifurcation thresholds, hH≈0.225 for the Hopf bifurcation, and hT≈0.129 for the Turing bifurcation can be obtained using (15) and (23), respectively (see Figure 1). In addition, due to Turing instability, three typical Turing patterns such as spotted, spot-stripe mixture, and stripelike patterns of prey (N) in system (6) can be observed accordingly as the harvesting rate h varies increasingly. Snapshots that represent such phenomena are illustrated in Figures 2, 3, and 4 when taking the harvesting rate h=0.13,0.18 and 0.22, respectively.Snapshots of contour pictures of the time evolution of prey in system (6) when D1=0.01 and hT<h=0.13<hH: (a) 0 iteration; (b) 400000 iterations; (c) 1000000 iterations.
(a)
(b)
(c)Snapshots of contour pictures of the time evolution of prey in system (6) when D1=0.01 and hT<h=0.18<hH: (a) 0 iteration; (b) 10000 iterations; (c) 1000000 iterations.
(a)
(b)
(c)Snapshots of contour pictures of the time evolution of prey in system (6) when D1=0.01 and hT<h=0.22<hH: (a) 0 iteration; (b) 10000 iterations; (c) 390000 iterations.
(a)
(b)
(c)To analyze these patterns by using the spatial frequency, we make use of the discrete Fourier transform. It is known that the discrete Fourier transform is one of useful mathematical tools to analyze a signal or an image by decomposing it into different periodic components. Especially, it has been widely used for spatial patterns (see [13, 26, 30]). Thus, we apply the spatial Fourier transform to the Turing patterns in Figures 2, 3, and 4 to get the spatial frequency as in Figure 5. The spatial frequency and direction of any component in the power spectrum are given in the length and direction, respectively, of a vector from the origin to the point on the circle (see [13, 26, 30]). We can infer from Figure 5 that Figures 2 and 3 have similar spatial frequency in the length of the space unit and presence of one mode with different wavelengths while Figure 4 has two modes with different wavelengths. Since the wavelength is defined by the inverse of the frequency, the small and large circles in Figure 5(c) have long and short wavelength, respectively.Figure 5
Spatial Fourier transformation (left-hand column) and radial average of the power spectrum (right-hand column) corresponding to the snapshots (c) of Figures2, 3, and 4 when D1=0.01. (a, d) h=0.13; (b, e) h=0.18; (c, f) h=0.22.If we selectD1=0.16 and h=0.23 in Domain III, which is the region of the pure Hopf instability, we can investigate influence of the Hopf instability phenomena on system (6). Of course, it is not easy to see the influence of the instability directly through snapshots as shown in Figure 6. However, an oscillatory pattern in time can be observed. From ecological point of view, this bifurcation gives rise to the coexistence of two species, prey and predator, in the long run. In order to investigate this phenomenon more clearly, it is needed to see a local phase portrait of system (6). For this, let us fix a spatial point (0.1196,0.1195); then the local phase portrait of system (6) for this point can be obtained, exhibited in Figure 7. Therefore, we can assert from Figure 7 that the stable limit cycle occurs by virtue of the Hopf instability. Also, we can estimate numerically the frequency of the periodic oscillations in time ωH≈40.6 and the corresponding wavelength λH≈0.1548. In fact, it is shown from (17) and (18) that the frequency of the periodic oscillations in time is ωH≈40.5578 and the corresponding wavelength is λH≈0.1549.Snapshots of contour pictures of the time evolution of prey in system (6) when D1=0.16 and hH<h=0.23<hT: (a) 0 iteration; (b) 50000 iterations; (c) 400000 iterations.
(a)
(b)
(c)Figure 7
Dynamical behavior of system (6) for a point (0.1892,0.2270) when D1=0.16 and h=0.23: (a) time series of prey; (b) time series of predator; (c) the stable limit cycle of system (6).


## 4. Conclusions and Discussions

In this study, we have investigated pattern formations of a ratio-dependent predator-prey system with linear harvesting rate. It has been shown from the bifurcation analysis that the system has Turing and Hopf bifurcations. For some fixed parameters, the numerical simulations illustrate that the spatial domain is divided into four regions by two bifurcation lines, Hopf and Turing bifurcation lines, as shown in Figure1. In domain II of Figure 1, as the linear harvesting rate h increases typical Turing patterns such as spotted, spot-stripelike mixtures, and stripelike patterns are emerging and in domain III, a stable limit cycle emerges due to the Hopf bifurcation. In addition, applying the spatial Fourier transform to such Turing patterns, we classify these patterns in terms of the spatial frequency.As mentioned above, we think about spatiotemporal dynamics of a predator-prey system with linear harvesting rate while Zhang et al. in [15] have dealt with a predator-prey system with constant harvesting rate. Especially, they reported that some spatial patterns including isolated groups such as stripelike and patch-like could be shown when the system has constant harvesting of the predator (see system (26)) whereas typical Turing patterns could emerge when the system has constant harvesting of the prey (see system (25)). Also Fasani and Rinaldi in [31] showed that Turing instability could not arise from prey harvesting:
(25)dNdT=N(1-N)-aPNN+P-h+DN∇2N,dPdT=P(-d+bPNN+P)+DP∇2P,(26)dNdT=N(1-N)-aPNN+P+DN∇2N,dPdT=P(-d+bPNN+P)-h+DP∇2P.
In this paper, we also dealt with harvesting of the prey, but the harvesting rate depends linearly on the density of the prey. We illustrated that this system has also typical Turing patterns such as spotted, stripelike, or spot-stripelike mixtures patterns which are similar phenomena to those when Wang et al. in [13] probed into a ratio-dependent predator-prey system without any harvesting rate. Moreover, we can infer from Figures 2, 3, and 4 that the presence of linear harvesting rate h for preys has dramatical effect on changing the spatial patterns from spotted patterns to the stripes pattern or coexistence of spotted pattern and stripes pattern of the prey. Biologically we have a conclusion that if the prey is harvested at a rate of its population, it has an opportunity to be centralized or assembled. From this, we could surmise that the prey has a tendency to take group defence against the harvest of it and the lower the harvesting rate is, the stronger the group defence of the prey is. Thus, we can deduce that the prey population is decreasing accordingly as the harvesting rate of the prey is increasing. On the contrary, in case of constant harvesting studied by Zhang et al. in [15], as the constant harvesting rate of the prey is increasing, the prey population is not decreasing but rather increasing (See Figure 8.) Such phenomena are out of range of general thinking that the harvesting species brings about reducing of its population.Snapshots of contour pictures of the time evolution of prey in system (25) when a=0.6, b=0.4, d=0.25,  DN=0.02, and DP=0.2: (a) h=0.148; (b) h=0.14872; (c) h=0.15.
(a)
(b)
(c)Therefore if one can monitor population of prey and predator, it is clear that the harvesting control method of prey or predator considering their populations is more realistic and reasonable than the constant harvesting control method. Thus, this research will be useful for studying the dynamic complexity of ecosystems or physical system.



---

*Source: 625973-2014-04-07.xml*