# Complex Dynamical Behavior of Holling–Tanner Predator-Prey Model with Cross-Diffusion

Authors: Caiyun Wang; Yongyong Pei; Yaqun Niu; Ruiqiang He
Journal: Complexity
 (2022)

---

## Abstract

Spatial predator-prey models have been studied by researchers for many years, because the exact distributions of the population can be well illustrated via pattern formation. In this paper, amplitude equations of a spatial Holling–Tanner predator-prey model are studied via multiple scale analysis. First, by amplitude equations, we obtain the corresponding intervals in which different kinds of patterns will be onset. Additionally, we get the conclusion that pattern transitions of the predator are induced by the increasing rate of conversion into predator biomass. Specifically, pattern transitions of the predator between distinct Turing pattern structures vary in an orderly manner: from spotted patterns to stripe patterns, and finally to black-eye patterns. Moreover, it is discovered that pattern transitions of prey can be induced by cross-diffusion; that is, patterns of prey transmit from spotted patterns to stripe patterns and finally to a mixture of spot and stripe patterns. Meanwhile, it is found that both effects of cross-diffusion and interaction between the prey and predator can lead to the complicated phenomenon of dynamics in the system of biology.


---

## Body

## 1. Introduction

Analysis of the dynamics of the predator-prey model is one of the most interesting topics in mathematics as well as in ecology. A variety of predator-prey models [1–7] that provide deep insight into the dynamical behaviors among interacting multiple species have already been investigated by many scholars in the past several decades. Pattern formation has already been analyzed in ecosystems [8–11] and epidemics for many years [12, 13]. In particular, pattern formation of the spatiotemporal predator-prey model has been avidly studied in recent decades [14–20] due to the following importance of the spatial patterns. First, the exact distributions of the population in both the spatial scale and time dimension can be well depicted by pattern formation. Second, further information on the evolution rules of individuals can be provided by spatial patterns. Finally, the influence of individual mobilities of one species on the other, such as stability and oscillatory dynamics, can be illustrated by spatial patterns.One of the predator-prey models is Leslie predator-prey model which was first introduced by Leslie [21] and takes the form(1)dNdτ=ϕNN−ψNP,dPdτ=θP1−hPN,where N and P stand for the population of prey and predator, respectively, at time τ. θ is the intrinsic growth rate of predator. h is the conversion factor of prey into predators. ϕN describes the specific rate of the prey if there is no predator. ψN is the predator function response to prey. The predator’s grow obeys logistic law, where the term N/h means the carrying capacity of predator’s environment which is proportional to the prey density. Robert May developed the Leslie predator-prey model by incorporating Holling type functional response [22, 23] to depict the predation rate. This model is known as Holling–Tanner model. Holling–Tanner model is suitable for ecological systems such as mite/spider mite, lynx/hare, sparrow/sparrow hawk, and so on [24–26]. In addition, it is assumed that the functional response is expressed by Holling type III response function [27, 28], i.e., ψN=mN2P/aN2+1, which is more suitable for the population of vertebrates. When the prey grows logistically with growth rate r and carrying capacity K in the absence of predator, i.e., ϕN=r1−N/K, we obtain the following Holling–Tanner model with Holling type III functional response:(2)dNdτ=rN1−NK−mN2PaN2+1,dPdτ=θP1−hPN,where a is the half-saturation constant. m is the maximum number of the prey that predator can capture per unit time.It is well known that spatial motion in predator–prey models includes both self-diffusion and cross-diffusion. Self-diffusion denotes the random individual mobility and indicates the movement of individuals from a higher- to lower-concentration region [29]. The concept of cross-diffusion was first proposed by Kerner [30] and was then introduced into the competitive species system [31, 32]. Taking the prey–predator model as an example, cross-diffusion terms explain the following biological meaning: predator species will move towards various directions and affect the density of different prey species at distinct places, and vice versa [33]. The mutual migration among species is typical cross-diffusion. Several phenomena cannot be well explained by only incorporating the self-diffusion term [34]; for example, there is no Turing instability in the Lotka–Volterra competitive model when only self-diffusion has been incorporated. Once the cross-diffusion term is taken into account, the above model can induce Turing instability under certain conditions. Numerous works are relevant to spatial motion [35–40]. Peng R demonstrated that cross-diffusion can induce stationary effects [41]. Oeda studied a predator-prey model with cross-diffusion and proposed a protection zone for the prey and concluded that the prey will survive because of the benefits of effects of cross-diffusion on the prey when several conditions are satisfied [42]. An attempt is made in this paper to understand the effect of cross-diffusion on the prey-predator model as well.It is clear that both intrinsic interaction between the prey and predator and the cross-diffusion affect the dynamics of the population in the predator-prey model with cross-diffusion. Our aim is to explore such effects. Moreover, most predator-prey models are nonlinear systems in order to reflect the real world [43, 44]. The nonlinear system should be approximated more accurately than in linear stability analysis. Based on Q. Ouyang’s work [45], multiple scale analysis is used to study a ratio-dependent predator-prey model with spatial motion in the present paper. Let N and P stand for the density of the prey and predator, respectively. The model is expressed as(3)∂N∂τ=rN1−NK−mN2PaN2+1+D11∇2N+D12∇2P,y,τ∈Ω×0,∞,∂P∂τ=θP1−hPN+D21∇2N+D22∇2P,y,τ∈Ω×0,∞,∂Ny,τ∂n=∂Py,τ∂n=0,y,τ∈∂Ω×0,∞,Ny,0=N0>0,y∈Ω,Py,0=P0>0,y∈Ω,where Ω represents the bounded domain in R2, and ∂Ω is the smooth bound of Ω with n being its external unit normal vector. The meanings and values of the parameters are given in Table 1.Introducing the dimensionless variablesu=N/K,v=mKP/r, and t=rτ, then system (3) can be transformed into(4)∂u∂t=u1−u−u2veu2+v+d11∇2u+d12∇2v,y,t∈Ω×0,∞,∂v∂t=vη1−γvu+d21∇2u+d22∇2v,y,t∈Ω×0,∞,∂uy,t∂n=∂vy,t∂n=0,y,t∈∂Ω×0,∞,uy,0=u0>0,y∈Ω,vy,0=v0>0,y∈Ω,where η=θ/r,e=aK2,γ=hr/mK2,d11=KD11/r,d12=rD12/m2K3,d21=mK3D21/r2, and d22=D22/mK.Table 1
The meaning of parameters in model (3).
ParameterValueMeaningrPositivePrey’s intrinsic growth ratemPositiveCapturing rateaPositiveHalf-saturation constantθPositivePredator’s intrinsic growth rateKPositiveCarrying capacityhPositiveConversion rate of biomass from prey to predatorD11PositiveDiffusion coefficient of preyD22PositiveDiffusion coefficient predatorD12NegativeCross-diffusion coefficient of preyD21PositiveCross-diffusion coefficient of predatorThe rest of this paper is organized as follows. In Section2, several conditions for the onset of the Turing instability are derived by bifurcation analysis. To show the results, an example is also provided. In Section 3, the amplitude equations near the Turing bifurcation point are obtained by virtue of weakly nonlinear analysis. To illustrate the theoretical analysis, a numerical simulation is conducted in Section 4. Moreover, pattern transitions of both the predator induced by parameter γ and prey induced by parameter d12 are obtained by simulations. Pattern transitions are obtained both by the biomass from prey to predator and by the cross-diffusion of the prey. The results reveal the complicated mutual effects of cross-diffusion and the intrinsic mechanism, e.g., biomass on the predator-prey system.


