 Review
 Open Access
 Published:
Numerical simulations of multicomponent ecological models with adaptive methods
Theoretical Biology and Medical Modellingvolume 13, Article number: 1 (2016)
Abstract
Background
The study of dynamic relationship between a multispecies models has gained a huge amount of scientific interest over the years and will continue to maintain its dominance in both ecology and mathematical ecology in the years to come due to its practical relevance and universal existence. Some of its emergence phenomena include spatiotemporal patterns, oscillating solutions, multiple steady states and spatial pattern formation.
Methods
Many timedependent partial differential equations are found combining loworder nonlinear with higherorder linear terms. In attempt to obtain a reliable results of such problems, it is desirable to use higherorder methods in both space and time. Most computations heretofore are restricted to second order in time due to some difficulties introduced by the combination of stiffness and nonlinearity. Hence, the dynamics of a reactiondiffusion models considered in this paper permit the use of two classic mathematical ideas. As a result, we introduce higher order finite difference approximation for the spatial discretization, and advance the resulting system of ODE with a family of exponential time differencing schemes. We present the stability properties of these methods along with the extensive numerical simulations for a number of multispecies models.
Results
When the diffusivity is small many of the models considered in this paper are found to exhibit a form of localized spatiotemporal patterns. Such patterns are correctly captured in the local analysis of the model equations. An extended 2D results that are in agreement with Turing typical patterns such as stripes and spots, as well as irregular snakelike structures are presented. We finally show that the designed schemes are dynamically consistent.
Conclusion
The dynamic complexities of some ecological models are studied by considering their linear stability analysis. Based on the choices of parameters in transforming the system into a dimensionless form, we were able to obtain a wellbalanced system that is biologically meaningful. The accuracy and reliability of the schemes are justified via the computational results presented for each of the diffusive multispecies models.
Background
The study of reactiondiffusion problems in ecological context have gained a huge amount of scientific interest, due to their practical relevance and emergence of some interesting phenomena such as spatial patterns, oscillating solutions, phase planes, chaotic behaviours and multiple steady states to mention a few. The most popular and wellknown predatorprey model is named after the two scientists, Alfred Lotka (1880–1949) and Vito Volterra (1860–1940). Lotka and Volterra in their earlier research work apply the model to address the interacting population systems called predator–prey. Our numerical study in this paper is aimed at reflecting the types of interactions which we describe as predation (a process where one species of organisms called predator depends solely on the other known as the prey, for survival), competition (a situation whereby two or more different species of organisms struggle for the available resources; definitely, we expect the growth rate of each population to decrease) and lastly, the mutualism or symbiosis (organisms coexist without negatively affecting each other and hence, species growth rate is increased) [2, 11, 12, 23, 32, 33, 37, 39, 43, 48, 49].
A lot of research attention has been devoted to the study of population dynamics with regards to ecological interactions over the past few decades. Such studies include the predatorprey system that describes the situation in which the existence of the species called the predator depends solely on the other species called the prey. The predatorprey system has received tremendous attraction over the years but has been represented mainly in terms of ordinary differential equations, which modelled the species distribution. The Dynamics of the LotkaVolterra predatorprey model are quite interesting. However, this model is structurally unstable since a small perturbation of the equations often results to a drastic change in the dynamical system. For this reason, the presence of diffusion mechanism takes place though it changes the behavior of the whole model to coupled partial differential equations, called reactiondiffusion system. With the introduction of diffusion, the analysis of the whole system remain tactical in the literature [46, 47], and therefore, numerical approximations are quite often used to simulate these types of models.
Predatorprey systems have been studied by many researchers in various forms. For instance, in bacteria ecology, computer simulations of complex spatiotemporal patterns [4, 11] of Bacillus subtilis based on stochastic models [22] and deterministic models [29], Allee effect of patchy invasion on predatorprey dynamics [1, 3, 5, 7, 13, 39]. The diffusive predatorprey systems have also been studied extensively, see, [11, 17, 26, 27, 31, 38, 43]. Moreover, Wang et al. [46] investigated the spatial pattern formation of a predatorprey system with preydependent functional response of Ivlev type and reactiondiffusion whereas the analysis of predatorprey systems showing the Holling type II functional response is examined in Garvie and Trenchea [12].
Another type of interspecies interaction is given by the competition. Competitive species models or community models further describe a situation where consumers share some resources that can affect their rate of production. Many ecologist, however, put greater weight on competition which was thought to play a predominant role over the years in structuring ecological communities. Notwithstanding, there is a classical model of competition due to Lotka [24, 25] and Volterra [44, 45]. The LotkaVolterra competition model is an interference model where two species are assumed to diminish each other's per capita growth rate by direct interference. It is usually assumed in this model that each species has a different population of different sizes that grow logistically in absence of each other and that each has a per capita growth rate that decreased linearly with the population size with their own intrinsic growth rate and carrying capacity. Mathematically, the simplest and instructive case is described by a system of two coupledreaction diffusion equations. The system of two competing species in just onedimensional space has attracted a lot of attentions, see [10, 14, 48] for details. Some of the evolution processes here are characterized owing to the fact that certain moments of time they experience a sudden change of state. To this end, we additionally consider a general case of $n$ competing species that is less investigated and still poorly understood for case n ≥ 2. Among the few works done when n > 2 include [37, 40, 43].
In mutualistic systems, organisms are found to evolve together. The existence of one has no negative effect on the other, each is part of the other's environment and coexist, and they make use of each other in such a way that both organisms are benefited. Mutualism has not been given as much attention as predation and competition. Readers are referred to [20, 23] for a thorough review of the natural history of mutualism. Community invasion models have an issue of significant importance in the contemporary study of biological and ecological systems which have drawn the attention of both theorists and ecologists since the foundation work of Holt [18]. Despite a considerable achievements recorded in the field of population dynamics modeling the interaction of a multispecies community, so many challenging and open problems that are of great ecological importance are yet to be addressed.
Mathematical analysis of the main equations
In this work, our major attention is on the twovariable reactiondiffusion systems. We shall adapt linear stability analysis method to discuss the general two species dynamics. Let u and v be the variables representing the two species of the LotkaVolterra predatorprey type. In the convention here, v is the predator, while u represents the prey.
The most relevant and general twospecies reactiondiffusion system is formulated as
subject to zeroflux boundary conditions on a closed interval, say [0, L]
We assume that the point \( \left(\widehat{u},\widehat{v}\right) \) is stable equilibrium state of the homogeneous system
that is \( f\left(\widehat{u},\widehat{v}\right)=0,\kern0.48em g\left(\widehat{u},\widehat{v}\right)=0 \). Stability of the steady states for general twovariable system can be represented by the Jacobian
Which leads to characteristic equation of the form λ ^{2} − trJλ + det J = 0, where
To examine the stability of the uniform steady state \( \left(\widehat{u}(x),\widehat{v}(x)\right)=\left(\widehat{u},\widehat{v}\right) \), we carry out the linear stability analysis in the spirit of Allen [2], Mendez et al. [28] and Murray [32, 33], we obtain
where λ _{ k }, the growth rates and the modes cos(kx) are the roots of polynomial
which corresponds to a polynomial
representing the dispersion relation, with
Known from the stability conditions in (5) that trJ < 0, thus
Which shows that the uniform steady state of (1) cannot undergo an oscillatory instability (or wave bifurcation) to a standing wave pattern.
A Turing instability corresponds to \( {\lambda}_{k_{trJ}}=0 \) for k _{ trJ } ≠ 0. That is, with Φ _{2} = 0, results to (J _{11} − D _{ u } k ^{2})(J _{22} − D _{ v } k ^{2}) − J _{12} J _{21} = 0 or
For the roots of (9) to be positive,
is a necessary but not sufficient condition for the Turing instability to occur. With reference to conditions in (5), Turing instability can occur if the diffusion coefficient D _{ u } ≠ D _{ v } and if the matrix elements J _{11} and J _{22} have opposite sign. So, Turing instabilities occur only in either pure or cross activatorinhibitor dynamical system.
The system (3) is of the pure LotkaVolterra type if the Jacobian agrees with the structure of the form
Again, it is of the cross LotkaVolterra type if the Jacobian has the structure
Clearly from systems (11) or (12), we have J _{11} > 0, J _{22} < 0 which together with D _{ v } J _{11} + D _{ u } J _{22} > 0 indicates that Turing instability can occur only if J _{22} > J _{11} since trJ < 0 and J _{12} J _{21} < 0 for det J > 0. If we let Θ _{ D } = D _{ v }/D _{ u } be the diffusion coefficients ratio, we can easily obtain from (10) that Θ _{ D } > J _{22}/J _{11} > 1. The indication here is that, for Turing instability to take place, the inhibitor must diffuse faster than the activator.
By rewriting (9) in the form
where
we have the roots of equation (9) given by
provided Ψ _{1} < 0 and condition (10) is satisfied. So, Turing instability occurs for (10) to have a double root, that is, if Ψ _{1} ^{2} − 4Ψ _{2} = 0.
In conclusion, the uniform steady state of the reactiondiffusion system (1), \( \left(\widehat{u}(x),\widehat{v}(x)\right)=\left(\widehat{u},\widehat{v}\right) \), that satisfy the stability conditions in (5) will be unstable in the presence of diffusion (called diffusion driveninstability) if Ψ _{1} < 0, that is, D _{ v } J _{11} + D _{ u } J _{22} > 0, and Ψ _{1} ^{2} − 4Ψ _{2} > 0, that is, (D _{ v } J _{11} + D _{ u } J _{22})^{2} > 4D _{ u } D _{ v } det J, with the band of unstable modes
A good research focus has to be given to the numerical simulation of multispecies dynamics in more than one dimensional space which has received little attention in the literature. We simulated a class of biological systems that lead to the evolution of traveling waves and formation of chaotic and spatiotemporal patterns arising in the context of mathematical ecology. Though, simulations that are based on the use of conventional methods in twodimensions are found to be time consuming. As a result, consideration is given to the design and method of implementing a viable numerical scheme that can handle a class of multicomponent reactiondiffusion problems efficiently.
Numerical methods
Many systems of nonlinear time dependent reactiondiffusion problems of partial differential equations that are of physical interest are written in the compact form
where D > 0 is the diffusion coefficient, ∇^{2} w = (∂^{2} w/∂x ^{2} + ∂^{2} w/∂y ^{2}) is the twodimensional Laplacian operator that represents the linear term, and N is nonlinear function of w and t. By following [34–36], we discretize the spatial domain by mesh (x _{ i }, y _{ j }) = (T _{1} + i × h _{ x }, T _{1} + j × h _{ y }) where h _{ x } = (T _{2} − T _{1})/(N _{ x } + 1), h _{ y } = (T _{2} − T _{1})/(N _{ y } + 1), for 0 ≤ i ≤ N _{ x } + 1, 0 ≤ j ≤ N _{ y } + 1. We approximate the secondorder derivatives by the fourthorder central difference operators
and
for i = 1, 2, …, N _{ x } and j = 1, 2, …, N _{ y } + 1. The discretized form of (13) lead to a coupled system ordinary differential equations (ODEs)
where L _{1}, L _{2} ∈ L and w = w(u, v).
Exponential integrators separate the linear term involving L, which is solved exactly by a matrix exponential, from the nonlinear term. The theory of numerical methods for the time integration of semilinear problems has been proposed by the application of the exponential methods. Cox and Matthews [6] presented derivation of exponential time differencing (ETD) methods. Few years later, a modification of the ETD RungeKutta methods of Cox and Matthews was made by Kassam and Trefethen [21], and it is from their paper that we present some details of the scheme. A new algorithm for the implementation of the exponential methods has been discussed in [9], that the algorithm evaluates the operator by the exponential methods with a quadrature formula that converges. Hochbruck and Ostermann [15] discussed further on the class of explicit multistep exponential and explicit exponential RungeKutta methods. In this paper, we use both the fourthorder exponential time differencing RungeKutta method [6, 21] and the fourthorder exponential multistep method of Adamstype [6, 16].
Exponential time differencing method
We present a brief introduction to the derivation of exponential time differencing RungeKutta and multistep methods of Adamstype along with their stability regions. Details of their derivations can be found in [6, 16] and the references therein. The exponential time differencing idea, applied here for the u component, involves the use of the integrating factor, e ^{− L t}. We multiply equation (14) by this factor and then integrate it over a timestep to obtain
This equation is known to be exact [6, 21]. Various ETD schemes come from how one approximates the integral on the right hand side in the above equation. The article by Cox and Matthews [6] contains various approximations to the integral. They first presented the sequence of recurrence formulae that give higherorder approximations of a multistep method of Adamstype. In their work, a general ETD scheme of orders was proposed as
The coefficients g_{ j } are obtained by the recurrence relation
By setting s = 4 in the explicit integrating formula (16), we obtain the fourthorder ETD scheme of Adamstype
where
denoted in this paper as ETDADAMS4.
Similarly, Cox and Matthews derived a set of ETD schemes that are based on RungeKutta timestepping, which they call ETDRK schemes. We only use the fourthorder scheme which we denoted as ETDRK4 in this paper. On setting s = 4 again in (16), we have the ETDRK4 formula
where
To circumvent the pronounced vulnerability of error cancelations in the higherorder ETDADAMS4 and ETDRK4 schemes [21], and to generalize the ETD schemes to nondiagonal problems, modified schemes are proposed with the aid of complex contour integral
to evaluate the coefficients of these schemes. Further details on derivations and applications of ETD Adamstype and ETD RungeKutta methods can be found in [6, 21, 41].
Stability analysis
We investigate the stability of the fourthorder ETDADAMS4 (18) and ETDRK4 (19) schemes by linearizing the nonlinear autonomous system [8, 19]
with N(w(t)) the nonlinear part. We suppose that there exists a fixed point w _{ 0 } such that Lw _{ 0 } + N(w _{ 0 } ) = 0. Linearizing about this fixed point, we obtain
where w(t) is now the perturbation of u _{ 0 } and λ = N′(w _{0}) is a diagonal or a block diagonal matrix containing the eigenvalue of N. In an attempt to keep the fixed point u0 stable, we require that Re(L + λ) < 0, for all λ. It is naturally important for a numerical method to satisfy this property so as to cover as much dynamics as possible.
When applying ETDADAMS4 (18) to the linearized problem (22), a polynomial equation of the orderfour in r is obtained in the form
where
In the real (x, y) plane, the righthand boundary for ETDADAMS4 scheme corresponds to substituting r = 1 in equations (23) is the line x + y = 0. The corresponding lefthand boundary for substituting r = −1, also in (23), is given by the curve
as displayed in Fig. 1.
In a similar fashion, the application of ETDRK4 method (19) to the linearized problem (22) leads to a recurrence relation
where
where x = λh, y = Lh. We can define the amplification factor for ETDRK4, r(x, y) for y > 0. If y = 0, the amplification factor becomes 1 − x + x ^{2}/2 − x ^{3}/6 + x ^{4}/24. Hence, we can see that the stability curve of ETDRK4 at y = 0 coincides with that of the classical fourthorder RungeKutta method, Fig. 2(a). We also see that lim_{ x,y → 0}∂_{ x } r(x, y) = − 1 and lim_{ x,y → 0}∂_{ y } r(x, y) = − 1. Hence, the absolute value of the amplification factor is given as r(x, y) ≤ 1.
The boundary of the stability region can be determined by setting r = e ^{iθ}, for θ ∈ [0, 2π]. We plot the stability region in the complex x plane and displayed in Fig. 2, where the horizontal and vertical axes represent the real and imaginary of x, respectively.
Numerical examples and results
In this section, numerical methods we discussed above are now applied to the three major classes of the LotkaVolterra twospecies models. In addition, comparison with other adaptive methods are made to justify the effectiveness and accuracy of the present method. A possible extension to two space dimensions is considered, since it is in higher dimensions that most of the ideas reported are of serious value.
Predatorprey system
It is clear from our introduction that predatorprey models are similar in description to both parasite and parasitoid models. A typical example of predatorprey model [11, 32] is the reactiondiffusion system
where U and V are the densities of the prey and predator respectively, D _{1} > 0 and D _{2} > 0 are diffusion coefficients for the prey and predator. α, β, γ, δ, h and K are positive parameters. The term αU(1 − U/K) represents the logistic growth, α is the intrinsic growth rate, and K the carrying capacity. The term γV is the percapita prey reduction due to consumption by the predator, and β describes the intensity of predation.
To reduce the number of parameters in (26), we nondimensionalize the model by rescaling the variables as
to yield
For the linear stability, we have to analyze the stability criteria of the nondiffusive system [17, 31, 42]. The spatial model (28) has the corresponding nondiffusive systems
with just three parameters μ > 0, ψ > 0 and φ > 0. There are other choices for the change of variables to put the system in dimensionless form, but we opt for the choice that suits our purpose since the dimensionless groupings used here give relative measures of the effect of dimensional parameters. For instance, ' now becomes the ratio of the linear growth rate of the predator to that of the prey, for ψ < 1. We expect the prey to reproduce faster than the predator otherwise the system will go into extinction.
At equilibrium, \( f\left(\widehat{u},\widehat{v}\right)=g\left(\widehat{u},\widehat{v}\right)=0 \), since the steady state populations û and \( \widehat{v} \) are solutions of du/dt = dv/dt = 0. Hence,
Naturally, for the dynamical system under consideration to be biologically meaningful, we should have both u ≥ 0, v ≥ 0 at all times. We observe from (30) that the system (28) has three positive steady states \( \left(\widehat{u},\widehat{v}\right) \), the two trivial states or saddle points are at point (0, 0) which describes complete extinction of both prey and predator and point (1, 0), which shows that the predator is absent leading to unbounded logistic growth of the prey species. The stationary point \( \left(\widehat{u},\widehat{v}\right) \) corresponding to the existence of predator and prey, bearing in mind that for the system under consideration to be biologically meaningful, the parameters must be strictly restricted to the positive quadrants, gives
The stability of the steady or equilibrium states are the singular points in the phase plane of (28). To determine them, we let
where A is regarded as the community matrix with eigenvalues given by
For stability, we require that Reλ < 0. Hence, the necessary and sufficient conditions for linear stability become
On substituting û in equation (31) provides the stability conditions in terms of the positive parameters μ, ψ and φ.
Results in onedimension for system (26)
Numerical results of the predatorprey system are shown in onedimension. The initial data and parameter values are given in the figure caption. The initial data are chosen as a result of small perturbations of the steady state solutions û and \( \widehat{v} \) of the spatially homogeneous system. By varying the choice of parameters lead to different spatial patterns, such as oscillatory smooth, intermittent structure and spatiotemporal patterns. It should be noted that other onedimensional spatial structures that are not captured here are possible, depending on the choice of the parameter values and initial data.
Figures 3 and 4 represent the unrealistic and realistic population dynamics of the predatorprey systems. The system with nonlinear part as described in Garvie [11] is quite unrealistic due to the choices of parameters used in transforming the system into a dimensionless form. This shortcoming actually motivates us to choose some appropriate parameters since it is always helpful to write the system in nondimensional form. Nondimensionalisation plays an important role when carefully considered because it reduces the number of parameters by grouping them in a more meaningful manner. So, the system described in Fig. 3 is totally unrealistic as it is prone to danger of extinction of the prey species that would in turn results to total breakdown of the ecosystem since all the predators will die out in absence of food. In Fig. 4, spatiotemporal oscillations arise and population oscillations are transient and regular. It should be noted that due to the formation of spatial pattern, the two species can dynamically coexist.
Results in twodimension for system (26)
We intend to mimic the twodimensional results obtained for the predatorprey system in [11, 27], we experiment with the same initial data
so as to induce a nontrivial spatiotemporal dynamics of the homogeneous stationary states û and \( \widehat{v} \). In Fig. 5, numerical simulations was done on a square domain size [0, 700] × [0, 700], with parameter values D = 0.1, μ = 0.2, ψ = 2, φ = 0.5 at nontrivial state \( \left(\widehat{u},\widehat{v}\right) \) =(6/35, 116/245). As simulation time is increased from t = 200 to t = 500, the spiral patterns in (a, b) are disjointed and spreads out in the domain to form a stripelike structures with emergence of some spots underneath. It should be mentioned that if the simulation time is further increased, say to t = 1500 and above, there is every tendency of getting a Turing and more complicated spatiotemporal patterns. In addition, we realized that the choice of initial conditions can influence the type of spatiotemporal dynamics of a reactiondiffusion problem in ecosystems.
A close look at the first and second columns in Fig. 5 have revealed that both predator and prey species have a similar distribution. As a result, our pattern formation analysis is henceforth restricted to only one distribution. We also observe in our experiments that increase in domain size actually results to increase in computational time. Henceforth, we choose to simulate with a smaller square domain of size [0, 250] × [0, 250].
Competitive system
Competition model describes a situation in which two or more species compete for the same (sufficient or insufficient) resources like food, territory or in some way inhibit each other of growth. For simplicity, and by following the approach we used for the predatorprey model, we consider here the twospecies LotkaVolterra competition model
with species U and V having logistic growth in the absence of the other. The parameters α _{1} and α _{2} represent their linear birth rates, β _{1} and β _{2} measure the competitive effect of V on U and vice versa, δ _{1} and δ _{2} stand for the diffusion coefficients of species U and V, and K _{1} and K _{2} are their respective carrying capacities.
Again, we nondimensionalize (36) by introducing a set of carefully selected dimensionless variables
As suggested by Medvinsky et al. [27] and Garvie [11], the local stability analysis will always grant a deeper understanding and will provide important information on the choice of parameters for numerical integration. Like the previous case, we continue with the local stability analysis in the absence of diffusion. Using (37) in (36), we obtain
For the linear stability analysis, we consider the case of spatially homogeneous solutions, in which the spatial model (38) is equivalent to the system of ordinary differential equations
Here, we regard the steady states and phase plane singularities, û and \( \widehat{v} \) as the solutions of f(u,v) = g(u, v) = 0. This gives four positive equilibrium states,
The good thing is that, all the four steady states exist in the positive quadrant which make the whole process meaningful in the biological and ecological contexts.
The first three stationary states are trivial whereas the last one is nontrivial. The state (0,0) corresponds to total washout state of the two species, the second state (1,0) stands for the existence and extinction of species u and v respectively and the third trivial state (0; 1) indicate that only species v will exist. It is obvious that none of the three trivial states could give a meaningful interpretation about the competition model, therefore, there is the need to explore further the nontrivial equilibrium state \( \left(\widehat{u},\widehat{v}\right) \). The points (0,0), (1,0) and (0,1) are all unstable (0,0) is an unstable node, (1,0) and (0,1) are saddle point equilibria. From (39), for f = g = 0, we have that (u − u ^{2} − φuv) = 0, it follows that either u = 0 or 1 − u − − φv = 0 and also from the second equation, μ(v − v ^{2} − φuv) = 0 which implies, μv = 0 and 1 − v − − φu = 0.
Now the Jacobian or community matrix for this system evaluated at \( \left(\widehat{u},\widehat{v}\right) \) is
The point (0, 0), is unstable since the eigenvalues λ obtained from
are λ _{1,2} = (1, μ). At the point (1, 0), the community matrix A gives
Hence, λ _{1,2} = (−1, μ(1 − ψ)). Therefore, the steady state \( \left(\widehat{u},\widehat{v}\right)=\left(1,0\right) \) is stable if ψ > 1 and unstable otherwise. In the same manner, we can see that the steady state (0, 1) has the community matrix A satisfying
The corresponding eigenvalues are λ _{1,2} = (−μ, (1 − φ)). This means that the steady state \( \left(\widehat{u},\widehat{v}\right) \) = (0, 1) is stable if φ > 1 and unstable if φ < 1.
For the fourth steady states, we have matrix,
The eigenvalues in this case are
Clearly, the stability of the steady state depends on the size of the positive parameters μ, φ and ψ subject to various cases such as; (φ > 1, ψ > 1), (φ > 1, ψ < 1), (φ < 1, ψ > 1) or (φ < 1, ψ < 1). A biological interpretation of Fig. 6(b) suggests that because the carrying capacity of species u is so high, this species is not limited by the resources to the extent at which species v seems to be. Stable coexistence occurs when the isoclines are arranged as in Fig. 6 (a) for K _{1} < K _{2}/ψ and K _{2} < K _{1}/φ. The populations converge on the intersection of the isoclines regardless of the initial population densities. The intersection point of the two lines gives the positive steady state as in (a) where the point (1.4, 1.4) corresponds to (1/φ, 1/ψ). The locations of the isoclines in (b) dictate that species u outcompetes species v, the point (1/φ, 1/ψ) corresponds to the value (1.6, 0.6) of species u and v, respectively. Clearly, on rearranging, we can see that ψ < K _{2}/K _{1} and φ < K _{1}/K _{2}, and these competition coefficients must be made as small as possible relative to the ratio of its carrying capacity to that of other species. These conditions must hold for both species simultaneously, and this is possible only if the carrying capacities of the two species are similar in such a way that their ratio is close to one. Figure 7 (a) describes the species declining population density associated with the competitive system (38), panels (b,c) refer to the time series solution, and (d) corresponding to the species phase plane diagram.
Two dimensional results for model (38)
We also carry out a twodimensional numerical simulations of the spatially extended competitive model (38). We employed the initial conditions (35) and the zeroflux boundary conditions on a square domain size of [0, 250] × [0, 250] with timestep Δt = 0.005 and grid width Δh = 0.25. Here the parameter values are set as
In Fig. 8, we show three typical Turing patterns obtained at (a) \( \left(\widehat{u},\widehat{v}\right) \) =(11/45, 110/253) for t = 300 and (b) \( \left(\widehat{u},\widehat{v}\right) \) =(0.05, 0.062) for t = 500. In both panels, we noticed the formation of Turing spots pattern emanating from the center of the domain, as a result, we fixed the parameter values as in (b) and increase the simulation time to t = 700. A pattern containing the mixture of spots and moonlike stripe patterns emerged in (c). From (ac) one can observed that irregular patterns prevail in the entire domain. However, the three patterns are essentially different from one another, because of their different wavelengths. We believe the possibility of getting other Turing dynamical structures depending on the choice of initial data and the length of simulation.
Mutualism system
This is a type of association in theoretical ecology in which the existence of one species has no negative influence on the other. This type of model receives little attention and has not been studied as others, even though its importance is comparable to that of preypredator and competition models. To start with, we shall analyze briefly the twospecies model
where F(U, V) = α _{1} U(1 − U/K _{1} + β _{1} V/K _{1}) and G(U, V) = α _{2} U(1 − V/K _{2} + β _{2} U/K _{2}) are the nonlinear reaction terms for the two species U and V, respectively. And σ _{1}, σ _{2}, α _{1}, α _{2}, β _{1}, β _{2}, K _{1}, K _{2} are all positive parameters. This system looks similar to equation (36), with exception that β ' s are treated positive in this case. We then nondimensionalize using the parameters
which on substitution in (43) yields
Again, by following the linear stability analysis, we study the stability criteria for the nondiffusive system
It is not difficult to see that the steady states \( \left(\widehat{u},\widehat{v}\right) \) for this system are
The Jacobian or community matrix for this system is
Proceeding in a similar manner like those for the previous cases, we can easily show that the points (0, 0), (1, 0) and (0, 1) are all unstable; the point (0, 0) is unstable node while (1, 0) and (0, 1) are the saddle point equilibria, whereas the fourth steady state for 1 − φψ > 0 (located in the positive quadrant) is a stable equilibrium. Mutual display of the species is reflected in Fig. 9, panel (a) shows linear behaviour of species u and v. Each of the species experienced an unbounded population growth since the existence of one has no effect on the other and their relationship is linear as in (b).
Two dimensional results for model (45)
Following [34], we take the boundary conditions
subject to the axisymmetric initial conditions
where ς ^{2} = x ^{2}/2 + y ^{2}. We perform some numerical simulations of the dynamical model (45) on the domain size [0, 250] × [0, 250] with timestep Δt = 0.05 and grid width Δh = 0.5, \( \left(\widehat{u},\widehat{v}\right) \) =(0.06125, 0.25). We fixed other parameters as in (42) to obtain Fig. 10. In the simulations at t = 500, the pattern structures start appearing like a cluster of stripes right from the domain center. It spreads out into irregular stripes as simulation time increased to t = 1000. Later, with further increase in time, the long stripes break into spots at t = 1500 as in (c). In panel (d) at t = 2000, spot patterns have covered the entire domain. Pure Turing spots pattern is achievable if the simulation time is further increased.
In order to justify the suitability and accuracy of the ETDADAMS4 and ETDRK4 schemes, we carried out numerical experiments on the three dynamical systems considered in this paper that is, the prey–predator system (28), competitive system (39), and the mutualism or symbiosis system (45). The performance of ETDRK4 and ETDADAMS4 are investigated and compared with the family of exponential time differencing multisteps schemes of order four, five and six which we denoted in this paper for brevity as ETDM4, ETDM5 and ETDM6 respectively.
It would go beyond the scope of this paper to give a complete classification of exponential integrators used for comparison. We focus on exponential time differencing method of Adamstype and exponential time differencing RungeKutta method, and we have mentioned earlier how they can be treated in the common framework of explicit exponential integrators. Details of these schemes are well classified in [16, 41] and references therein, with historical survey offered by Minchev and Wright [30].
We report the maximum relative errors of the solution defined by
where û _{j} is a goldstandard run computed with the schemes at Δt = 1/2048 and u _{j} is computed values of the solution u at point j, and N is the number of interior points defined on the collocation interval
Figure 11 (a) shows the performance of the schemes when applied to the preypredator system (28) at parameter values t = 1, μ = 0.1, ψ = 0.08, φ = 0.01, δ = 0.01 for N = 200. Panel (b) is obtained with parameters t = 1, μ = 0.5, ψ = 0.15, φ = 0.15, δ = 0.5 and N = 200 for the competitive system (39). The performance of the schemes when applied to the mutualism system (45) at parameter values t = 1, μ = 0.5, ψ = 0.5, φ = 0.5, δ = 0.1, N = 200 is shown in panel (c). We compute the relative errors using a goldstandard run obtained with the schemes using Δt = 1/2048 and compare with various time steps 1/2^{ρ}, ρ = 1, …, 10 [34, 36].
It is obvious from the results presented in Fig. 11 that the ETDRK4 has a better convergence when compared to other exponential time differencing methods for each of the problems considered in this paper. Due to the similarity and the choices of parameters used in the simulations of the competitive and the mutualism systems, one observes that the schemes have similar behaviour. The difference is noticeable in their amplitudes. The ETDADAMS4 competes very well with ETDRK4 when applied to the dynamical systems but the ETDRK4 appears to have the overall credit.
The following experiment in Table 1 was performed in onedimension with predatorprey system (28) in a smaller domain size (0; 100) and the computation was terminated at final time t = 1,…, 4. The parameter values are: μ = 0.4, ψ = 0.08, φ = 0.05, Δt = 0.25 for N = 200. We use the builtin Matlab tic  toc to check the computational time of the schemes. Both schemes runs in seconds. Our numerical experiments in onedimension demonstrate a strong case for abandoning the ETDM4, ETDM5 and ETDM6 schemes. In obtaining the 2D results in Fig. 5, it was observed that the ETDRK4 timestepping scheme performed about two times faster than the ETDADAMS4 scheme. That is, the computational time required for ETDADAMS4 is about 48 % more than that of the ETDRK4. As a result, we carried out the 2D experiments with the ETDRK4 scheme.
Conclusions
In this paper, firstly, the dynamic complexities of the ecological models consisting of preypredator, competitive and mutualism reactiondiffusion dynamics are studied by considering their linear stability analysis in the absence of diffusion, and secondly by the numerical approach with the presence of diffusion. We discretized the governing models in space using a fourthorder central finite difference scheme and integrate the resulting ODEs with the exponential time differencing schemes whose formulations were based on the RungeKutta and multistep methods of Adamstype. We investigate the stability of the schemes and plots their stability regions. We present the results in both one and two dimensions to unveil their pattern formations. The numerical experiments in 2D reveal some of the typical patterns such as stripes and spots, as well as irregular snakelike patterns. Further, we compared the results obtained with both ETDADAMS4 and ETDRK4 for each of the dynamics, with their exponential fourth, fifth and sixthorders counterparts denoted as ETDM4, ETDM5 and ETDM6, respectively, and observed that the ETDRK4 is most reliable and computationally promising in terms of efficiency and accuracy when compared to other methods used in this paper. It worth mentioning that the methodology presented in this work can be extended to higher dimensional practical problems.
Abbreviations
 ETDRK4:

