Keep, break and breakout in food chains with two and three species

1.   Introduction

An ecosystem is a web of complex interactions among species. Food chains can be regarded as fundamental building blocks of the web. It has been studied broadly in ecological science, applied mathematics, economics and engineering science ^[1,2]. Food chain refers to the sequence of events in an ecosystem, where one organism eats others and then is eaten by another organism. A food chain links at least two species but can involve more species interacting in a more complex way. When one of the links in a food chain is no longer present the food chain will break. As examples of mechanisms for the food chain break we have a species becoming extinct ^[3,4] or when a feral animal takes over ^[5]. Sometimes, this can cause other species in the food chain to disappear or the whole ecosystem can become imbalanced or even collapse.

Since the use of mathematical models to describe the ecological system understanding of spatial and temporal behaviors of the species' population has been a central issue. One great interest for a model with multispecies' interactions is whether the involved species can keep or even stabilize at a coexistence steady state. In this case, the ordinary differential equation (ODE) system considering the species homogeneously distributed and partial differential equation (PDE) system considering the species in spatial inhomogeneous distribution are applied. The earliest ODE system of two-species food chain can date back to Lotka ^[6] and Volterra ^[7] and the relative study about coupling in predator-prey dynamics ^[8], parameterizing variable assimilation efficiency ^[9] and restructuring by recognizing prey-dependent conversion efficiency and mortality rates ^[10] are developed. Meanwhile in the tri-trophic ODE food chains it is reported that there exist rich local and global dynamical behaviors (i.e., the stability of equilibrium, local and global bifurcations, limit cycles, peak-to-peak dynamics) ^{[11,12,13,14]}. But how does the population evolve with time and space and what is the propagating behavior? It should be retrieved from the PDE systems.

A food chain is not only a relation of who eats whom but also includes the property of length of the chain ^[15,16], and can shape the community pattern ^[17,18,19]. Due to the non-uniform distribution of resources, the species will move (diffuse) to search for food and hence interact with each other. This movement of species may have an impact on the structure of the food chain. The population dynamics and stability of the interacting species concerning spatial diffusion has recently become intensive research in evolutionary game theory ^[20], replicator dynamics ^[21], and in general reaction-diffusion systems. The use of reaction-diffusion modeling interacting species can encompass diverse interests, such as spatial diffusion leading to various spatiotemporal pattern formations ^[22,23], spiral waves in complex Ginzburg-Landau equations ^[24], and cooperation patterns in prisoners' dilemma dynamics ^[25,26]. Besides these observations of spatial structure and pattern in two or three dimensions, a direct study of the effect of species' diffusion on population dynamics and wave propagation in one dimension has been developed ^{[27,28,29,30,31]}. It has been shown that propagating wavefront can exist for the competitive ^[29] and mutualistic ^[30] and even general ^[31] two-species Lotka-Volterra system with diffusion. Such a propagating wavefront represents a progressive replacement of one equilibrium (ahead of the front) by another (behind the front). So in food chains it would be of interest to find out other possible population dynamics and their properties of propagation.

In this paper, based on the traditional Lotka-Volterra predator-prey model we will consider the more realistic Rosenzweig-MacArthur model ^[32,33]. We first analyze the ODE model of two species and get the situations of cyclic coexistence and stationary coexistence. Then considering spatial diffusion we investigate the population dynamics and steady wavefront propagation in these two different coexistence situations. On one hand, we can analyze and simulate the phase behavior with time in ODE system. On the other hand, considering spatial diffusion of species we can observe the population evolution and wavefront propagation in the food chains with two or three species and get a speed diagram and then give directions to the real world.

2.   ODE model of two species

The use of mathematical models to describe the time evolution of the ecological system has been studied broadly. In the classical Lotka-Volterra predator-prey model, the number of prey is consumed and the density of the prey population is a linearly increasing function. The functional response in ecology is the intake rate of a consumer as a function of food density. It is associated with the numerical response, which is the reproduction rate of a consumer as a function of food density. Following C. S. Holling ^{[34,35,36,37]} the functional responses in the classical Lotka–Volterra predator-prey model is defined as Type Ⅰ. In Type Ⅱ the functional response is given by the following function $f(R) = \frac{{aR}}{{1 + ahR}}$, where f denotes intake rate and R denotes food (or resource) density. The rate at which the consumer encounters food items per unit of food density is called the attack rate a. The average time spent on processing a food item is called the handling time h. Rosenzweig and MacArthur ^[32] promoted a more realistic model than Lotka-Volterra model, considering the Type Ⅱ functional response instead of the Type Ⅰ functional response. Moreover this model has been broadly chosen because of its adaptability to a great variety of food chains and the richness of its behavior ^[33,37].