## 2. Bifurcation Analysis

We will show the Turing domain through bifurcation analysis in this section; for details of linear stability analysis of the existence, local stability, and the types of the positive equilibrium of model (4), the reader is referred to [46], while in the following, a concrete example is given to show the existence and stability of positive equilibrium. For convenience, let(5)u1−u−u2veu2+1≜fu,v,vη1−γvu≜gu,v.Without loss of generality, letE∗u∗,v∗ be the positive equilibrium of model (4). System (4) can be linearized at E∗u∗,v∗,(6)∂A^∂t=JA^+DLA^,where(7)A^=u−u∗v−v∗,J=fufvgugv,D=d11d12d21d22,L=∇2∇2∇2∇2,where(8)fu=−2eu∗3+eu∗2−1eu∗2+1,fv=−u∗2eu∗2+1,gu=ηγ,gv=−η,and obviously,(9)trJ=fu+gv,detJ=fugv−fvgu.We will give a concrete value of the equilibriumE∗u∗,v∗ as an example of the existence and stability of the equilibrium with the other values of parameters shown in Table 2.Table 2
A concrete example to illustrate the existence and stability of the equilibriumE∗u∗,v∗.
ηγeu∗v∗trJdetJ0.50.0014556040769000.25171.75−0.0260087340.1381004367Next, expand Â in Fourier space, that is, substituting(10)A^=ck1ck2eλkt+ik⋅r,into (6), where k=kx,ky, r=x,y. Following characteristic equation is obtained:(11)λk2−trkλk+Δk=0,where(12)trk=fu+gv−k2d11+d22,Δk=fugv−fvgu−k2d22fu+d11gv−d12gu+d21fv+d11d22−d12d21k4,and the eigenvalue of (11) is given as(13)λk=trk±trk2−4Δk2.Hopf bifurcation occurs when the conditionsReλk=0,Imλk≠0,trJ=0, and detJ>0, at k=0, are satisfied. When conditions such as Reλk=0 and Imλk=0, which require Δk=0 at k=kT, are satisfied, Turing bifurcation appears, and(14)kT2=fugv−fvgud11d22−d12d21.Substituting (14) into ΔK=0, the expression of γT, the value of which is the critical value for the appearance of Turing bifurcation, is obtained:(15)γT=−d11d12ηu2e−2d12d22u3e+d12d22u2e+2d11d22u2−d12d21u2−d11d12η−d12d22+22d11d122e2ηu5d22−2d123e2ηu5d21−d11d122e2ηu4d22+d123e2ηu4d21−d112d12eηu4d22+d11d122eηu4d21−2d11d12eu5d222+2d11d122eηu3d22+d11d12eu4d222−2d123eηu3d21−d122eu4d21d22+2d12eu5d21d22+d112u4d222−d11d12u4d21d22−d112d12ηu2d22+d11d122ηu2d21−d11d12u2d222+d122u2d21d22+d11d122ηd22−d123ηd211/2ηeu2+1d112e2u4η2−4d11e2ηu5d22+8d12ηu5ε2d21+4e2u6d222+2d11e2ηu4d22−4d12e2ηu4d21−4e2u5d222−2d11eηu4d21+e2u4d222−4eu5d21d22+2d112eη2u2−4d11eηu3d22+8d12eηu3d21+2eu4d21d22+4eu3d222+u4d212−2d11ηu2d21−2eu2d222+d112η2−2u2d21d22−2d11ηd22+4d12ηd21+d222.When the values of the parameters are set asη=0.8,e=130,d11=0.1,d22=8.5,d12=0.008, and d21=−0.08, let Turing domain denote the area where Turing patterns will be observed. Turing domain and the real part of λk of model (4) can be plotted as shown in Figure 1 and Figure 2, respectively.Figure 1
Turing domain of model (4).Figure 2
Dispersion relation of system (4), which shows the real part of eigenvalues λk of model (4).