Fourthorder exponential time differencing RungeKutta method
 ETDM:

Exponential time differencing multistep method
 ETDADAMS:

Exponential time differencing method of Adamstype
References
 1.
Allee WC. The Social Life of Animals. New York: Norton; 1938.
 2.
Allen LJS. An Introduction to Mathematical Biology. New Jersey: Pearson Education, Inc.; 2007.
 3.
Amarasekare P. Interactions between local dynamics and dispersal: insights from single species models. Theor Popul Biol. 1998;53:44–59.
 4.
Baek H, Jung DI, Wang Z. Pattern formation in a semiratiodependent predatorprey system with diffusion. Discr Dyn Natur Soc. 2013;2013(657286):14. doi:10.1155/2013/657286.
 5.
Berryman AA. Population Systems: A General Introduction. New York: Plenum Press; 1981.
 6.
Cox SM, Matthews PC. Exponential time differencing for stiff systems. J Comput Phys. 2002;176:430–55.
 7.
Dennis B. Allee effects: population growth, critical density, and the chance of extinction. Nat Res Model. 1989;3:481–538.
 8.
Du Q, Zhu W. Analysis and applications of the exponential time differencing schemes and their contour integration modifications. BIT Numer Math. 2005;45:307–28.
 9.
LopezFernandez M, Palencia C. On the numerical inversion of the Laplace transform of certain holomorphic mappings. Appl Numer Math. 2004;51:289–303.
 10.
