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

Authors: Lei Zhang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Predator-prey models describe biological phenomena of pursuit-evasion interaction. And this interaction exists widely in the world for the necessary energy supplement of species. In this paper, 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, that is, the evolution process of the system near the coexistence equilibrium point(u2*,v2*,w2*), and find that the model dynamics exhibits complex pattern replication. For fixed parameters, on increasing the control parameter c1, the sequence “holes → holes-stripe mixtures → stripes → spots-stripe mixtures → spots” pattern is observed. And in the case of pure Hopf instability, the model exhibits chaotic wave pattern replication. Furthermore, we consider the pattern formation in the case of which the top predator is extinct, that is, the evolution process of the system near the equilibrium point (u1*,v1*,0), and find that the model dynamics exhibits stripes-spots pattern replication. Our results show that reaction-diffusion model is an appropriate tool for investigating fundamental mechanism of complex spatiotemporal dynamics. It will be useful for studying the dynamic complexity of ecosystems.


---

## Body

## 1. Introduction

Predator-prey models are studied in detail in the focus on equilibria, stability, asymptotic behavior, persistence, bifurcation, chaos, and so on [1–8]. In the past 40 years, with the idea of Turing [9], spatial extended models, in which not only the species evolve through time but also distribute in space, and pattern formation are one of the hot spots [6, 8, 10–20].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 predator-prey 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:(1)dudt=ug(u)-c1f1(u,v)v,dvdt=(m1f1(u,v)-q1)v-c2p2(v,w)w,dwdt=(m2f2(v,w)-q2)w,
where u stands for prey density, w is the top predator density, and v—the density of the intermediate predator—describes the predator of u and the prey of w; g(u) is the per capita rate of increase of the prey in the absence of predation. And all coefficients are positive constants, c1 and c2 are the maximum ingestion rates of intermediate predator and top predator, m1 is the conversion factor of prey to intermediate predator, m2 is the conversion factor of intermediate predator to top predator, q1 is the food-independent death rate of the intermediate predator, and q2 is the food-independent death rate of the top predator. fi 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. cifi is the amount of prey consumed per predator per unit time, mifi is the predator production per capita with predation.In this paper, we focus on the following ratio-dependent food chain model [28]:
(2)dudt=u(1-u)-c1uu+vv,dvdt=(m1uu+v-q1)v-c2vv+ww,dwdt=(m2vv+w-q2)w.
The necessary condition of the persistence of v and w is m1>q1 and m2>q2, respectively.When all the species distribute randomly in the space, model (2) can be rewritten with a supplement:
(3)dudt=u(1-u)-c1uu+vv+d1∇2u,dvdt=(m1uu+v-q1)v-c2vv+ww+d2∇2v,dwdt=(m2vv+w-q2)w+d3∇2w,
where d1,  d2, and d3 are the diffusion coefficients of the three species, respectively, ∇2=∂/∂x2+∂/∂y2 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:
(4)u(x,y,0)≥0,v(x,y,0)≥0,w(x,y,0)≥0,(x,y)∈Ω⊂ℝ2,∂u∂n=∂v∂n=∂v∂n=0,(x,y)∈∂Ω.
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.


## 2. Linear Stability Analysis