## 3. Multiple Scale Analysis

To obtain the intervals of the control parameter for different types of patterns, the amplitude equations for Turing patterns of system (4) are analyzed via multiple scale analysis [45, 47] near the branch point in this section, with γ and γT being the control parameter and bifurcation threshold, respectively. Near γ=γT, the solution of system (4) is(16)uv=∑j=13AjuAjveikj⋅r+c.c…,where kj=kT, Aju, Ajv stand for the amplitude of the pattern, and c.cstands for the conjugate item.(17)∑j=13A¯juA¯jve−ikj⋅r.System (4) can be written at E∗u∗,v∗ as(18)∂X∂t=LX+N˜,where(19)L=fu+d11∇2fv+d12∇2gu+d21∇2gv+d22∇2,N˜=16f111u3+16f222v3+16f112+f121+f211u2v+16f122+f212+f221uv2+12f12+f21uv+12f11u2+12f22v216g111u3+16g222v3+16g112+g121+g211u2v+16g122+g212+g221uv2+12g12+g21uv+12g11u2+12g22v2.γ,X,t, and N then can be expanded around the threshold γT as follows. The specific expressions of h2 and h3 are(20)γT−γ=εγ1+ε2γ2+…,X=uv=εu1v1+ε2u2v2+ε3u3v3+oε4,t=T0+εT1+ε2T2,N˜=ε2h2+ε3h3+oε4,where the specific expressions of h2 and h3 are(21)h2=12f11u12+12f12+f21u1v1+12f22v1212g11u12+12g12+g21u1v1+12g22v12,h3=16f111u13+16f222v13+16f112+f121+f211u12v1+16f122+f212+f221u1v12+f11u2u1+f22v2v1+12f12+f21u2v1+u1v216g111u13+16g222v13+16g112+g121+g211u12v1+16g122+g212+g221u1v12+g11u2u1+g22v2v1+12g12+g21u2v1+u1v2.Meanwhile, the operatorL can be expanded by Taylor expansion and can be written as(22)L=LT+γT−γM,where(23)LT=fu+d11∇2fv+d12∇2gu∗+d21∇2gv+d22∇2,gu∗=ηγT,M=00ηγT20.Let(24)Aj=AjuAjv.∂Aj/∂T0=0 is satisfied, because the amplitude is a slow variable. For the derivation of the amplitude, one has(25)∂Aj∂t=ε∂Aj∂T1+ε2∂Aj∂T2+….Substituting (20)–(25) into (18) and collecting ε,ε2, and ε3, one has the linear systems ε:(26)LTu1v1=0.ε2:(27)LTu2v2=∂∂T1u1v1−d11Mu1v1−h2≜FuFv.ε3:(28)LTu3v3=∂∂T1u2v2+∂∂T2u1v1−d11Mu2v2−d22Mu1v1−h3≜HuHv,where(29)FuFv=∑j=13FujFvjeikj⋅r+c.c..(30)HuHv=∑j=13HujHvjeikj⋅r+c.c..Solving equation (26) in the light of [36] one obtains the solution of equation (26):(31)u1v1=l1W1eik1⋅r+W2eik2⋅r+W3eik3⋅r+c.c.,where l=d12kT2−fv/fu−d11kT2, c.cdenotes complex conjugate, and Wii=1,2,3 is the amplitude of the pattern for the first-order perturbation with mode eikj⋅r, the form of which is determined by higher-order terms.On the basis of the Fredholm solvability condition, (31) is substituted into equation (27) to obtain the eigenvector of the operator of LT+,(32)1l2e−ikj⋅r+c.c.,j=1,2,3,where LT+ is the eigenvector of the zero eigenvalue of the adjoint operator of LT, and(33)l2=kT2d12−fvgv−d22kT2.The orthogonality condition is given by(34)1,l2FuiFvi=0i=1,2,3,and from (34), the following equalities are given:(35)l+l2∂W1∂T1=γ1l2lηγT2W1+f11l2+f12+f21l+f22+l2g11∗l2+g12∗+g21∗l+g22∗W¯2W¯3,l+l2∂W2∂T1=γ1l2lηγT2W2+f11l2+f12+f21l+f22+l2g11∗l2+g12∗+g21∗l+g22∗W¯1W¯3,l+l2∂W3∂T1=γ1l2lηγT2W3+f11l2+f12+f21l+f22+l2g11∗l2+g12∗+g21∗l+g22∗W¯1W¯2,where gij∗ denotes the outcome of γ in gij replaced by γT .Assume that the following expression is the solution of equation (27),(36)u2v2=U0V0+∑j=13UjVjeikj⋅r+∑j=13UjjVjjei2kj⋅r+U12V12eik1−k2⋅r+U23V23eik2−k3⋅r+U31V31eik1−k1⋅r+c.c.,.Substituting (36) into (27) and collecting the coefficients of exp0,expikj⋅r,expi2kj⋅r, and expikj−ki⋅r, the following expressions are obtained:(37)U0V0=u0v0W12+W22+W32,Uj=lVj,UjjVij==u11v11Wj2,UmnVmn=umnvmnWmW¯n,where(38)u0v0=1fugv−fvgu∗⋅−gvf11l2+f12+f21l+f22+fvg11∗l2+g12∗+g21∗l+g22∗−fug11∗l2+g12∗+g21∗l+g22∗+gu∗f11l2+f12+f21l+f22,u11v11=12⋅1fu−4kT2d11gv−4kT2d22−fv−4kT2d12gu∗−4kT2d21⋅−gv−4KT2d22f11l2+f12+f21l+f22+fv−4kT2d12g11∗l2+g12∗+g21∗l+g22∗−fu−4kT2d11g11∗l2+g12∗+g21∗l+g22∗+gu∗−4kT2d21f11l2+f12+f21l+f22,umvmn=12⋅1fu−4kT2d11gv−4kT2d22−fv−4kT2d12gu∗−4kT2d21⋅−gv−4KT2d22f11l2+f12+f21l+f22+fv−4kT2d12g11∗l2+g12∗+g21∗l+g22∗−fu−4kT2d11g11∗l2+g12∗+g21∗l+g22∗+gu∗−4kT2d21f11l2+f12+f21l+f22.Using the Fredholm solvability condition, one obtains(39)1,l2HujHvj=0j=1,2,3.Equation (39) yields(40)l+l2∂V1∂T1+∂W1∂T2=ll2ηγT2γ1V1+γ2W1+A˜V¯2W¯3+V¯3W¯2+B1W12+B2W22+W32W1,l+l2∂V2∂T1+∂W2∂T2=ll2ηγT2γ1V2+γ2W2+A˜V¯1W¯3+V¯3W¯1+B1W22+B2W12+W32W2,l+l2∂V3∂T1+∂W3∂T2=ll2ηγT2γ1V3+γ2W3+A˜V¯1W¯2+V¯2W¯1+B1W32+B2W12+W22W3,where(41)A˜=lm1+n1+l2lm2+n2,B1=am1+a1n1+P12+l2am2+a1n2+P22,B2=bm1+b1n1+P1+l2bm2+b1n2+P2,a=u0+u11,b=u0+umn,a1=v0+v11,b1=v0+vmn,m1=f11l+12f12+f21,n1=f22+12f12+f21l,m2=g11∗l+12g12∗+g21∗,n2=g22∗+12g12∗+g21∗l,P1=f111l3+f112+f121+f211l2+f122+f212+f221l+f222,P2=g111∗l3+g112∗+g121∗+g211∗l2+g122∗+g212∗+g221∗l+g222,where gijk∗ denotes the outcome of γ in gijk replaced by γT.Through the above expression, together with (13) and the following formula(42)Ai=Aiu=LAiv,i=1,2,3,one has the amplitude equations for patterns of system (4):(43)τ0∂A1∂t=μA1+hA¯2A¯3−g1A12+g2A22+A32A1,τ0∂A2∂t=μA2+hA¯1A¯3−g1A22+g2A12+A32A2,τ0∂A3∂t=μA3+hA¯1A¯2−g1A32+g2A12+A22A3,where(44)τ0=l+l2γTll2η,μ=γ1T−γ1γT,h=AγTl2l2η,g1=−B1γTl3l2η,g2=−B2γTl3l2η.Based on [45] and the above analysis, we get the intervals for the pattern selection, i.e., the corresponding intervals for different kinds of pattern being observed. The results presented in Table 3 are obtained.Table 3
Pattern selection of model (4).
μPatternStabilityμ2<μ=0.0068256080<μ3Hπ spot patternUnstableμ2<μ=0.0290608556<μ3Mixture of cold spots and stripe patternsUnstableμ2<μ=0.1031783475<μ3Stripe patternUnstableμ3<μ=0.1254135951<μ4Mixture of hot spots and stripe patternsBistableμ3<μ=0.2217663346<μ4H0 spot patternUnstable