Gakkhar S, Naji RK. Order and chaos in s food web consisting of a predator and two independent preys. Commun Nonl Sci Numer Simul. 2005;10:105–20.
 11.
Garvie M. Finitedifference schemes for reactiondiffusion equations modeling predatorpray interactions in MATLAB. Bullet Math Biol. 2007;69:931–56.
 12.
Garvie M, Trenchea C. Spatiotemporal dynamics of two generic predatorprey models. J Biol Dyn. 2010;4:559–70.
 13.
Gyllenberg M, Hemminki J, Tammaru T. Allee effects can both conserve and create spatial heterogeneity inpopulation densities. Theor Popul Biol. 1999;56:231–42.
 14.
Harmon JP, Andow DA. Indirect effects between shared prey, predictions for biological control. Biol Control. 2004;49:605–25.
 15.
Hochbruck M, Ostermann A. Exponential RungeKutta methods for parabolic problems. Appl Numer Math. 2005;53:323–39.
 16.
Hochbruck M, Ostermann A. Exponential multistep methods of Adamstype. BIT Numer Math. 2011;51:889–908.
 17.
Holmes EE, Lewis MA, Banks JE, Veit RR. Partial differential equations: Spatial interactions andpopulation dynamics. Ecology. 1994;75:17–29.
 18.
