Spatiotemporal Patterns in a Ratio-Dependent Food Chain Model with Reaction-Diffusion

and Applied Analysis 3

Food web models describe the same phenomena as predator-prey models, but the former description is more actual than the latter since our real world is so complex.Until recently, food webs models are widely studied as predatorprey models [12-14, 17, 21-29].But as far as we know, spatially extended models seem rare and not regarded.In fact, we live in a spatial world, and the spatial component of ecological interactions has been identified as an important factor in how ecological communities are shaped.Understanding the role of space is challenging both theoretically and empirically [30].And the issue of spatial and spatiotemporal pattern formation in biological communities is probably one of the most exciting problems in modern biology and ecology [31,32].
And the food web models with spatial distribution will do better job than the classical models.
In general, a classical food chain model with the nondimensional form can be written as follows: where  stands for prey density,  is the top predator density, and V-the density of the intermediate predator-describes the predator of  and the prey of ; () is the per capita rate of increase of the prey in the absence of predation.And all coefficients are positive constants,  In this paper, we focus on the following ratio-dependent food chain model [28]: ( The necessary condition of the persistence of V and  is  1 >  1 and  2 >  2 , respectively.When all the species distribute randomly in the space, model (2) can be rewritten with a supplement: where  1 ,  2 , and  3 are the diffusion coefficients of the three species, respectively, ∇ 2 = / 2 +/ 2 is the usual Laplacian operator in two-dimensional space, and other parameters have the same definitions as those above.Model (3) is to be analyzed under the nonzero initial condition and Neumann, or zero flux, boundary conditions: In the above, 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 [17].This paper is organized as follows.In the next section, we give a local stability analysis of model (3).Then, we present the pattern formation of model (3) via numerical simulations, which is followed by Section 2. Finally, we give some discussions in Section 4.

Linear Stability Analysis
There are two equilibria (steady states) in model (2), which correspond to spatially homogeneous equilibria of model (3): corresponding to top-predator extinction when  1 ( 1 − 1)/ 1 <  1 .Consider corresponding to coexistence of prey and predators when or or where V (10) And in the presence of diffusion, set  =  * + ũ, V = V * + Ṽ,  =  * + w, and the standard linear analysis predicts exponentially growing solutions of model ( 3) in the form where ⃗ ⋅ ⃗  =  2 , , and  are the wave-number and frequency, respectively.
And the eigenvalue equation then reads where the diffusion matrix  = diag( 1 ,  2 ,  3 ), and  the Jacobian matrix Then we can obtain the eigenvalues () as functions of the wave number  as the roots of where   And one type of bifurcation will break one type of symmetry of a system; that is, in the bifurcation point, two equilibrium states intersect and exchange their stability.Biologically speaking, this bifurcation corresponds to a smooth transition between equilibrium states [33].The reactiondiffusion systems have led to the characterization of two basic types of symmetry-breaking bifurcations-Hopf and Turing bifurcation, which are responsible for the emergence of spatiotemporal patterns.
The onset of Hopf instability corresponds to the case when a pair of imaginary eigenvalues cross the real axis from the negative to the positive side.And this situation occurs only when the diffusion vanishes.Mathematically speaking, the Hopf bifurcation occurs when Re(( 2 )) = 0, Im(( 2 )) ̸ = 0 at the wavenumber  = 0.For unstable steady states to heterogeneous perturbations leading to Turing patterns, the real part of the eigenvalue, Re(( 2 )), has to be greater than zero.Mathematically speaking, the Turing bifurcation occurs when I(()) = 0, R(()) = 0 at the wavenumber  ̸ = 0. Here, we take  1 as the bifurcation parameter; linear stability analysis yields the bifurcation diagram with  1 = 1.5,  2 = 2,  2 = 1,  2 = 0.5,  1 = 0.01,  2 = 0.1,  3 = 1, and  1 is a variational parameter (c.f., Figure 1).In Figure 1, the spotted curve is critical state in which above the spotted curve, the three species cannot both be positive; under the spotted curve, they are both positive.The  1 - 1 bifurcation diagram shows the two bifurcation curves separate the coexistence space into four domains.In domain I, located below all two bifurcation lines, the uniform steady state is the only stable solution of the model.Domain II is the region of pure Turing instability.Domain III is the region of pure Hopf instability.When the parameters correspond to domain IV, which is located above all two bifurcation lines, both Hopf instability and Turing instability occur.
In Figure 1, the stationary state in the parameter domains II and IV (sometimes called the "Turing space") is unstable only to a nonuniform perturbation.As expected, this domain exists only when the inhibitor species (for predator-prey system, predator ) diffuses faster than the activator species (for predator-prey system, prey V, ) and the area of this Turing space increases with  3 >  2 >  1 .

Pattern Formation
In this section, we perform extensive numerical simulations of the spatially extended model (3) in two-dimensional  3) is integrated initially in two-dimensional space from the homogeneous steady state; that is, we start with the unstable uniform solution  * 2 = ( * 2 , V * 2 ,  * 2 ) with small random perturbation superimposed; in each, the initial condition is always a small amplitude random perturbation (±5 × 10 −4 ), using an explicit Euler method for the time integration with a time stepsize of Δ = 0.01.We use the standard five-point approximation for the Laplacian operator with the Zero-flux boundary conditions and the system size is 50 × 50 space units with space stepsize (lattice constant) ℎ = Δ = Δ = 0.25, discretized through  → ( 0 ,  1 , . . .,   , . . .,  200 ) and  → ( 0 ,  1 , . . .,   , . . .,  200 ).The form of the Laplacian operator is taken as follows: The concentrations ( +1 , , V +1 , ,  +1 , ) at the moment ( + 1) at the mesh position (, ) are given by When the evolution processes reached steady state, we took a snapshot with white corresponding to the high value of prey  while black corresponding to the low one.
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 prey , for instance.
From the bifurcation diagram in the above section (cf., Figure 1), the results of numerical simulations show that the type of the system dynamics is determined by the values of  1 and  1 .And for different sets of parameters, the features of the spatial patterns become essentially different if  1 exceeds the bifurcation curves which depend on  1 .
First, we consider the pattern formation for the parameters ( 1 ,  1 ) located in domain II (c.f., Figure 1); the region of pure Turing instability occurs while Hopf stability occurs.As an example, we show the time evolution of three typical patterns when  1 = 0.6.With the parameters set, one can conclude that the critical value of Hopf bifurcation is  1 = 1.81365 and Turing bifurcation value is  1 = 2.12057.So, the values of  1 that we adopt are between 1.81365 and 2.12057.
As an example, in Figure 2, we show the time evolution of holes pattern of prey  at 0, 20000, 60000, and 200000 iterations for ( 1 ,  1 ) = (0.6, 1.84).In this case, one can see that for model (1), the pattern takes a long time to settle down, starting with a homogeneous state ( * , V * ,  * ) = (0.20267, 0.15498, 0.15498) (c.f., Figure 5(a)), and the random initial distribution leads to the formation of regular holes (c.f., Figure 2(d)).This pattern (c.f., Figure 2(d)) consists of black (minimum density of ) hexagons on a white (maximum density of ) background, that is, isolated zones with low population densities.Baurmann et al. [33] called this type pattern "cold spots" and von Hardenberg et al. [34] called it "holes." In this paper, we adopt the name "holes." When increasing  1 to 1.87, a few of stripes emerge, and the remainder of the holes pattern remains time independent (Figure 3(a)).And while increasing  1 to 1.95, model dynamics exhibits a transition from stripe-hole growth to stripes replication; that is, holes decay and the stripes pattern emerges (Figure 3(b)).
Next, we consider the pattern formation in domain IV in Figure 1; both Hopf and Turing instability occur in this domain.We adopt  1 = 0.6 and 2.12057 <  1 < 2.30769the maximized value of the coexistence of prey and their predators.The model dynamics exhibits two typical pattern formations.
In Figure 4, with the increasing of  1 to 2.1, a few of white hexagons (i.e., spots, associate with high population densities) fill in the stripes; that is, the stripes-spots pattern emerges (c.f., Figure 4(a)).And while increasing  1 to 2.21, model dynamics exhibits a transition from stripe-spots growth to spots replication; that is, stripes decay and the spots pattern emerges(c.f., Figure 4(b)).
From Figures 2-4, one can see that, with fixed parameters, on increasing the control parameter  1 , the sequence "holes → holes-stripes mixtures → stripes → spots-stripes mixtures → spots" pattern is observed.
In addition, we consider the pattern formation when ( 1 ,  1 ) locates in domain III in Figure 1, pure Hopf instability occurs.Figure 5 shows the evolution of the chaotic wave pattern of prey  at 0, 50000, 100000, and 200000 iterations with ( 1 ,  1 ) = (1.07, 5.15).With these fixed parameters, the critical value of Hopf bifurcation is  1 = 5.13108 and the Turing bifurcation values equal  1 = 5.24545.In order to make it clearer, in Figure 6, we show oscillate time series plots of , V,  (c.f., Figures 6(a), 6(b), and 6(c)), respectively.And phase portrait (c.f., Figure 6(d)) shows that there exhibits the "local" phase plane of the system obtained in a fixed point  * 2 = (0.38200, 0.05209, 0.05209) inside the region invaded by the irregular spatiotemporal oscillations.
Furthermore, we restrict our attention to the case when the top predator vanishes.Extinction of the top predator is studied by Chiu and coworkers; they gave a criterion for the extinction of top predator [35].Here, we will illustrate the pattern formation about this case.
According to food chain model,  * 1 = ( * 1 , V * 1 , 0) describes extinction of the top predator.With the same method and the same parameters in Section 2, the bifurcation diagram is shown in Figure 7.In Figure 7, the spotted curve is critical state in which the domain above the spotted curve is noncoexistence space; the domain under the spotted curve is coexistence space.Only Turing curve intersects with the spotted curve, and it separates the coexistence space into two domains.When ( 1 ,  1 ) locates in domain I, under the Turing curve, the steady state is only stable solution of model (3); when ( 1 ,  1 ) locates in domain II in Figure 7, pure Turing instability occurs.That is to say, domain II is the "Turing space" only.