## 4. Numerical Simulations of Patterns

In this section, the value ofγ is chosen and numerical simulations are performed to verify the theoretical results. Model (4) is simulated via a finite difference approximation for spatial derivatives and the Euler method is employed. The zero-flux boundary condition with 200×200 lattice sites is employed, and the space step Δh=1 and time step Δt=0.01 are set. The simulation is run until the pattern features no longer vary. For the similar distributions on prey and predator, the predator’s pattern is chosen to illustrate details of the investigation.


### 4.1. Selection of the Predator’s Pattern

The aim here is to verify the theoretical analysis presented in Section3 and determine whether the patterns arise in their corresponding Turing domains as given in the preceding section. First, the other parameters in system (4), shown in Table 4, are fixed and the value of γ is varied. The values of parameters in amplitude equation (43) are then calculated in Table 5.Table 4
Parameters set in system (4).
ηed11d22d12d211.41000.123120.006−0.07Table 5
Parameters calculated in system (4).
γTg1g2hμ1μ2μ3μ40.013492091 7311.852514 8771.487989 955.892541 152−0.0560653277700.115719457 90.929 39Forγ=0.0134, one has μ2<μ=0.0068256080<μ3. On the basis of the results obtained in Section 3, the unstable spot pattern Hπ of the predator is obtained as shown in Figure 3(a). To show the instability of the Hπ spot pattern, by changing the value of γ slightly, namely, γ=0.0131 and γ=0.0121, where μ2<μ=0.0290608556<μ3 and μ2<μ=0.1031783475<μ3, it can be seen from Figure 3 that spot pattern Hπ turns into an unstable mixture of cold spots and stripe patterns in Figure 3(b) and stripe pattern in Figure 3(c). This mathematical analysis is consistent with numerical simulation.Figure 3
Unstable spot pattern, mixture of spots and stripes patterns, and stripe pattern of the predator. Other parameters are set as in Table4, and (a) γ=0.0134, (b) γ=0.0131, and (c) γ=0.0121.
(a)(b)(c)Forγ=0.0118,μ3<μ=0.1254135951<μ4, in light of the results presented in Table 3, a stable mixture of a hot spots and stripes pattern is obtained. Figure 4 shows the evolution of the mixture of hot spots and stripes of the predator. The numerical simulation verifies the theoretical analysis. Moreover, when γ=0.0105 is chosen, corresponding to μ3<μ=0.2217663346<μ4, the H0 spot pattern is obtained. Figure 5 shows the iterations of the hot spot pattern.Figure 4
Iterations of the mixture of hot spots and stripes of the predator withγ=0.0118; the other parameters are set as in Table 4. Iterations: (a) 0, (b) 15,000, (c) 50,000, and (d) 100,000.
(a)(b)(c)(d)Figure 5
Iterations of hot spot pattern of the predator withγ=0.0105; other parameters are set as in Table 4. Iterations: (a) 0, (b) 20,000, (c) 30,000, and (d) 100,000.
(a)(b)(c)(d)