Holt RD. Predation, apparent competition, and the structure of prey communities. Theor Popul Biol. 1977;12:197–229.
 19.
de la Hoz F, Vadilo F. An exponential time differencing method for the nonlinear schrodinger equation. Comput Phys Commun. 2008;179:449–56.
 20.
Janzen DH. The natural history mutualisms. In: Boucher DH, editor. Biol Mutual. Oxford: Oxford University Press; 1985. p. 44–99.
 21.
Kassam AK, Trefethen LN. Fourthorder timestepping for stiff PDEs. SIAM J Sci Comput. 2005;26:1214–33.
 22.
Kawasaki K, Mochizuki A, Matsushita M, Umeda T, Shigesada N. Modeling spatiotemporal patternsgenerated by Bacillus subtilis. J Theor Biol. 1997;188:177–85.
 23.
Kot M. Elements of Mathematical Ecology. United Kingdom: Cambridge University Press; 2001.
 24.
Lotka AJ. The Elements of Physical Biology. Baltimore: Williams and Wilkins; 1925.
 25.
Lotka AJ. The growth of mixed populations, two species competing for a common food supply. J Washington Acad Sci. 1932;22:461–9.
 26.
Malchow H. Spatiotemporal pattern formation in nonlinear nonequilibrium plankton dynamics. Proc Roy Soc London B. 1993;251:103–9.
 27.
