# Pattern Formation in a Cross-Diffusive Holling-Tanner Model

Authors: Weiming Wang; Zhengguang Guo; R. K. Upadhyay; Yezhi Lin
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

We present a theoretical analysis of the processes of pattern formation that involves organisms distribution and their interaction of spatially distributed population with self- as well as cross-diffusion in a Holling-Tanner predator-prey model; the sufficient conditions for the Turing instability with zero-flux boundary conditions are obtained; Hopf and Turing bifurcation in a spatial domain is presented, too. Furthermore, we present novel numerical evidence of time evolution of patterns controlled by self- as well as cross-diffusion in the model, and find that the model dynamics exhibits a cross-diffusion controlled formation growth not only to spots, but also to strips, holes, and stripes-spots replication. And the methods and results in the present paper may be useful for the research of the pattern formation in the cross-diffusive model.


---

## Body

## 1. Introduction

A fundamental goal of theoretical ecology is to understand the interactions of individual organisms with each other and with the environment and to determine the distribution of populations and the structure of communities. Empirical evidence suggests that the spatial scale and structure of environment can influence population interactions [1]. The study of complex population dynamics is nearly as old as population ecology, starting with the pioneering work of Lotka and Volterra, a simple model of interacting species that still bears their joint names [2, 3].And the predator-prey system models such a phenomenon, pursuit-evasion, predators pursuing prey and prey escaping from the predators [4]. In other words, in nature, there is a tendency that the preys would keep away from predators and the escape velocity of the preys may be taken as proportional to the dispersive velocity of the predators. In the same manner, there is a tendency that the predators would get closer to the preys and the chase velocity of predators may be considered to be proportional to the dispersive velocity of the preys [5]. Keeping these in view, cross-diffusion arises, which was proposed first by Kerner [6] and first applied in competitive population system by Shigesada et al. [7]. From the pioneering work of Turing [8], spatially continuous models formulated as reaction-diffusion equations have been intensively used to describe spatiotemporal dynamics and to investigate mechanisms for pattern formation [9]. And the appearance and evolution of these patterns have been a focus of recent research activity across several disciplines [10].In recent years, there has been considerable interest to investigate the stability behavior of a predator-prey system by taking into account the effect of self- as well as cross-diffusion [11–25]. Cross-diffusion expresses the population fluxes of one species due to the presence of the other species.Furthermore, in [23], the authors gave a numerical study of pattern formation in the Holling-Tanner model with self- and cross-diffusion as the form
(1.1)∂N∂t=N(1-N)-aNPN+b+∇2N+D12∇2P≜f(N,P)+∇2N+D12∇2P,∂P∂t=rP(1-PN)+D21∇2N+∇2P≜g(N,P)+D21∇2N+∇2P,
with the following nonzero initial conditions:
(1.2)N(x,y,0)>0,P(x,y,0)>0,(x,y)∈Ω=[0,Lx]×[0,Ly]
and zero-flux boundary conditions:
(1.3)∂N∂n=∂P∂n=0,(x,y)∈∂Ω,
where N(t), P(t) represent population densities of prey and predator, respectively, r is the intrinsic growth rate or biotic potential of the prey, a is the maximal predator per capita consumption rate, that is, the maximum number of prey that can be eaten by a predator in each time unit, and b is the number of preys necessary to achieve one-half of the maximum rate a [26]. And the nonnegative constants, D12 and D21, called cross-diffusion coefficients, express the respective population fluxes of the prey and predators resulting from the presence of the other species, respectively. Lx and Ly give the size of the system in the directions of x and y, respectively. In (1.3), n is the outward unit normal vector of the boundary ∂Ω, which we will assume is smooth. The main reason for choosing such boundary conditions is that we are interested in the self-organization of pattern; zero-flux conditions imply no external input [11, 27].Biologically, cross diffusion implies countertransport [11]. And the induced cross-diffusion rate D12 in (1.1) represents the tendency of the prey N to keep away from its predators P and D21 represents the tendency of the predator to chase its prey. The cross-diffusion coefficients, D12 and D21, may be positive or negative. Positive cross-diffusion coefficient denotes that one species tends to move in the direction of lower concentration of another species while negative cross diffusion expresses the population fluxes of one species in the direction of higher concentration of the other species [13].In [23], the authors found that, for the equal self-diffusion coefficients (i.e., the coefficients of ∇2N and ∇2P in (1.1) are both equal to 1), only spots pattern could be obtained. In addition, they indicated that cross diffusion may have an effect on the distribution of the species, that is, it may lead the species to be isolated.There comes a question: besides spots, does model (1.1) exhibit any other pattern replication controlled by cross-diffusion?In this paper, based on the results of [23], we mainly focus on the effect of self- as well as cross-diffusion on pattern formation in the two-species Holling-Tanner predator-prey model. In the next section, we give the sufficient conditions for the Turing instability with zero-flux boundary conditions. And, by using the bifurcation theory, we give the Hopf and Turing bifurcation analysis of the model. Then, we present and discuss the results of complex, not simple, pattern formations via numerical simulations, which are followed by the last section, that is, concluding remarks.