### 4.2. Pattern Transitions of the Predator Induced by Parameterγ

In this subsection, how the spatial distribution of the predator transforms as the conversion rate of the prey into predator biomass changes is investigated. The values of the other parameters are set in Table6, and the value of γ is increased gradually. Specifically, when γ is set as γ=0.007041633659, the hot spot pattern will dominate the entire spatial region in Figure 6(a). In Figure 6(b), when γ=0.008509501103 is chosen, then the mixture of hot spots and stripes arises. When γ increases to 0.009 332 994 621, only a stripe pattern can be simulated in Figure 6(c). When γ grows further to γ=0.01032204568, the mixture of cold spots and stripes occurs as shown in Figure 6(d). Finally, when γ becomes even larger, i.e., γ=0.01092704597, the cold spot pattern shown in Figure 6(e) manifests. Figure 6 shows not only the regular changes of the pattern transition, but also that the density of the predator increases from 31 to 31.5 gradually as the γ value increases. Since γ is the amount of the prey captured by every predator, the larger the value of parameter γ, the more food the predator will eat. From the biological perspective, the more food the predator eats, the larger the predator population will be. In other words, our results are consistent with the real case in the ecosystem.Table 6
Parameters set in system (4).
ηed11d22d12d210.81300.18.50.008−0.08Figure 6
Spatial patterns of the predator induced by conversion rate of prey into predator biomass; other parameters are set as in Table5. (a) γ=0.007041633659, (b) γ=0.008509501103, (c) γ=0.009332994621, (d) γ=0.01032204568, and (e) γ=0.01092704597.
(a)(b)(c)(d)(e)


