Pattern formation in three species food web model in spatiotemporal domain with Beddington – DeAngelis functional response

Abstract. A mathematical model is proposed to study a three species food web model of preypredator system in spatiotemporal domain. In this model, we have included three state variables, namely, one prey and two first order predators population with Beddington-DeAngelis predation functional response. We have obtained the local stability conditions for interior equilibrium and the existence of Hopf-bifurcation with respect to the mutual interference of predator as bifurcation parameter for the temporal system. We mainly focus on spatiotemporal system and provided an analytical and numerical explanation for understanding the diffusion driven instability condition. The different types of spatial patterns with respect to different time steps and diffusion coefficients are obtained. Furthermore, the higher-order stability analysis of the spatiotemporal domain is explored.


Introduction
The foundation of population ecology was laid by animal ecologist in the first half of this century.A group of organism of the same species, which survive together in one ecological area at the same time is called a population.Within a population, all the individuals capable of reproduction have the opportunity to reproduce with other mature members of the group.Populations are always changing, hence, they are dynamic in nature [1][2][3][4][5].
The increasing population of the world could attract the attention of not only the ecologist but also of behavioral scientists.It is one of the important issues in ecology to identify some general properties about the structure (for instance, see [6] and its references).The length of food chain is one of the important features interesting for such theoretical studies.This is one of the method to estimate the energy (or a certain material) c Vilnius University, 2014 is transferred from a primary producer to a consumer.The average number of links from each producer or each top predator is regarded as the length of food chain [7].Although the network of energy flows in a food web is in general rather complex, it could be theoretically simplified to a linear chain of energy flows using the method discussed by Murray [8] and Britton [9] and along their theory, we could resolve and reconstruct the network of energy flows into a linear chain for a food web.
Many researchers have used two species model for pattern formation based on coupled reaction diffusion equations.The necessary and sufficient condition for diffusion driven instability, which leads to the formation of spatial patterns, has been derived and very interesting patterns have also been observed from the numerical simulations [10][11][12][13][14].
Recently researchers have studied the formation of patterns for different three species interacting discrete or continuous systems [1,5,[15][16][17].Most of the authors have considered a food chain model with diffusion and investigated the diffusion driven instability in the spatial system.
Keeping in view the above discussion, we will study a food web model for one prey and two predator system with Beddington-DeAngelis functional response (see Fig. 1).
The main contribution in this paper is to study the effect of diffusion on the three species food web model with Beddington-DeAngelis functional response.We have obtained analytically as well numerically the diffusion driven instability condition for the spatial system.In addition, we have obtained the different types of spatial patterns.Finally, the higher order stability analysis for the three species prey-predator system has explored.
The rest of this paper is organized as follows.In Section 2, we have proposed spatiotemporal mathematical model, Section 3 presents the local stability analysis and Hopfbifurcation for the temporal system.In Section 4, we have derived the analytical conditions for diffusion driven instability and the numerical simulation are performed.In Section 5, the higher order stability analysis is discussed.Finally, conclusion is given in last Section 6.

The mathematical model
A mathematical model of one prey utilized by two predators has been studied in temporal domain with Beddington-DeAngelis functional response [2].Motivated from the work and to make more realistic one, we extend this model in spatiotemporal domain to study it's spatial dynamics.Let U , V and W are the population densities of prey and two predators, respectively, at time T and spatial location (X, Y ).Similar as in [2], the trophic function between prey species and predator species has been described by a Beddington-DeAngelis functional response with a maximum grazing rate A i and fixed half saturation value B i , i = 1, 2, for two predators.The factors M i denotes mutual interference of predators of the same species.The parameter C i represents the conversion ratios of prey to respective predators.Moreover, the prey species grow logistically with a intrinsic growth rate R and carrying capacity K.The factors D i , i = 1, 2, are the death rates of predators species V and W , respectively.Finally, D a , D b and D c are diffusivity coefficients for prey and predators population, respectively.Thus the mathematical model governing the spatiotemporal dynamics of the three interacting species in the prey-predator community can be described by the following system of reaction-diffusion equations: with non-negative initial conditions U (0), V (0), W (0) 0. It is assumed that all parameters are positive constants.The subject to the no-flux boundary conditions and known positive initial distribution of populations are described by Here To reduce the number of parameters of the system (1)-(3), we use the following transformations: After, these substitutions we get the following system of equations (for more details, see [2]), where The positive initial distributions and no-flux boundary conditions become u(x, y, 0) > 0, v(x, y, 0) > 0, w(x, y, 0) > 0, (x, y) ∈ Ω, In equations ( 6)-( 8), Here u, v and w stand for the population densities of one prey and two predators at any instant of time t and at any point (x, y) ∈ Ω.Also the no-flux boundary conditions are used.
In the next section, we will study the dynamics of the corresponding temporal system of the spatiotemporal system ( 6)-( 8).

Analysis of temporal system
In this section, we will study the dynamical behavior of system ( 6)- (8) in the absence of diffusion, (i.e., taking diffusion coefficients D a , D b and D c equal to zero) and populations are homogeneous through out the space.Hence, the corresponding temporal system is given by Here all the parameters are strictly positive constants.There are five biologically feasible steady states: E 0 = (0, 0, 0), E 1 = (1, 0, 0), E 2 = (û, v, 0), E 2 = (ũ, 0, w) and E * = (u * , v * , w * ) for system (11)-( 13) (for more details about the steady states and stability of system, see [2]).Here we mainly focus on interior equilibrium and the unique positive u * is given by Hence, the interior equilibrium E * exists only when u * (a . The general variation matrix corresponding to system (11)-( 13) is given by where The characteristic equation for the equilibrium The constant coefficients A 1 , A 2 and A 3 can be easily calculate from the above general variation matrix.From Routh-Hurwitz criteria it follows that E * is locally stable if

Existence of Hopf-bifurcation
Now, we will study the Hopf-bifurcation of above system, taking m 1 (i.e., mutual interference of predators) as the bifurcation parameter.Again, the necessary and sufficient conditions for the existence of the Hopf-bifurcation if there exists m 1 = m 10 such that: , where u i is the real parts of the eigenvalues of the characteristic equation ( 14) of the form Now, we will verify the condition (iii) of Hopf-bifurcation.Put λ = u + iv in ( 14), we get On separating the real and imaginary parts and eliminating v between real and imaginary parts, we get Now, we have u(m 10 ) = 0 as Here differentiating (16) with respect to m 1 , then we get Now, since at m 1 = m 10 , u(m 10 ) = 0, we obtain The parametric values are 6, which ensures that the above system has a Hopf-bifurcation.It is shown graphically in Fig. 2.
2 ) > 0 of the reaction diffusion system ( 6)- (8).Obviously, the interior equilibrium point E * for the non-spatial system is a spatially homogeneous steady-state for the reaction-diffusion system ( 6)-( 8).We assume that E * is the non-spatially homogeneous equilibrium is stable with respect to spatially homogeneous following perturbation: where , η and ρ are chosen to be small and k 2 = (k 2 x + k 2 y ) is the wave number.Substituting ( 17)-( 19) into ( 6)-( 8), linearizing the system around the interior equilibrium E * , we get the characteristic equation as follows: with The eigenvalues are the solutions of the characteristic equation The coefficients p 2 , p 1 and p 0 can be found from the above matrix J k .
According to the Routh-Hurwitz criterion, all the eigenvalues have negative real parts if and only if the following conditions hold: This is best understood in terms of the invariants of the matrix and of its inverse matrix where Here matrix M ij is the adjunct of J k .We obtain the following conditions of the steady-state stability (i.e., stability for any value of k): (i) all diagonal cofactors of matrix J k must be positive; (ii) all diagonal elements of matrix J k must be negative.The two above condition taken together are sufficient to ensure stability of a give steady state.It means that instability for some k > 0 can only be observed if at least one of them is violated.Thus we arrive at the following necessary condition for the Turing instability [15,18]: (i) the largest diagonal element of matrix J k must be positive and/or (ii) the smallest diagonal cofactor of matrix J k must be negative.By the Routh-Hurwitz criteria, instability takes place if and only if one of the conditions ( 22)-( 24) is broken.We consider (23) for instability condition According to Routh-Hurwitz criterion, p 0 (k 2 ) < 0 is sufficient condition for matrix J k being unstable.Let us assume that M 33 < 0. If we choose D a = 0, D b = 0, then where k c > 0 is the small critical wave number.Hence, in this system, diffusion-driven instability occurs, i.e., the small spatio-temporal perturbations around the homogeneous steady-state are unstable and, hence, of generation of spatio-temporal pattern is justified.Now, we obtained the eigenvalues of the characteristic equation ( 21) numerically of the spatial system (6)-( 8).Here we choose some parametric values of a 1 = 0.

Spatiotemporal pattern formation
It is well known that the analytical solution of the coupled reaction diffusion system is not always possible.Hence, one has to use numerical simulations to solve them.The spatiotemporal system ( 6)-( 8) is solved numerically in two-dimensional space using a finite difference method for the spatial derivatives.In order to avoid numerical artifacts, the values of the time and space steps have been chosen sufficiently small.For the numerical simulations, the initial distributions of the species are considered as small spatial perturbation of the uniform equilibrium.All the numerical simulations use the zero-flux boundary condition in a square habitat of size 200 × 200 and 100 × 100.The iterations are performed for different step sizes in time.
The spatial distributions of prey-predator system in the time evaluation are given in Figs.4-7.By varying coupling parameters, we observed that if one parameter value   changes, then spatial structure changes over the times of the spatial system.In Figs.4-7, we observed well organized structures for the spatial distribution of populations, and also when time T increase, the densities of different classes of population become uniform throughout the space.The parametric values are used same as in Fig. 3.In Figs. 4 and 6, we have shown the pattern formation with different time steps.It can be observed that a stationary "mixtures → stripe-spot mixtures → spots" patterns are time-dependent, as similar in [19].From Figs. 5 and 7, we observed that the pattern sequence are mixtures of "spots → stripe" and "spots-stripe → spots", respectively, as in [20][21][22][23].Finally, From Fig. 7, it is observed that the higher diffusivity coefficients stabilized the spatial system.

Higher order stability analysis
In this subsection, we will determine the instability conditions by the higher-order spatiotemporal perturbation terms [24].We choose a general two non-dimensional reactiondiffusion system.System ( 6)-( 8) is recalled with specific choice of parameter values.The three dimensional reaction diffusion systems are described as follows: with no-flux boundary conditions and initial distribution of population within 2D bounded domain.The interior equilibrium point E * for the non-spatial system corresponding to the system is a spatially homogeneous equilibrium for system (27)-(29).We consider E * is locally asymptotically stable equilibrium for the temporal model.Taking the spatial perturbations u(t, x, y), v(t, x, y) and w(t, x, y) on the steady states u * , v * , w * defined by u = u * + n(t, x, y), v = v * + p(t, x, y), w = w * + m(t, x, y) and then expanding the temporal part in Taylor series up to second order around the steady state, we find following three expressions: Now, taking spatial perturbation in the form of with no-flux boundary condition leads to the following three system of equations: It is clear from above three equations that the growth or decay of first-order perturbation terms depends upon the second-order perturbation terms.Further, we need the dynamical equations for second-order perturbation terms involved in (33)-(35).Multiplying each term of equation ( 33) by 2u and neglecting the contribution of third-order perturbation terms, we find the dynamical equation for u 2 as and proceeding in a similar fashion, the dynamical equations for remaining second-order perturbations are given by The truncation of third and higher-order terms in Taylor series expansion and neglecting of third and higher-order perturbation terms during derivation of dynamical equations ( 33)-(41) leads us to a closed system of equations for n, p, m, n 2 , p 2 , m 2 , np, pm, nm.
Otherwise, one cannot avoid infinite hierarchy of dynamical equations for perturbation terms.Truncation of higher-order terms does not affect the understanding of the role of leading-order non-linearity.Applicability and significance of the analysis can be justified with the perturbation terms up to order three for system ( 6)-( 8) with suitable choice of parameter values.Consideration of third-and higher-order perturbation terms may be required for this type of analysis use in other system.It also depends upon the nonlinearity involved.The dynamical equations ( 33)-( 41) can be written into a compact matrix form as follows: where X = [n, p, m, n 2 , p 2 , m 2 , np, pm, nm] T and   6)-( 8) around the steady state (u * , v * , w * ) = (0.0783169, 0.51213, 0.209704).We found that one eigenvalue having positive real part for a range of values of k in Fig. 8, we have plotted largest Re{λ(k)} defined as linear obtained by solving (21) along with largest Re{λ(k)} defined as higher order computed numerically for the characteristic equation (43) for a range of wavelengths.It is clear that linear and higher order are positive for k ∈ (0, 0.5) over the entire range (see Fig. 8).