There are two equilibria (steady states) in model (2), which correspond to spatially homogeneous equilibria of model (3):
(5)E1*=(u1*,v1*,0)=(--m1+c1m1-c1q1m1,-(-m1+c1m1-c1q1)(m1-q1)m1q1,0),
corresponding to top-predator extinction when m1(c1-1)/c1<q1. Consider
(6)E2*=(u2*,v2*,w2*),
corresponding to coexistence of prey and predators when
(7)c1<1,m2c2-q2c2m2<q1<m2c2-q2c2m2-m1+m1c1,
or
(8)c1=1,q1<m2c2-q2c2m2,
or
(9)c1>1,m2c2-q2c2m2-m1+m1c1<q1<m2c2-q2c2m2,
where
(10)u2*=m1m2-c1m1m2+c1m2q1+c1m2c2-c1q2c2m1m2,v2*=-(m1m2-c1m1m2+c1m2q1+c1m2c2-c1q2c2)×(-m1m2+m2q1+m2c2-q2c2)×(m1m2(m2q1+m2c2-q2c2))-1,w2*=-(m1m2-c1m1m2+c1m2q1+c1m2c2-c1q2c2)×(-m1m2+m2q1+m2c2-q2c2)(m2-q2)×(m1m2(m2q1+m2c2-q2c2)q2)-1.And in the presence of diffusion, setu=u*+u~, v=v*+v~, w=w*+w~, and the standard linear analysis predicts exponentially growing solutions of model (3) in the form
(11)u~(r,t)~eλtteik→·r,v~(r,t)~eλtteik→·r,w~(r,t)~eλtteik→·r,r=(x,y),
where k→·k→=k2, k, and λ are the wave-number and frequency, respectively.And the eigenvalue equation then reads(12)|λI+k2D-J|=0,
where the diffusion matrix D=diag⁡(d1,d2,d3), and J the Jacobian matrix(13)J=[1-2u-c1vu+v+c1uv(u+v)2-c1uu+v+c1uv(u+v)20m1vu+v-m1uv(u+v)2m1uu+v-m1uv(u+v)2-q1-c2wv+w+c2vw(v+w)2-c2vv+w+c2vw(v+w)20m2wv+w-m2vw(v+w)2m2vv+w-m2vw(v+w)2-q2](u2*,v2*,w2*)=[J11J12J13J21J22J23J31J32J33].Then we can obtain the eigenvaluesλ(k) as functions of the wave number k as the roots of
(14)λ3+p(k2)λ2+q(k2)λ+r(k2)=0,
where
(15)p(k2)=(d1+d2+d3)k2-(J11+J22+J33),q(k2)=(d1d2+d1d3+d2d3)k4-(d1J22+d1J33+d2J11+d2J33+d3J11+d3J22)k2+(J11J22+J11J33-J12J21+J22J33-J23J32),r(k2)=d1d2d3k6-(d1d2J33+d1d3J22+d2d3J11)k4+(d1J22J33-d1J23J32+d2J11J33+d3J11J22-d3J12J21)k2+(-J11J22J33+J11J23J32+J12J21J33).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 reaction-diffusion 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ℜ𝔢(λ(k2))=0, ℑ𝔪(λ(k2))≠0 at the wavenumber k=0. For unstable steady states to heterogeneous perturbations leading to Turing patterns, the real part of the eigenvalue, ℜ𝔢(λ(k2)), has to be greater than zero. Mathematically speaking, the Turing bifurcation occurs when ℑ(λ(k))=0,  ℜ(λ(k))=0 at the wavenumber k≠0.Here, we takec1 as the bifurcation parameter; linear stability analysis yields the bifurcation diagram with m1=1.5, m2=2, q2=1, c2=0.5, d1=0.01, d2=0.1, d3=1, and q1 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 c1-q1 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.Figure 1
c
1
-
q
1 bifurcation diagram for model 3 with m1=1.5,  m2=2,  q2=1,  c2=0.5,  d1=0.01,  d2=0.1,  d3=1,  and  q1  a variational parameter. Hopf and Turing bifurcation curves separate the coexistence parameter space into four domains. …: the dividing line of coexistence and noncoexistence of prey and their predators.In Figure1, 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 u) diffuses faster than the activator species (for predator-prey system, prey v,w) and the area of this Turing space increases with d3>d2>d1.


## 3. Pattern Formation