### 4.3. Pattern Transitions of Prey Induced by Parameterd12

The influence of cross-diffusion on system (4) is studied here. Several values of the parameters are given in Table 7. The value of the cross-diffusion parameter d12 is decreased. When fixing d12 as 0.01, the spot pattern shown in Figure 7(a) can be simulated. When d12 decreases to 0.006, the mixture of hot spots and stripes emerges, as shown in Figure 7(b). In Figure 7(c), a stripe pattern arises as d12 further decreases to 0.004. Finally, the entire spatial domain is filled in with cold spots and stripes when d12 is chosen only as 0.001. Figure 7 shows the regular changes of the pattern transition with the density of the prey increasing from 0.5 to 0.6 gradually as the value of d12 decreases. 12 is the cross-diffusion rate of the prey when the prey species encounters the predator. In other words, the larger the value of 12, the greater the chance that the prey will encounter the predator; that is, the prey will be at greater risk of being eaten by the predator as the prey’s population decreases. This result has a clear biological meaning.Table 7
Parameters set in system (4).
ηeγd11d22d210.85300.002048014560.212−0.074Figure 7
Spatial patterns of prey induced by cross-diffusion. (a)d12=0.01, (b) d12=0.006, (c) d12=0.004, and (d) d12=0.001.
(a)(b)(c)(d)