Medvinsky AB, Petrovskii SV, Tikhonova IA, Malchow H, Li BL. Spatiotemporal complexity of plankton and fish dynamics. SIAM Rev. 2002;44:311–70.
 28.
Mendez V, Fedotov S, Horsthemke W. ReactionTransport Systems: Mesoscopic Foundations, Fronts, and Spatial Instabilities. Berlin Heidelberg: Springer; 2010.
 29.
Mimura M, Sakaguchi H, Matsushita M. Reactiondiffusion modelling of bacterial colony patterns. Physica A. 2000;282:283–303.
 30.
Minchev BV, Wright WM. A review of exponential integrators for first order semilinear problems, Technical Report NTNU. Department of Mathematical Sciences, Norwegian University of Science and Technology, (2005), Preprint.
 31.
Murray JD. Mathematical Biology. Berlin: Springer; 1989.
 32.
Murray JD. Mathematical Biology I: An Introduction. New York: Springer; 2002.
 33.
Murray JD. Mathematical Biology II: Spatial Models and Biomedical Applications. Berlin: Springer; 2003.
 34.
Owolabi KM. Robust IMEX schemes for solving twodimensional reactiondiffusion models. Int J Nonlinear Sci Numer Simul. 2015;16:271–84.
 35.
Owolabi KM, Patidar KC. Numerical solution of singular patterns in onedimensional GrayScottlike models. Int J Nonlinear Sci Numer Simul. 2014;15:437–62.
 36.