Our model equations are described as follows:

Here r₁, r₂ represent the growth rates of prey and predator. The environmental carrying capacity of prey is taken to be 1 and the intraspecific inhibition of prey is assumed to be 1. Considering a predator that needs time to kill and eat prey, the Holling Type Ⅱ functional response is introduced to the model, here α_{1, 2} denotes the prey's growth rate and capacity in absence of a predator, and k is the conversion efficiency of predator after consuming prey.

In this model we find when the nullclines of prey will not change monotonically and there is a symmetry axis as shown by the dotted line in Figure 1(a), which separate the phase diagram into two regions. In Figure 1(b) the intersection of the two nullclines is on the left of the symmetry axis and the trajectory of the flow closes a cycle. While in Figure 1(c) the intersection is on the right of the symmetry axis and the trajectory will flow to this point. The population evolutions of predator and prey are shown in Figures 2 and 3, respectively.

The phase diagram and population evolution with time of predator and prey can be drawn by mathematics as shown in Figures 2 and 3. According to region Ⅰ of Figure 1(a) there exists only one limit cycle, that is predator and prey can coexist and their population will change periodically when the prey is enough to bear the oscillation and satisfy the predators as shown in Figure 2, so it is called cyclic coexistence. While in region Ⅱ as shown in Figure 3 there is no limit cycle in the phase diagram and the simple closed curves will screw to a stable fixed point. And the population of predator and prey will fluctuate with time at the beginning and finally tend to the stable state, which shows a stationary coexistence. Here the initial condition is related to the population evolution, that is the oscillation of population in cyclic coexistence or fluctuation of the population in stationary coexistence will be stronger if the initial population of predator and prey are further away from the intersection point.

The phase behaviors can be further verified by theoretical analysis. There are always two fixed points (N₁*, N₂*) (in mathematics the fixed point does not change with time and space) through setting the partial to time and space is equal to 0 in Eq (1): O = (0, 0), A = (1, 0) and the third point B = ($\frac{1}{{k{\alpha _1} - {\alpha _2}}}, {\rm{ (}}1 - \frac{1}{{k{\alpha _1} - {\alpha _2}}})(1 + \frac{{{\alpha _2}}}{{k{\alpha _1} - {\alpha _2}}})/{\alpha _1}$) only appears when satisfying the physical meaning.

The stabilities of the fixed points in these ODEs (1) can be analyzed as:

The Jacobian matrix of Eq (2) is described by:

At every fixed point the eigenvector λ is determined by:

The sign of the eigenvector determines the stability of the fixed point. When λ are both negative this fixed point will be stable and both positive corresponding to an unstable fixed point.

Such as at the point O = (0, 0):

r₁, r₂ denote the growth rates of two species and we define them are both positive parameters so λ₁ is positive, λ₂ is negative and O is a saddle point.

Similarly for the fixed point A = (1, 0):

λ₁ is negative and λ₂ must be positive to satisfy the nonnegative populations, which leading A to be a saddle point.

While for the fixed point B = (N₁*, N₂*) = ($\frac{1}{{k{\alpha _1} - {\alpha _2}}}, {\rm{ (}}1 - \frac{1}{{k{\alpha _1} - {\alpha _2}}})(1 + \frac{{{\alpha _2}}}{{k{\alpha _1} - {\alpha _2}}})/{\alpha _1}$):

In Figure 2, ${r_1} = {r_2} = 1, {\alpha _1} = 1, {\alpha _2} = 2, k = 5$, B = (N₁*, N₂*) = (1/8, 35/32), $\left| {\begin{array}{*{20}{c}} {0.05}&{ - 0.1}\\ 7&0 \end{array}} \right|$, the trace is positive and λ₁, λ₂ are complex, So B is an unstable spiral. While in Figure 3 ${r_1} = {r_2} = 1, {\alpha _1} = 1, {\alpha _2} = 2, k = 2$, B = (N₁*, N₂*) = (1/2, 1), $\left| {\begin{array}{*{20}{c}} { - 0.25}&{ - 0.25}\\ {\rm{1}}&0 \end{array}} \right|$ here the trace is negative and λ₁, λ₂ are complex, So B is a stable spiral.