## 2. Linear Stability Analysis

In the absence of diffusion, model (1.1) has two equilibrium solutions in the positive quadrant. One equilibrium point is given by (N0,P0)=(1,0). In the (N,P) phase plane, (1,0) is a saddle [28]. Another equilibrium point E*=(N*,P*) depends on the parameters a and b and is given by
(2.1)N*=P*=1-a-b+M2,
where M=(1-a-b)2+4b.And in the presence of diffusion, we will introduce small perturbationsU1=N-N* and U2=P-P*, where |U1|,|U2|≪1. The zero-flux boundary conditions (1.3) imply that no external input is imposed from outside. To study the effect of diffusion on the model system, we have considered the linearized form of the system as follows:
(2.2)∂U1∂t=J11U1+J12U2+∇2U1+D12∇2U2,∂U2∂t=J21U1+J22U2+D21∇2U1+∇2U2,
where
(2.3)J11=∂f∂N|(N*,P*)=(M3+(2-3a+b)M2+(1-4a+3a2-b2)M-a3+(2-b)a2+(-1-b2+2b)a-b-2b2-b3)(1-a-b+M)2,J12=∂f∂P|(N*,P*)=a(-1+a+b-M)(1-a+b+M),J21=∂g∂N|(N*,P*)=r,J22=∂g∂P|(N*,P*)=-r.SetJ=(J11J12J21J22). And define D=(1D12D211) as the diffusion matrix, and the determinant is det(D)=1-D12D21.Following Malchow et al. [29], we can know that any solution of the system (2.2) can be expanded into a Fourier series so that
(2.4)U1(r,t)=∑n,m=0∞‍unm(r,t)=∑n,m=0∞‍αnm(t)sinkr,U2(r,t)=∑n,m=0∞‍vnm(r,t)=∑n,m=0∞‍βnm(t)sinkr,
where r=(x,y) and 0<x<Lx, 0<y<Ly. k=(kn,km) and kn=nπ/Lx, km=mπ/Ly are the corresponding wavenumbers.Having substitutedunm and vnm with (2.2), we obtain
(2.5)dαnmdt=(J11-k2)αnm+(J12-D12k2)βnm,dβnmdt=(J21-D21k2)αnm+(J22-k2)βnm,
where k2=kn2+km2.A general solution of (2.5) has the form C1exp(λ1t)+C2exp(λ2t), where the constants C1 and C2 are determined by the initial conditions (1.2) and the exponents λ1 and λ2 are the eigenvalues of the following matrix:
(2.6)D~=(J11-k2J12-D12k2J21-D21k2J22-k2).Correspondingly,λ1 and λ2 arise as the solution of the following equation:
(2.7)λ2-tr(D~)λ+det(D~)=0,
where
(2.8)tr(D~)=J11+J22-2k2,det(D~)=(1-D12D21)k4+(-J11-J22+D12J21+D21J12)k2+det(J).Theorem 2.1.
If0<D12,D21<1, and b≥a/(1+a), then the uniform steady state E* of model (1.1) is globally asymptotically stable.Proof.
For global stability of nonspatial model of (1.1), we select a Liapunov function:
(2.9)V(N,P)=∫N*N‍ξ-N*ξϕ(ξ)dξ+1r∫P*P‍η-P*ηdη,
where ϕ(N)=aN/(N+b). Then
(2.10)dVdt(N,P)=N-N*Nϕ(N)dNdt+P-P*rPdPdt.
Substituting the value ofdN/dt and dP/dt from the nonspatial model of (1.1), we obtained
(2.11)dVdt=N-N*N(N(1-N)ϕ(N)-P*)-1N(P-P*)2.
Noting thatP*=(1/a)(1-N*)(N*+b), one can obtain
(2.12)dVdt=-(N-N*)2aN(N+N*+b-1)-1N(P-P*)2.
If b≥a/(1+a) holds, dV/dt≤0.
Next, we select the Liapunov function for model (1.1) (two dimensional with diffusion case):
(2.13)V2(t)=∬ΩV(N,P)dxdy,
so
(2.14)dV2dt=∬ΩdVdtdxdy+∬Ω(∇2N+D12∇2P)∂V∂Ndxdy+∬Ω(D21∇2N+∇2P)∂V∂Pdxdy.
Using Green's first identity in the plane(2.15)∬ΩF∇2Gdxdy=∫∂ΩF∂G∂nds-∬Ω(∇F·∇G)dxdy.
And considering the zero-flux boundary conditions (1.3), one can show that
(2.16)∬Ω(∇2N+D12∇2P)∂V∂Ndxdy=-∬Ω∂2V∂N2[(∂N∂x)2+(∂N∂y)2+D12(∂N∂x∂P∂x+∂N∂y∂P∂y)]dxdy,∬Ω(D21∇2N+∇2P)∂V∂Pdxdy=-∬Ω∂2V∂P2[D21(∂N∂x∂P∂x+∂N∂y∂P∂y)+(∂P∂x)2+(∂P∂y)2]dxdy,
then
(2.17)dV2dt=I1+I2,
where
(2.18)I1=∬ΩdVdtdxdy,I2=-∬Ω[∂2V∂N2((∂N∂x)2+(∂N∂y)2)+∂2V∂P2((∂P∂x)2+(∂P∂y)2)]dxdy-∬Ω(D12∂2V∂N2+D21∂2V∂P2)(∂N∂x∂P∂x+∂N∂y∂P∂y)dxdy.
Obviously, I1≤0 and I2≤0. The proof is complete.And from [19], we know that the equilibrium E* is Turing unstable if it is an asymptotically stable equilibrium of model (1.1) without self and cross diffusion but is unstable with respect to solutions of model (1.1). Hence, Turing instability sets in when the condition either tr(D~)<0 or det(D~)>0 is violated, which subject to the conditions J11+J22<0 and det(J)>0. It is evident that the condition tr(D~)<0 is not violated when the requirement J11+J22<0 is met. Hence, only violation of condition det(D~)>0 gives rise to Turing instability. Then the condition for Turing instability is given by
(2.19)H(k2)=det(D~)≡(1-D12D21)k4+(-J11-J22+D12J21+D21J12)k2+det(J)<0.Thus, a sufficient condition for Turing instability is thatH(k2)min is negative. Therefore,
(2.20)H(k2)min=det(J)-(J11+J22-D12J21-D21J12)24(1-D12D21)<0.Equation  (2.20) leads to the following final criterion for Turing instability:
(2.21)(J11+J22-D12J21-D21J12)2>4(1-D12D21)det(J).
That is to say, if H(k2)min is negative, then a range of modes k12<k2<k22 will grow Turing instability (see Figure 1).Figure 1
Modes with wavenumbers lying between the zeros,k12, and k22, of H(k2) grow in the Turing instability. Parameters: a=0.8,b=0.1,r=0.1,D12=-0.08,D21=0.8, hence, k12=0.12440, k22=0.38086, and Hmin=-0.01749.Summarizing the previous discussions, we can obtain the following theorem.Theorem 2.2.
The equilibriumE* of model (1.1) is Turing instability if J11+J22>D12J21+D21J12 and J11+J22-D12J21-D21J12>2(1-D12D21)det(J).It is easy to see that the minimum ofH(k2) occurs at k2=km2, where
(2.22)km2=J11+J22-D12J21-D21J122(1-D12D21)>0.
The critical wavenumber kc of the first perturbations to grow is found by evaluating km from (2.22).Figure2 shows that the linear stability analysis yields the bifurcation diagram with a=0.8, b=0.1, and D21=0.8. Turing and Hopf lines intersect at the Turing-Hopf bifurcation point (D12,r)=(-0.00489,0.12572) and separate the parametric space into four domains. On domain I, located above all two bifurcation lines, the steady state is the only stable solution of the model. Domain II is the region of pure Turing instability. In domain III, which is located above all two bifurcation lines, both Hopf and Turing instability occur. And domain IV is the region of pure Hopf instability.Figure 2
Bifurcation diagram for model (1.1) with parameters a=0.8,b=0.1, and D21=0.8. The Turing bifurcation line is rT=(707983156D122-(7436338235+50000A)D12+8189199110 − 62500A)/(109(D12+1)(4D122-D12-5)), where A=109(14.63-2.336D122-8.784D12). Hopf bifurcation line is rH=0.12572. Turing and Hopf lines intersect at the Turing-Hopf bifurcation point (D12,r)=(-0.00489,0.12572) and separate the parametric space into four domains.