In this section, we perform extensive numerical simulations of the spatially extended model (3) in two-dimensional spaces, and the qualitative results are shown here. The parameters are m1=1.5,  m2=2,  q2=1,  c2=0.5,  d1=0.01,  d2=0.1, and d3=1. Model (3) is integrated initially in two-dimensional space from the homogeneous steady state; that is, we start with the unstable uniform solution E2*=(u2*,v2*,w2*) 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 Δt=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) h=Δx=Δy=0.25, discretized through x→(x0,x1,…,xi,…,x200) and y→(y0,y1,…,yj,…,y200). The form of the Laplacian operator is taken as follows:
(16)Δhui,jn=ui+1,jn+ui-1,jn+ui,j+1n+ui,j-1n-4ui,jnh2.The concentrations(ui,jn+1,vi,jn+1,wi,jn+1) at the moment (n+1)τ at the mesh position (i,j) are given by
(17)ui,jn+1=ui,jn+τd1Δhui,jn+τf1(ui,jn,vi,jn,wi,jn),vi,jn+1=vi,jn+τd2Δhvi,jn+τf2(ui,jn,vi,jn,wi,jn),wi,jn+1=wi,jn+τd3Δhwi,jn+τf3(ui,jn,vi,jn,wi,jn).When the evolution processes reached steady state, we took a snapshot with white corresponding to the high value of preyu 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 preyu, for instance.From the bifurcation diagram in the above section (cf., Figure1), the results of numerical simulations show that the type of the system dynamics is determined by the values of c1 and q1. And for different sets of parameters, the features of the spatial patterns become essentially different if c1 exceeds the bifurcation curves which depend on q1.First, we consider the pattern formation for the parameters(q1,c1) 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 q1=0.6. With the parameters set, one can conclude that the critical value of Hopf bifurcation is c1=1.81365 and Turing bifurcation value is c1=2.12057. So, the values of c1 that we adopt are between 1.81365 and 2.12057.As an example, in Figure2, we show the time evolution of holes pattern of prey u at 0, 20000, 60000, and 200000 iterations for (q1,c1)=(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 (u*,v*,w*)=(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 u) hexagons on a white (maximum density of u) 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.”Holes pattern of the preyu obtained with model (3) with c1=1.84 and q1=0.6. Iterations: (a) 0, (b) 20000, (c) 60000, and (d) 200000.
(a)
(b)
(c)
(d)When increasingc1 to 1.87, a few of stripes emerge, and the remainder of the holes pattern remains time independent (Figure 3(a)). And while increasing c1 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)).Pattern formation of the preyu of model (3) with q1=0.6,  c1=1.87 (a) and c1=1.95 (b). Iterations: (a) 100000 and (b) 20000.
(a)
(b)Next, we consider the pattern formation in domain IV in Figure1; both Hopf and Turing instability occur in this domain. We adopt q1=0.6 and 2.12057<c1<2.30769—the maximized value of the coexistence of prey and their predators. The model dynamics exhibits two typical pattern formations.In Figure4, with the increasing of c1 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 c1 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)).Pattern formation of the preyu of model (3) with q1=0.6,  c1=2.1 (a) and c1=2.21 (b). Iterations: (a) 100000 and (b) 20000.
(a)
(b)Chaotic wave pattern of the preyu obtained with model (3) with c1=5.15 and q1=1.07. Iterations: (a) 0, (b) 50000, (c) 100000, and (d) 200000.
(a)
(b)
(c)
(d)From Figures2–4, one can see that, with fixed parameters, on increasing the control parameter c1, the sequence “holes → holes-stripes mixtures → stripes → spots-stripes mixtures → spots” pattern is observed.In addition, we consider the pattern formation when(c1,q1) locates in domain III in Figure 1, pure Hopf instability occurs. Figure 5 shows the evolution of the chaotic wave pattern of prey u at 0, 50000, 100000, and 200000 iterations with (q1,c1)=(1.07,5.15). With these fixed parameters, the critical value of Hopf bifurcation is c1=5.13108 and the Turing bifurcation values equal c1=5.24545. In order to make it clearer, in Figure 6, we show oscillate time series plots of u,v,w (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 E2*=(0.38200,0.05209,0.05209) inside the region invaded by the irregular spatiotemporal oscillations.Dynamical behaviors of model (3). (a) Time-series plot of u, (b) time-series plot of v, (c) time-series plot of w; (d) phase portrait. The parameters are the same as those in Figure 5.
(a)
(b)
(c)
(d)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,E1*=(u1*,v1*,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 (c1,q1) locates in domain I, under the Turing curve, the steady state is only stable solution of model (3); when (c1,q1) locates in domain II in Figure 7, pure Turing instability occurs. That is to say, domain II is the “Turing space” only.Figure 7
c
1
-
q
1 bifurcation diagram for model (3) with m1=1.5,  m2=2,  q2=1, c2=0.5,  d1=0.01,  d2=0.1,  d3=1 and q1; the corresponding steady state is (u1*,v1*,0)=(0.04667,0.30333,0). …: the dividing line of coexistence and noncoexistence of prey and their predators.Figure8 shows the evolution of the spatial pattern of prey at 0, 10000, 100000, and 300000 iterations with c1=1.1 and q1=0.2; that is, (c1,q1) point locates in domain II in Figure 7. The random initial distribution around the steady state E1*=(0.04667,0.30333,0) leads to the formation of stripes-spots pattern (c.f., Figure 8(d)).Spatiotemporal pattern of the preyu of model (3) with c1=1.07 and q1=0.2; the initial condition is the random initial distribution around the steady state (u1*,v1*,0)=(0.04667,0.30333,0). Iterations: (a) 0, (b) 10000, (c) 100000, and (d) 300000.
(a)
(b)
(c)
(d)


## 4. 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 pointE2*=(u2*,v2*,w2*), we find that the model dynamics exhibits complex pattern replication, such as holes, holes-stripes, stripes, spots-stripes, spots, and chaotic wave pattern. And for the extinction of the top predator equilibrium point E1*=(u1*,v1*,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 preyu 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 w.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 asMatlab, 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.



---

*Source: 130851-2014-04-17.xml*