Compare these two different coexistences we found that in cyclic coexistence the population oscillation of predator and prey are stronger and stronger with time and the food chain will finally break if the oscillation cannot be well controlled enough, we can call this break as a temporal break in food chain. While in stationary coexistence the trajectory of the flow will fast screw to the third nonnegative fixed point and the population fluctuation also quickly approach to be calm, so it is easier to keep the food chain than in cyclic coexistence, and we regard it as temporal keep of food chain. Here we find the conversion efficiency k takes a very important role in determining the way of coexistence and the temporal keep or break of the predator-prey food chain system.

3.   PDE model of two species

When considering the limit of resources, such as water, food of prey, the mobility or migration of species cannot be ignored and it would produce important effects on the ecology of the system. On one hand the stability of the system can be changed by introducing species' diffusion. Some papers that studied stability in Lotka-Volterra systems concluded that the addition of diffusion can lead to different outcomes. On one hand Vance and Allen showed that dispersal did not always promote the stability of the population ^[38,39]. Takeuchi and Hastings suggested that diffusion does not affect the system's stability ^[40,41]. On the other hand the effects of species' mobility modeled by spatial diffusion can result in possible wavefronts giving rise to interesting spatiotemporal dynamics. Referring to our previous papers about competitive and mutualistic Lotka-Volterra system with diffusion ^[28,29,30], we know that in the presence of diffusion the ODEs can resemble the well-studied Fisher-Kolmogorov equation ^[42,43], and then Fisher's wavefront ^[44,45] and not limited to this wave profiles are observed. Such a propagating wavefront represents a progressive replacement of one equilibrium by another, including the intermediate ones ^{[27,28,29,30,31]}. Apart from stability analysis, these smooth traveling waves can also be used to provide information about population dynamics ^[46]. Here it would be of interest to find out the effect of different coexistence on the propagating behavior of this realistic predator-prey model and then study the temporal and spatial behaviors in the food chains.

The model equations with spatial diffusion are described as follows:

Firstly through numerically integrating the PDEs in Eq (3) we can observe the population evolution directly. Zero-flux Neumann boundary conditions are used. Forward Euler and upwind difference methods are applied to calculate time and space differentiation respectively. And the time is in units of 1/r₁ and space is in units of $\sqrt {{D_1}/{r_1}}$ in the numerical results. We choose two different diffusion coefficients in two coexistence regions as shown in Figures 4 and 5.

From Figure 2 we know that for cyclic coexistence there are oscillations in the population of both predator and prey and this oscillation will bring the risk for keeping the food chain with time. But it is worth to mind that considering species' diffusion this temporal oscillation of population in ODE translate to the fluctuation of the population at the beginning of the wave profiles as shown in Figure 4. Furthermore for cyclic coexistence there is still a spatial step-up oscillation at the end of the traveling wavefronts connecting from X₂, but finally the population oscillation can be controlled and more stable wave propagation is demonstrated than in ODE without diffusion, which is consistent with the demonstration of Barbera et al. ^[46]. While in stationary coexistence there is no oscillation at all except a burst before the ending of the wave profiles as shown in Figure 5, which follows the characters of the wave propagation in the competitive and mutualistic Lotka-Volterra models ^[28,29,30], where there is only stationary coexistence.

The population evolution of the prey and predator in the food chain is shown in Figures 4 and 5. In Figures 4(a) and 5(a) the predator and prey coexist in the same space and propagate with the same speed, which can be regarded as spatial keep of the food chain. This similar population dynamics has also been demonstrated in the hyperbolic reaction-diffusion model for the aquatic food chain, where considering the limit of space just the wave process at finite velocity is studied ^[46]. But here in Figures 4(b) and 5(b) the prey will perform division in space, some of the preys keep a slow wave speed to satisfy the predator and the other preys can escape from predation by propagating with a faster wavefront in the new space, which can be regarded as a spatial breakout of prey in the food chain.

The appearance of a spatial breakout of the population can also be traced by theoretical analysis. Assuming local plane wavefronts with n₁(x, t) = U₁(x −c₁t) and n₂(x, t) = U₂(x −c₂t), Eq (2) can be expressed as a four-dimensional first-order ODE dynamical system:

This new ODE system always has three fixed points X₀ ≡ (0, 0, 0, 0), X₁ ≡ (1, 0, 0, 0) and another fixed point X₂ ≡ ($\frac{1}{{k{\alpha _1} - {\alpha _2}}}, {\rm{ (}}1 - \frac{1}{{k{\alpha _1} - {\alpha _2}}})(1 + \frac{{{\alpha _2}}}{{k{\alpha _1} - {\alpha _2}}})/{\alpha _1}, {\rm{ }}0, {\rm{ }}0$) emerges when both populations are nonnegative. Nonlinear dynamical analysis is employed to analyze these fixed points one by one.

For steady state X₀ = (0, 0, 0, 0), the Jacobian matrix of eigenvalues is:

The eigenvalues of the Jacobian matrix obey $[\lambda (\lambda + \frac{{{c_1}}}{{{D_1}}}) + \frac{{{r_1}}}{{{D_1}}}][\lambda (\lambda + \frac{{{c_2}}}{{{D_2}}}) - \frac{{{r_2}}}{{{D_2}}}] = 0$. It is clear that the eigenvalues decided by species 2 are always real. For the eigenvalues related to species 1, there is a possibility to result in a spiral. But the physical requirement demands that these two eigenvalues must be real so it leads to ${c_1} \geqslant 2\sqrt {D{}_1{r_1}}$. In this case the fixed point X₀ is a saddle.

For steady state X₁ = (1, 0, 0, 0), the Jacobian matrix of eigenvalues is:

The eigenvalues obey $[\lambda (\lambda + \frac{{{c_1}}}{{{D_1}}}) - \frac{{{r_1}}}{{{D_1}}}][\lambda (\lambda + \frac{{{c_2}}}{{{D_2}}}) + \frac{{{r_2}}}{{{D_2}}}(\frac{{k{\alpha _1}}}{{1 + {\alpha _2}}} -1)] = 0$. Applying a similar analysis we get ${c_2} \geqslant 2\sqrt {D{}_2{r_2}(\frac{{k{\alpha _1}}}{{1 + {\alpha _2}}} - 1)}$. And the fixed point X₁ is also a saddle.

For steady state X₂, the eigenvalues cannot be expressed analytically and it does not set a limit to the wave speed.

From the above analysis we find traveling speed c is decided by both diffusion coefficient (D) and proliferation rate (r) of the species, it determines the steady propagating wave profiles and the population distribution.

From Figure 6 we find the phase diagram of two different wave properties is separated by the boundary c₁ = c₂. In region Ⅰ predator and prey propagate with the same speeds decided by the speed of prey c₁, it is no doubt that the two species coexist in the same space. Accordingly in the phase space the phase flow is X₂→ X₀. While in region Ⅱ the predator and prey will not propagate homogeneously. One can see the prey escapes from the control of the predator and survives in a new space without a predator. The dynamics can be described by the flow of fixed points: X₂→X₁→X₀, the dynamics of the system is first attracted to the intermediate steady state (saddle fixed point), leading to a slow wavefront of N₂ and waveback of N₁ in the wave profiles. The two different wave propagations are driven by the wave speed parameters.

Similar results have been observed in the competitive and mutualistic Lotka-Volterra model of two species ^[28,29,30], where there are no cyclic dynamics in the corresponding ODE system. But here are two different coexistences: stationary and cyclic (as shown in the ODE model) encountering the steady wave propagation. From the above results we find that on one hand diffusion can translate population oscillation and approach to stable population evolution in cyclic coexistence, on the other hand the presence of species' diffusion can lead to steady wavefront propagation and realize a spatial breakout in the food chain. Next we will mainly focus on the speed parameters and population distribution of spatial keep and breakout in the food chain with three species.

4.   PDE model of three species

The ecosystem that we wish to model is based on the Rosenzweig-MacArthur model, which assumes the prey is logistic and that the predator and top-predator have a Holling Type Ⅱfunctional response. And it is a linear three-species food chain where the lowest-level prey N₁ is preyed upon by a mid-level species N₂, which, in turn, is preyed upon by a top-level predator N₃. There are two links in this food chain: the first link is N₂ predating N₁, the second link is N₃ predating N₂. Examples of such three-species ecosystems include mouse-snake-owl and worm-robin-falcon ^[47]. Employing bifurcation analysis the dynamics of this Rosenzweig-MacArthur food chain of ODE system have recently been classified into four sets (stationary, cyclic at low frequency, cyclic at high frequency, and chaotic) ^[48].