Owolabi KM, Patidar KC. Higherorder timestepping methods for timedependent reactiondiffusion equations arising in biology. Appl Math Comput. 2014;240:30–50.
 37.
Petrovskii S, Kawasaki K, Takasu F, Shigesada N. Diffusive waves, dynamic stabilization and spatiotemporal chaos in a community of three competitive species. Japan J Industr Appl Math. 2001;18:459–81.
 38.
Petrovskii S, Malchow H. Wave of chaos: new mechanism of pattern formation in spatiotemporal population dynamics. Theor Popul Biol. 2001;59:157–74.
 39.
Petrovskii S, Morozov AY, Venturino E. Allee e_ect makes possible patchy invasion in a predatorprey system. Ecol Lett. 2002;5:345–52.
 40.
Satnoianu RA, Menzinger M, Maini PK. Turing istabilities in general systems. J Math Biol. 2000;41:493–512.
 41.
Schmelzer T, Trefethen LN. Evaluating matrix functions for exponential integrators via CaratheodoryFejer approximation and contour integrals. Elect Trans Numer Anal. 2007;29:1–18.
 42.
Sherratt J. Periodic travelling waves in cyclic predatorprey systems. Ecol Lett. 2001;4:30–7.
 43.
Volpert V, Petrovskii S. Reactiondiffusion waves in biology. Phys Life Rev. 2009;6:267–310.
 44.
