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

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 crossdiffusion and interaction between the prey and predator can lead to the complicated phenomenon of dynamics in the system of biology.


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][2][3][4][5][6][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][9][10][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][15][16][17][18][19][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 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. e 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. is 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][25][26]. In addition, it is assumed that the functional response is expressed by Holling type III response function [27,28], i.e., ψ(N) � mN 2 P/aN 2 + 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) � r(1 − N/K), we obtain the following Holling-Tanner model with Holling type III functional response: 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. Selfdiffusion denotes the random individual mobility and indicates the movement of individuals from a higher-to lowerconcentration region [29]. e 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]. e 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][36][37][38][39][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]. e 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. e model is expressed as where Ω represents the bounded domain in R 2 , and zΩ is the smooth bound of Ω with n being its external unit normal vector. e meanings and values of the parameters are given in Table 1.
Introducing the dimensionless variables u � N/K, v � mKP/r, and t � rτ, then system (3) can be transformed into 2 Complexity where η � θ/r, e � aK 2 , c � hr/mK 2 , d 11 � KD 11 /r, d 12 � rD 12 /m 2 K 3 , d 21 � mK 3 D 21 /r 2 , and d 22 � D 22 /mK. e rest of this paper is organized as follows. In Section 2, 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 c and prey induced by parameter d 12 are obtained by simulations. Pattern transitions are obtained both by the biomass from prey to predator and by the crossdiffusion of the prey. e results reveal the complicated mutual effects of cross-diffusion and the intrinsic mechanism, e.g., biomass on the predator-prey system.

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 Without loss of generality, let E * (u * , v * ) be the positive equilibrium of model (4). System (4) can be linearized at where where and obviously, We will give a concrete value of the equilibrium E * (u * , v * ) as an example of the existence and stability of the equilibrium with the other values of parameters shown in Table 2.
Next, expandÂ in Fourier space, that is, substituting into (6), where k � (k x , k y ), r � (x, y). Following characteristic equation is obtained: where and the eigenvalue of (11) is given as Positive Diffusion coefficient of prey D 22 Positive Diffusion coefficient predator D 12 Negative Cross-diffusion coefficient of prey D 21 Positive Cross-diffusion coefficient of predator Complexity 3 Hopf bifurcation occurs when the conditions Re(λ k ) � 0, Im(λ k ) ≠ 0, tr(J) � 0, and det(J) > 0, at k � 0, are satisfied. When conditions such as Re(λ k ) � 0 and Im(λ k ) � 0, which require Δ k � 0 at k � k T , are satisfied, Turing bifurcation appears, and Substituting (14) into Δ K � 0, the expression of c T , the value of which is the critical value for the appearance of Turing bifurcation, is obtained: When the values of the parameters are set as η � 0.8, e � 130, d 11 � 0.1, d 22 � 8.5, d 12 � 0.008, and d 21 � −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.

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 c and c T being the control parameter and bifurcation threshold, respectively. Near c � c T , the solution of system (4) is where |k j | � k T , A u j , A v j stand for the amplitude of the pattern, and c.cstands for the conjugate item.
System (4) can be written at E * (u * , v * ) as where where the specific expressions of h 2 and h 3 are Meanwhile, the operator L can be expanded by Taylor expansion and can be written as where Let zA j /zT 0 � 0 is satisfied, because the amplitude is a slow variable. For the derivation of the amplitude, one has Substituting (20)- (25) into (18) and collecting ε, ε 2 , and ε 3 , one has the linear systems ε: where Solving equation (26) in the light of [36] one obtains the solution of equation (26): 6 Complexity W 1 e ik 1 ·r + W 2 e ik 2 ·r + W 3 e ik 3 ·r + c.c., (31) where l � d 12 k 2 T − f v /f u − d 11 k 2 T , c.cdenotes complex conjugate, and W i (i � 1, 2, 3) is the amplitude of the pattern for the first-order perturbation with mode e ik j ·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 L + T , 1 where L + T is the eigenvector of the zero eigenvalue of the adjoint operator of L T , and e orthogonality condition is given by and from (34), the following equalities are given: where g * ij denotes the outcome of c in g ij replaced by c T . Assume that the following expression is the solution of equation (27), Substituting (36) into (27) and collecting the coefficients of exp(0), exp(ik j · r), exp(i2k j · r), and exp(i(k j − k i ) · r), the following expressions are obtained:

Numerical Simulations of Patterns
In this section, the value of c 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. e 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. e 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.

Selection of the Predator's Pattern.
e aim here is to verify the theoretical analysis presented in Section 3 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 c is varied. e values of parameters in amplitude equation (43) are then calculated in Table 5.

Complexity
For c � 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 c slightly, namely, c � 0.0131 and c � 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). is mathematical analysis is consistent with numerical simulation.
For c � 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. e numerical simulation verifies the theoretical analysis. Moreover, when c � 0.0105 is chosen, corresponding to μ 3 < μ � 0.2217663346 < μ 4 , the H 0 spot pattern is obtained. Figure 5 shows the iterations of the hot spot pattern.

Pattern Transitions of the Predator Induced by Parameter
c. In this subsection, how the spatial distribution of the predator transforms as the conversion rate of the prey into predator biomass changes is investigated. e values of the other parameters are set in Table 6, and the value of c is increased gradually. Specifically, when c is set as c � 0.007041633659, the hot spot pattern will dominate the entire spatial region in Figure 6(a). In Figure 6(b), when c � 0.008509501103 is chosen, then the mixture of hot spots and stripes arises. When c increases to 0.009 332 994 621, only a stripe pattern can be simulated in Figure 6(c). When c grows further to c � 0.01032204568, the mixture of cold spots and stripes occurs as shown in Figure 6(d). Finally, when c becomes even larger, i.e., c � 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 c value increases. Since c is the amount of the prey captured by every predator, the larger the value of parameter c, 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.

Pattern Transitions of Prey Induced by Parameter d 12 .
e influence of cross-diffusion on system (4) is studied here. Several values of the parameters are given in Table 7.
e value of the cross-diffusion parameter d 12 is Table 3: Pattern selection of model (4).
μ Pattern Stability μ 2 < μ � 0.0068256080 < μ 3 H π spot pattern Unstable μ 2 < μ � 0.0290608556 < μ 3 Mixture of cold spots and stripe patterns Unstable μ 2 < μ � 0.1031783475 < μ 3 Stripe pattern Unstable μ 3 < μ � 0.1254135951 < μ 4 Mixture of hot spots and stripe patterns Bistable H 0 spot pattern Unstable   decreased. When fixing d 12 as 0.01, the spot pattern shown in Figure 7(a) can be simulated. When d 12 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 d 12 further decreases to 0.004. Finally, the entire spatial domain is filled in with cold spots and stripes when d 12 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 d 12 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. is result has a clear biological meaning.  Table 4. Iterations: (a) 0, (b) 20,000, (c) 30,000, and (d) 100,000.  Figure 6: Spatial patterns of the predator induced by conversion rate of prey into predator biomass; other parameters are set as in Table 5.

Conclusion and Discussion
A prey-predator model of Holling-Tanner type with crossdiffusion 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. en, by choosing c as a bifurcation parameter, the amplitude equations were deduced and the Turing pattern selection of model (4) presented. e 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. e 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 crossdiffusion 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 parameter d 12 reveal that the more frequent the prey moves, the greater chance the prey will be at greater risk of population decreasing. is 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 Table 7: Parameters set in system (4). 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 d 12 . Similarly, we can study pattern transitions of prey induced by the other cross-diffusion parameter d 21 and compare the differences between pattern transitions induced by d 12 and d 21 in our future study. Besides, we chose c as the control parameter because c 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.
ese 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 c 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.

Data Availability
All data and models generated or used during the study appear in the submitted article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.