In recent years stationary pattern in the tri-trophic food chain through inducing diffusion has been studied extensively, and many important phenomena have been observed. But direct observing the population evolution of a three-species food chain in one dimension through searching the plane wave solution is little concerned. In the above PDE model of two species we have studied cyclic and stationary states when considering the motion of species. We found diffusion can realize stable population evolution as long as the PDEs can resemble the Fisher-Kolmogorov equation and there exist steady wave profiles connecting the equilibrium fixed points although it is in an unstable cyclic state in the ODE system. So we do not focus on dynamics in ODEs with three species and we will directly study spatial keep and breakout in this more complicated food chain in PDEs. We simplify the model as follows:

Here a, b, c, d are positive parameters and the meaning is referring to the ODE model of two species. Based on the experience of the PDE model of two species here we firstly give nonlinear stability analysis and get the theoretical phase diagram and then observe the population evolution map with definite parameters.

Similarly we assume local plane wave fronts with n₁(x, t) = U₁(x−c₁t), n₂(x, t) = U₂(x−c₂t) and n₃(x, t) = U₃(x–c₃t), Eq (5) are expressed as a six-dimensional first-order ODE dynamical system:

The fixed points in Eq (6) are P₀ = (0, 0, 0, 0, 0, 0), P₁ = (1, 0, 0, 0, 0, 0) and P₂ = ($\frac{1}{{b - 1}}, {\rm{ }}[1- \frac{1}{{{{(b -1)}^2}}}]/a, {\rm{ }}0, {\rm{ }}0, {\rm{ }}0, {\rm{ }}0$) appearing when b > 2, and P₃ = ($\sqrt {1 - \frac{a}{{d - 1}}}, \frac{1}{{d - 1}}, (\frac{{b\sqrt {d - a - 1} }}{{\sqrt {d - 1} + \sqrt {d - a - 1} }} - 1)*\frac{d}{{c(d - 1)}}, {\rm{ }}0, {\rm{ }}0, {\rm{ }}0$) appearing when d > a + 1.

For steady state P₀ = (0, 0, 0, 0, 0, 0), the Jacobian matrix of eigenvalues is:

So we can get: ${c_1} \geqslant 2\sqrt {D{}_1{r_1}}$

For steady state P₁ (1, 0, 0, 0, 0, 0), the Jacobian matrix of eigenvalues is:

So we get: ${c_2} \geqslant 2\sqrt {{D_2}{r_2}(\frac{b}{2} - 1)}$.

For steady stateP₂ = ($\frac{1}{{b - 1}}, {\rm{ }}\frac{{b(b - 2)}}{{a{{(b - 1)}^2}}}, {\rm{ }}0$, 0, 0, 0), the Jacobian matrix of eigenvalues is:$\begin{array}{|cccccc|} -\lambda & 0 & 0 & 1 & 0 & 0 \\ 0 & -\lambda & 0 & 0 & 1 & 0 \\ 0 & 0 & -\lambda & 0 & 0 & 1 \\ -\frac{r_{1}}{D_{1}}\left(1-2 N_{1}-\frac{a N_{2}}{\left(1+N_{1}\right)^{2}}\right) & \frac{r_{1} N_{1} a}{D_{1}\left(1+N_{1}\right)} & 0 & -\frac{c_{1}}{D_{1}}-\lambda & 0 & 0 \\ -\frac{r_{2} N_{2} b}{D_{2}\left(1+N_{1}\right)^{2}} & \frac{r_{2}}{D_{2}}\left(\frac{b N_{1}}{1+N_{1}}-1\right) & -\frac{r_{2} c N_{2}}{1+N_{2}} & 0 & -\frac{c_{2}}{D_{2}}-\lambda & 0 \\ 0 & 0 & \frac{r_{3}}{D_{3}}\left(\frac{d N_{2}}{1+N_{2}}-1\right) & 0 & 0 & -\frac{c_{3}}{D_{3}}-\lambda \end{array}$ due to${\rm{ }}\frac{{b{N_1}}}{{1 + {N_1}}} - 1 = 0$ and$1 - 2{N_1} - \frac{{a{N_2}}}{{{{(1 + {N_1})}^2}}} = \frac{{ - 2}}{{b(b - 1)}} < 0$, so there is only one speed constraint:${c_3} \geqslant 2\sqrt {{D_3}{r_3}(\frac{{db(b - 2)}}{{a{{(b - 1)}^2} + b(b - 2)}} - 1)}$.