Volterra V. Fluctuation in abundance of the species considered mathematically. Nature. 1926;118:558–60.
 45.
Volterra V. Variations and Flunctuations of the Numbers of Individuals in Animal and Species Living together, Reprinted in 1931 in R.N. Chapman, Animal Ecology. New York: McGrawHill, 1926.
 46.
Wang W, Liu QX, Jin Z. Spatiotemporal complexity of a ratiodependent predatorprey system. Phys Rev E. 2007;75:1539–3755.
 47.
Wang W, Zhang L, Wang H, Li Z. Pattern formation of a predatorprey system with Ivlevtype function response. Ecol Model. 2010;221:131–40.
 48.
Yu H, Zhong S, Agarwal RP. Mathematics and dynamic analysis of an apparent competition community model with impulsive effect. Math Comput Model. 2010;52:25–36.
 49.
Yu H, Zhong S, Agarwal RP, Xiong L. Species permanence and dynamical behavior analysis of an impulsively controlled ecological system with distributed time delay. Comput Math Applic. 2010;59:3824–35.
Acknowledgements
OK acknowledges the partial financial support from the Federal Government of Nigeria under the 2009 Education Trust Fund Academic Staff Training and Development (AST&D) Intervention. The research contained in this paper is also supported by the South African National Research Foundation (NRF).
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
Both OK and PK conceived the research work. Manuscript writing and simulation: OK. Model analysis: OK. Both authors proofread and approved the final manuscript.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Coexistence
 Competitive
 Exponential time differencing
 Predatorprey
 Mutualism
 Multistep methods
 Reactiondiffusion systems
 Stability analysis
AMS Subject Classification
 Primary 70K05
 65M20
 Secondary 70K20
 74H65