## 4.1. Selection of the Predator’s Pattern

The aim here is to verify the theoretical analysis presented in Section3 and determine whether the patterns arise in their corresponding Turing domains as given in the preceding section. First, the other parameters in system (4), shown in Table 4, are fixed and the value of γ is varied. The values of parameters in amplitude equation (43) are then calculated in Table 5.Table 4
Parameters set in system (4).
ηed11d22d12d211.41000.123120.006−0.07Table 5
Parameters calculated in system (4).
γTg1g2hμ1μ2μ3μ40.013492091 7311.852514 8771.487989 955.892541 152−0.0560653277700.115719457 90.929 39Forγ=0.0134, one has μ2<μ=0.0068256080<μ3. On the basis of the results obtained in Section 3, the unstable spot pattern Hπ of the predator is obtained as shown in Figure 3(a). To show the instability of the Hπ spot pattern, by changing the value of γ slightly, namely, γ=0.0131 and γ=0.0121, where μ2<μ=0.0290608556<μ3 and μ2<μ=0.1031783475<μ3, it can be seen from Figure 3 that spot pattern Hπ turns into an unstable mixture of cold spots and stripe patterns in Figure 3(b) and stripe pattern in Figure 3(c). This mathematical analysis is consistent with numerical simulation.Figure 3
Unstable spot pattern, mixture of spots and stripes patterns, and stripe pattern of the predator. Other parameters are set as in Table4, and (a) γ=0.0134, (b) γ=0.0131, and (c) γ=0.0121.
(a)(b)(c)Forγ=0.0118,μ3<μ=0.1254135951<μ4, in light of the results presented in Table 3, a stable mixture of a hot spots and stripes pattern is obtained. Figure 4 shows the evolution of the mixture of hot spots and stripes of the predator. The numerical simulation verifies the theoretical analysis. Moreover, when γ=0.0105 is chosen, corresponding to μ3<μ=0.2217663346<μ4, the H0 spot pattern is obtained. Figure 5 shows the iterations of the hot spot pattern.Figure 4
Iterations of the mixture of hot spots and stripes of the predator withγ=0.0118; the other parameters are set as in Table 4. Iterations: (a) 0, (b) 15,000, (c) 50,000, and (d) 100,000.
(a)(b)(c)(d)Figure 5
Iterations of hot spot pattern of the predator withγ=0.0105; other parameters are set as in Table 4. Iterations: (a) 0, (b) 20,000, (c) 30,000, and (d) 100,000.
(a)(b)(c)(d)


## 4.2. Pattern Transitions of the Predator Induced by Parameterγ

In this subsection, how the spatial distribution of the predator transforms as the conversion rate of the prey into predator biomass changes is investigated. The values of the other parameters are set in Table6, and the value of γ is increased gradually. Specifically, when γ is set as γ=0.007041633659, the hot spot pattern will dominate the entire spatial region in Figure 6(a). In Figure 6(b), when γ=0.008509501103 is chosen, then the mixture of hot spots and stripes arises. When γ increases to 0.009 332 994 621, only a stripe pattern can be simulated in Figure 6(c). When γ grows further to γ=0.01032204568, the mixture of cold spots and stripes occurs as shown in Figure 6(d). Finally, when γ becomes even larger, i.e., γ=0.01092704597, the cold spot pattern shown in Figure 6(e) manifests. Figure 6 shows not only the regular changes of the pattern transition, but also that the density of the predator increases from 31 to 31.5 gradually as the γ value increases. Since γ is the amount of the prey captured by every predator, the larger the value of parameter γ, the more food the predator will eat. From the biological perspective, the more food the predator eats, the larger the predator population will be. In other words, our results are consistent with the real case in the ecosystem.Table 6
Parameters set in system (4).
ηed11d22d12d210.81300.18.50.008−0.08Figure 6
Spatial patterns of the predator induced by conversion rate of prey into predator biomass; other parameters are set as in Table5. (a) γ=0.007041633659, (b) γ=0.008509501103, (c) γ=0.009332994621, (d) γ=0.01032204568, and (e) γ=0.01092704597.
(a)(b)(c)(d)(e)


## 4.3. Pattern Transitions of Prey Induced by Parameterd12