For steady state P₃, there is no speed limit.

Similarly the steady propagating wave profiles are determined by the speed parameters of three species, and the phase diagram of steady wave propagation is shown in Figure 7.

For a three-species food chain there are three speed limits in the wave propagation, so there are six groups of speed sorting. As shown in Figure 7 they can be divided into four regions of wave profiles. Combining with the observation of wave profiles we can further investigate the population distribution of this food chain.

Then we study the wavefront profiles and population distribution by numerically solving the PDEs in Eq (5). We set the parameter values as a = 2, b = 4, c = 3, d = 5 and get the population density in each steady state with P₀ (0, 0, 0), P₁ (1, 0, 0), P₂ (1/3, 4/9, 0), P₃ (0.707, 0.25, 0.27). Then we fix the growth rates r₁ = r₂ = r₃ = 1 and just change the diffusion coefficients according to the theoretical phase diagram to trace the wave propagation in four different regions.

4.1.   Keep of food chain (c₁ < c₂ < c₃ or c₁ < c₃ < c₂)

From Figure 8 we find that three species propagate at the same speeds and they coexist in the same space. The prey is just enough to satisfy the intermediate predator, which is also right to satisfy the top-predator. So the food chain can be kept stable forever. The flow of the dynamics can be described as P₃→P₀.

4.2.   Breakout of basal prey (c₂ < c₁ < c₃ or c₂ < c₃ < c₁)

As shown in Figure 9 we find there is a breakout in food chain. This situation originates from the division in the propagating speeds of the basal prey. The predator owns a slow speed, which results in the same slow speed as the top-predator, but the prey can escape by propagating with a faster wavefront while leaving a slower waveback of vanishing population. The waveback of the basal prey propagates with the same speed as the wavefront of intermediate predator and top-predator, reflecting the fact that they are consumed by the intermediate predator, which is also right consumed by top-predator. The dynamics of this system can be described as P₃→P₁→P₀. The faster speed of the basal preys is determined by the fixed point P₀ and the slow speed of intermediate predator and top-predator is determined by P₁. Comparing the two figures in Figure 9 we find however the top-predator travels faster or slower than the basal prey, they cannot pass the intermediate predator to connect each other. So the mid-level species is a key role to link three species in this linear food chain.

4.3.   Breakout of the group of basal prey and intermediate predator (c₃ < c₁ < c₂)

From Figure 10 we also find a breakout in the food chain, but it is different from the situation in Figure 9. Here intermediate predator and basal prey escape as a group from the top-predator. There is a simultaneous division of wave speeds both in basal prey and intermediate predator. Only the top-predator owns a slow wavefront, both basal preys and intermediate predators escape as a group to a new space with faster wavefronts and leaving wavebacks to keep the same slow speed as the top-predator. Moreover we find the waveback of the basal prey shows a decreasing population and that of intermediate predator shows an increasing population. It demonstrates that in the new space there is a new distribution of population and the intermediate predator gets a larger population. The dynamics of this system can be described as P₃→P₂→P₀. The faster speed of the basal preys is determined by the fixed point P₀ and the slow speeds of intermediate predator and basal prey are determined by P₂.

4.4.   Each breakout of basal prey and intermediate predator (c₃ < c₂ < c₁)

The most complicated wave profiles in the food chain system are shown in Figure 11. Here basal prey can live in three different regions. In the region of basal prey's slowest traveling speed where three species propagate with the same speed and the food chain keeps in space as a whole. In the middle region, the top-predator extinguishes but intermediate predator and basal prey remain and similarly intermediate predator gets a larger population. In the farthest region from the origin where basal prey travels with the fastest speed and survives freely. We find in each breakout there is a new distribution in the population. The flow of dynamics in this case can be described as P₃→P₂→P₁→P₀. The fixed points P₂, P₁ and P₀ respectively determine the values of wavefront speeds: c₃, c₂ and c₁.