Conclusions and Remarks
In summary, we have investigated a ratio-dependent spatially extended food chain model.Based on the bifurcation analysis (Hopf and Turing), we give the spatial pattern formation via numerical simulation.For the coexistence equilibrium point  * 2 = ( * 2 , V * 2 ,  * 2 ), we find that the model dynamics exhibits complex pattern replication, such as holes, holesstripes, stripes, spots-stripes, spots, and chaotic wave pattern.And for the extinction of the top predator equilibrium point  * 1 = ( * 1 , V * 1 , 0), we find that the model dynamics exhibits stripes-spots pattern replication.
In fact, in our world, every day, hundreds of species are extinct, and the extinction of a species is a fearful thing.And the top predator is extinct because there is a balance between the prey  and the intermediate predator V.In the case we considered, the density of the intermediate predator V is not small, but very big.The intermediate predator V is strong enough to fight back the top predator .
On the other hand, in the analysis of bifurcations (i.e., Hopf and Turing), we find that huge-sized computations are required, so we have to obtain more help via computers.In fact, computer-aided analysis is useful for nonlinear analysis.And computers have played an important role throughout the history of ecology.Today, numerical simulations also play an important role in spatial ecology.There are some international mathematical softwares, such as Matlab, Maple, and Mathematica, all of which have powerful function library and can provide scientific calculation and programming with friendly platform.We have finished all our symbolic computations in Maple and obtained our pattern snapshots (i.e., numerical simulations) in Matlab as Maple is more superior in symbolic computations while Matlab is more superior in numerical computations.

Figure 6 :
Figure 6: Dynamical behaviors of model (3).(a) Time-series plot of , (b) time-series plot of V, (c) time-series plot of ; (d) phase portrait.The parameters are the same as those in Figure 5.
1 is the food-independent death rate of the intermediate predator, and  2 is the food-independent death rate of the top predator.  is the functional response.The functional response is the prey consumption rate by an average single predator.It can be influenced by the prey consumption rate and the predator density.    is the amount of prey consumed per predator per unit time,     is the predator production per capita with predation.