## 3. Pattern Formation

In this section, we performed extensive numerical simulations of the spatially extended model (1.1) in 2-dimension spaces, and the qualitative results are shown here. All our numerical simulations employ the zero-flux boundary conditions with a system size of Lx×Ly, with Lx=Ly=400 discretized through x→(x0,x1,x2,…,xn) and y→(y0,y1,y2,…,yn), with n=200. Other parameters are fixed as a=0.8, b=0.1, and D21=0.8. The numerical integration of (1.1) was performed by means of forward Euler integration, with a time step of τ=0.05 and spatial resolution h=2 and using the standard five-point approximation for the 2D Laplacian with the zero-flux boundary conditions [30, 31]. More precisely, the concentrations (Ni,jn+1,Pi,jn+1) at the moment (n+1)τ at the mesh position (i,j) are given by
(3.1)Ni,jn+1=Ni,jn+τΔhNi,jn+τD12ΔhPi,jn+τf(Ni,jn,Pi,jn),Pi,jn+1=Pi,jn+τD21ΔhNi,jn+τΔhPi,jn+τg(Ni,jn,Pi,jn),
with the Laplacian defined by
(3.2)ΔhNi,jn=Ni+1,jn+Ni-1,jn+Ni,j+1n+Ni,j-1n-4Ni,jnh2.Initially, the entire system is placed in the stationary state(N*,P*)=(0.370156,0.370156), and the propagation velocity of the initial perturbation is thus on the order of 5×10-4 space units per time unit. And the system is then integrated for 100000 or 300 000 time steps and some images saved. After the initial period during which the perturbation spreads, the system goes either into a time dependent state, or to an essentially steady state (time independent).In the numerical simulations, different types of dynamics are observed, and it is found that the distributions of predator and prey are always of the same type. Consequently, we can restrict our analysis of pattern formation to one distribution. In this section, we show the distribution of preyN, for instance.Firstly, we show the pattern formation for the parameters(D12,r) located in domain II (cf. Figure 2); the region of pure Turing instability occurs while Hopf stability occurs. We have performed a large number of simulations and found that there only exhibits spots pattern in this domain. As an example, we show the time evolution of spots pattern of prey N at 0, 10000, 30000, and 300000 iteration for (D12,r)=(0.12,0.1278) in Figure 3. In this case, one can see that for model (1.1), the random initial distribution (cf. Figure 3(a)) leads to the formation of regular spots except for apparently stable defects (cf. Figure 3(d)).Spots pattern obtained with model (1.1) for (D12,r)=(0.12,0.1278). Iterations: (a) 0, (b) 10000, (c) 30000, and (d) 300000.
(a)
(b)
(c)
(d)Next, we show the pattern for the parameters(D12,r) located in domain III (cf. Figure 2), both Hopf and Turing instability occur. The model dynamics exhibits spatiotemporal complexity of pattern formation. In Figure 4, we show four typical types of patterns we obtained via numerical simulation.Four types of patterns obtained with model (1.1) for 100000 iterations. (a) Pattern α: (D12,r)=(-0.04,0.00574). (b) Pattern β: (D12,r)=(0.152,0.0882). (c) Pattern γ: (D12,r)=(0.12,0.03496). (d) Pattern δ: (D12,r)=(0.2,0.01).
(a)
(b)
(c)
(d)Patternα is time independent and consists primarily of stripes. Pattern β is time independent and consists of stripes and spots. Pattern γ is time dependent: a few of the stripes without apparent decay, but the remainder of the spots pattern remains time independent. Pattern δ is time independent and consists of black spots on a white background, that is, isolated zones with low-population densities. Baurmann et al. [32] called this type pattern “cold spots” and von Hardenberg et al. [33] called it “holes.” In this paper, we adopt the name “holes.”


## 4. Concluding Remarks

In summary, we have investigated a cross-diffusive Holling-Tanner predator-prey model with equal self-diffusive coefficients. Based on the bifurcation analysis (Hopf and Turing), we give the spatial pattern formation via numerical simulation, that is, the evolution process of the system near the coexistence equilibrium point(N*,P*).In contrast to the results in [23], we find that the model dynamics exhibits a cross-diffusion controlled formation growth not only to spots (in [23], Sun et al. claimed that the spots pattern is the only pattern of the model), but also to stripes, holes, and stripes-spots replication. That is to say, the pattern formation of the Holling-Tanner predator-prey model is not simple, but rich and complex.On the other hand, in the predator-prey model, predators will tend to gravitate toward higher concentrations of prey while prey will preferentially move toward regions where predators are rare. Models of predator-prey systems with cross diffusion have been extensively analyzed in the literature, though often with respect to their mathematical properties rather than to provide insight into the kinds of patterns that can emerge [24]. And the methods and results in the present paper may be useful for the research of the pattern formation in the cross-diffusive predator-prey model.



---

*Source: 828219-2012-12-30.xml*