5.   Conclusions

This linear food chain model can be promoted to more species. The intermediate species are both predator of lower-level species and prey of high-level species and their traveling speeds play key roles in spatial keep and breakout in the whole food chain. When a predator travels faster than prey, its speed cannot be unlimitedly promoted due to the constrain of prey and the food chain has to keep. But when the prey propagates faster than the predator, it can escape from predation and survive in the new space without a predator and realize a breakout in the food chain.

This linear food chain is only a branch in complicated food webs, so we can compare our results with the relative literature of food webs and get some extension. Borvall et al. ^[49] found that the risk of additional species extinction decreases with the number of species per functional group, which states that the more species there are in a food web, the system is less susceptible to cascading extinctions. From our paper we can indeed find that there is more probability to reorganize the population distribution so the system could be less extinct. Especially there are reports that when a top predator has been removed the risk of cascading extinction will decrease ^[50] and the decline and disappearance of the top predator will affect the distribution and abundance of the basal prey ^[51], which is consistent with our result that the lower-order species can survive in new space without predator and realize a breakout. McCann et al. ^[52] study the similar Rosenzweig-MacArthur model and discuss the dynamics of spatially coupled food webs. They also get a result that mobile consumers are important factors to the food web stability, and may display contradictory effects depending on the spatial structure. However in our model we found the mobile lower-order organisms could actively escape from predation and realize complex population dynamics, though for simplicity we did not consider the limit of space and resource here. In future work we may study the effect of time delay and spatial limit and even consider omnivory in a more realistic linear food chain in this framework.

For the Lotka-Volterra model with diffusion, many mathematicians have studied the existence and stability of traveling wave solutions and the parameter dependence of wave speed in theoretical analysis and mathematical proof ^{[42,43,44,45]}, but little is combined with numerical verification and give further discussion in the light of possible biological relevance and ecological implications. In studying the model on cell movement and growth we have calculated in a previous paper the wave profile analytically, which agrees with the numerical integrating the reaction-diffusion equations, and the wave speed can also be calculated by theory and verified by numerical measurement ^[53]. Especially this propagating speed of wavefront can be compared with experiments of Escherichia coli growth ^[54] so we can use the cell model to qualitatively explain wound healing and tumor growth in biology. Moreover in studying the hyperbolic reaction-diffusion model for the aquatic food chain, Barbera et al. also observed the similar wave propagation through numerically integrating the governing equations ^[46], and they just consider a finite velocity due to the limit of space and did not concern the wavefront connecting the intermediate state but what is the general population behavior in food chain? That is the motivation that we studied the population dynamics in this linear food chain through combining theoretical analysis and numerical measurement and we indeed observed a spatial breakout with velocity division beyond finite velocity ^[46].

Food chain contains a large number of transfers of energy or nutrients from the base to the top of a food web but they are extremely difficult to analyze. The simple model of two or three species has been broadly studied in theory and simulation, theoretical analysis mainly aims to get phase behavior and stability while simulation in two or three dimensions is to get pattern formation and transition. But little research combines ODE with PDE and performs theoretical analysis and numerical verification. In this paper we execute a systematical study on the population dynamics in Rosenzweig-MacArthur model from two species to three species, from ODE model to PDE model. Now we summarize our results in at least two aspects.

On one hand considering species' diffusion in food chain model, most studies focus on the population pattern formation and transition in two or three dimensions. We concern about the existence of the plane wave solution in this reaction-diffusion model and perform nonlinear stability analysis and then verify it through observing wave propagation in one dimension. We analyze the effect of diffusion on system stability and get insights into the population dynamics and wave propagation. Especially for unstable cyclic coexistence in ODEs, we find the presence of diffusion can translate the population oscillation and realize more stable population evolution.

On the other hand linear food chain is only one branch in the food webs, where the interactions among species are much more complicated and the decline or even extinction of one species could lead to dramatic changes in species composition and ecosystem process. In our simple model we have found the temporal keep and break of the food chain strongly depend on the transition efficiency of the predator, moreover considering species' diffusion there occurs spatial breakout in the food chain and the lower-level species has a chance to occupy new space and realize new distribution. This breakout of species may alter the population distribution and bring a new direction in food-chain research.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Acknowledgment

This work has been supported by NSFC of China under Grant No. 11204132.