The influence of cross-diffusion on system (4) is studied here. Several values of the parameters are given in Table 7. The value of the cross-diffusion parameter d12 is decreased. When fixing d12 as 0.01, the spot pattern shown in Figure 7(a) can be simulated. When d12 decreases to 0.006, the mixture of hot spots and stripes emerges, as shown in Figure 7(b). In Figure 7(c), a stripe pattern arises as d12 further decreases to 0.004. Finally, the entire spatial domain is filled in with cold spots and stripes when d12 is chosen only as 0.001. Figure 7 shows the regular changes of the pattern transition with the density of the prey increasing from 0.5 to 0.6 gradually as the value of d12 decreases. 12 is the cross-diffusion rate of the prey when the prey species encounters the predator. In other words, the larger the value of 12, the greater the chance that the prey will encounter the predator; that is, the prey will be at greater risk of being eaten by the predator as the prey’s population decreases. This result has a clear biological meaning.Table 7
Parameters set in system (4).
ηeγd11d22d210.85300.002048014560.212−0.074Figure 7
Spatial patterns of prey induced by cross-diffusion. (a)d12=0.01, (b) d12=0.006, (c) d12=0.004, and (d) d12=0.001.
(a)(b)(c)(d)


## 5. Conclusion and Discussion

A prey-predator model of Holling–Tanner type with cross-diffusion was studied. Pattern selection and formation of the model were analyzed via weakly multiple scale analysis. First, the condition for the Turing instability of the model was obtained. Then, by choosingγ as a bifurcation parameter, the amplitude equations were deduced and the Turing pattern selection of model (4) presented. The theoretical analysis was verified by numerical results, as follows: An unstable Hπ spot pattern, unstable mixture of cold spots and stripes, and unstable stripe pattern appear, when μ satisfies μ2<μ=0.0068256080<μ3,μ2<μ=0.0290608556<μ3 and μ2<μ=0.1031783475<μ3, respectively. The stable mixture of spots and stripes arises when μ3<μ=0.1254135951<μ4 holds. Finally, it was found that the conversion of biomass from prey to predator, as well as cross-diffusion, can also result in the transition of the patterns. Compared with [9], we not only computed the interval of pattern selection mathematically, but also obtained the Turing pattern corresponding with the interval. Moreover, we studied the prey-predator model of Holling–Tanner type with cross-diffusion explosively, from aspects of pattern selection and the influence of the cross-diffusion and the intrinsic factor of conversion rate of biomass from prey to predator with rich numerical simulations.Cross-diffusion terms explain the influence of the population density of one species on the movement of other species apart from random diffusion (i.e., self-diffusion) of both species. Our numerical results about pattern transitions of prey induced by parameterd12 reveal that the more frequent the prey moves, the greater chance the prey will be at greater risk of population decreasing. This helps us further understand the cross-diffusion induced change of dynamics. Besides, together with our results about pattern transitions of the predator induced by parameter conversion rate of biomass from prey to predator, which is an intrinsic factor of prey-predator model, our study indicates that the real intrinsic mutual influence between prey and predator, along with the inherent mechanism of cross-diffusion, dominates the persistence or extinction of the entire biological system [48]. Our results are general in nature and can be used to study the effect of cross-diffusion on other prey-predator models both analytically and numerically. However, we mainly consider pattern transitions of prey induced by parameter d12. Similarly, we can study pattern transitions of prey induced by the other cross-diffusion parameter d21 and compare the differences between pattern transitions induced by d12 and d21 in our future study. Besides, we chose γ as the control parameter because γ is the parameter that stands for apparent relationship between the predator and prey. We can also choose another parameter as the control parameter to investigate the pattern transition in our future work.Multiple scale analysis around the Turing bifurcation is used to analyze the stability of different Turing patterns via the amplitude equations. Our study gets the intervals for the pattern selection. These results have been validated numerically and they matched the theoretical analysis to a large extent. Linear analysis is to analyze Turing bifurcation and Turing domain. Compared with the linear analysis, multiple scale analysis is to analyze the dynamical behavior of our model in detail (i.e., pattern selection). However, one should be aware that pattern selection is only meaningful when the value of the controlled parameter is selected to be close to the threshold ofγT. Moreover, the method of the multiple scale analysis may be applied in spatial epidemic models in the future to gain further insight into the extinction and persistence of a disease [49, 50]. Furthermore, we can take into account the imprecise value of the parameter and study the Turing patterns of an imprecise prey-predator model in future study.



---

*Source: 8238384-2022-01-10.xml*