Conclusion
In this paper, a one prey and two predators system with Beddington-DeAngelis functional response is considered.It is shown that there exists Hopf-bifurcation with respect to mutual interference of predator.In the qualitative analysis, we studied dynamical behavior of the temporal system.It is established that when the rate of mutual interference of predator, i.e., M 1 , crosses its threshold value, i.e., M 1 = M 10 , then prey, first predator and Nonlinear Anal.Model.Control, 2014, Vol. 19, No. 2, 155-171 second predator populations start oscillating around the interior equilibrium.The above result has been shown numerically in Fig. 2 for different values of M 1 .In particular, in Fig. 2(a), we observe that the interior equilibrium is stable, when M 1 = 0.95, but when it crosses the threshold value of M 1 = 0.95, the above system shows Hopf-bifurcation, as shown in Figs.2(b)-(d).Furthermore, we have observed that the diffusion instability occurs in the spatial system (see Fig. 3).We have observed the nature of spatial patterns with respect to time (see Figs. 4-6).Also, we observed that if a diffusivity coefficient increases, then the population densities become uniform and spotted pattern observed (see Fig. 7).All these spatial patterns show that the qualitative changes lead to spatial density distribution of the spatial system for each species.Furthermore, we have analyzed the stability of linear and non-linear systems with the help of higher order stability analysis and also observed that the linear and non-linear systems change the behavior of the system from instability to stability (see Fig. 8).Our result shows that the modeling by reactiondiffusion equation is an appropriate tool for investigating fundamental mechanisms of spatiotemporal dynamics in the real world food web system.

Fig. 2 .
Fig. 2. The phase representation of three species around the interior equilibrium.

Fig. 3 .
Fig. 3. Plot of max Re(λ(k)) against k.The other parametric values are given in text.

Fig. 8 .
Fig. 8. Plot of maximum value of Re(λ(k)) versus k: black line for linear and red line for non-linear systems.The parametric values are given in the text.(Online version in colour.) 2is the Laplacian operator in two-dimensional cartesian coordinate system, Ω is 2D bounded rectangular domain with boundary ∂Ω, ∂/∂n is the outward drawn normal derivative on the boundary, D a , D b , D c are positive constant diffusion coefficients for one prey and two predators population respectively.
u, v, w are new scaled measures of population size; D a , D b , D c are