Abstract: Mathematical modelling plays an important role in the study of prey-predator interactions. The prey-predator interactions can be modelled by using the basic Lotka-Volterra equations, which describe the dynamics of biological systems in which two species i.e., the prey and predator interact. However, the presence of diffusion mechanism, changes the behaviour and nature of the whole model. It is now a reaction-diffusion system, which takes the form of partial differential equations and very difficult to be solved analytically. Thus, most researchers have turned to numerical simulations so as to study the behaviour of the system. In this study, the numerical simulations technique, namely finite-difference method, is employed so as to study the rich spatio-temporal dynamical structure of diffusive prey-predator model. The formations of numerous patterns are observed i.e., spiral waves, patchy structures, spiral defect chaos and spatio-temporal chaos due to Turing instability in our prey-predator model. These two-dimensional patterns are in fact very beautiful and very interesting to be observed and interpreted accordingly from the ecological point of views.
INTRODUCTION
Mathematical modelling plays an important role in the study of population dynamics and ecological interactions. Ecological systems are characterized based on the interactions between multiple species and the natural environmental factor. There are several types of interaction between species, namely, competition, mutualism and prey-predator. The prey-predator interaction has received great attention since the early days of ecological science discipline. The basic prey-predator interaction is normally described by using a system of ordinary differential equations, which modelled the spatial distribution of species as the time evolves.
The presence of diffusion mechanism however changes the behaviour and nature of the whole model. It is now a partial differential equation, which can be categorized as a reaction-diffusion system. The diffusive prey-predator model has been studied extensively by Murray (1989), Levin et al. (1993), Malchow (1993), Medvinsky et al. (2002) and Holmes et al. (1994).
The inclusion of diffusion terms however has made the prey-predator model tend to be complicated and it becomes very difficult to analyze and solve analytically. According to Wang et al. (2007) the non-uniform stationary state of prey-predator model, which corresponds to the spatio-temporal patterns, cannot be found analytically. Thus, many scientists resort to numerical simulations so as to study the behaviour of the system. If the model parameters are selected properly, the numerical simulations will give rise to rich spatio-temporal dynamical structure and amazingly beautiful two-dimensional patterns, viz., homogeneous distributions, stationary patterns (stripe-like or spotted or coexistence of both), emergence of spatio-temporal chaos and spiral waves. Refer to Wang et al. (2007) for more discussion on this matter.
Apart from this, many researchers have performed analytical and numerical investigations of diffusive prey-predator model from different aspects. They have studied the problem of biological invasion in prey-predator model and the dynamics of population through non-Turing pattern formation with specified initial conditions. Refer to Sherratt (2001), Shigesada and Kawasaki (1997) and Petrovskii and Malchow (2001) for the discussion on this matter.
In general, this study was concerned with another situation. The main objective was to investigate the spatio-temporal pattern of diffusive prey-predator model and the emergence of irregular chaotic pattern as a result of prey-predator interaction. Plus, this study also intended to explore the occurrence of diffusion-induced instability (Turing instability) and its effect to the dynamics of prey-predator model.
MATHEMATICAL MODEL
A necessary condition to observe Turing instability is that predator must diffuse
faster than prey i.e. α<<β, whereby α and β are the
diffusion constants for our prey-predator model, respectively. Theoretically,
defining
Now, consider a system of reaction-diffusion equation which takes the form:
| (1) |
where, u and v represent the prey and predator population densities, whilst α and β are diffusion constants of both prey and predator, respectively. For two dimensional case:
The reaction-diffusion Eq. 1 is subjected to the following initial condition:
| (2) |
and Neumann boundary condition (i.e., zero-flux boundary condition):
| (3) |
where,
In order to study the spatio-temporal dynamics of diffusive prey-predator model, let presume that the local growth of the prey species is logistic and the predator community displays the Holling type II functional response. So, F (u,v) and G (u,v) are specified as below:
| (4) |
where, φ is the maximal growth rate of prey population, γ is the carrying capacity for the prey population and H is the half-saturation density for the prey population.
Introducing the dimensionless variables and parameters:
| (5) |
and this will result in the following reaction-terms:
| (6) |
So, the diffusive prey-predator model becomes:
| (7) |
By using model (7), the local stability analysis is conducted so as to determine the equilibrium states of prey-predator model and then the nature of these equilibrium points are investigated for further analysis.
LOCAL STABILITY ANALYSIS
The local stability analysis is performed so as to see the local system dynamics (without diffusion i.e., α = β = 0) and later the spatio-temporal pattern formation is analyzed. According to Medvinsky et al. (2002) and Garvie (2007), it is important to consider the local dynamics of the system as this will provide guidelines on the suitable choice of parameters for numerical simulations.
Firstly, the steady states of the prey-predator model are determined by setting F (u*, v*) = G (u*, v*) = 0 in the reaction terms:
| (8) |
From Eq. 8, there are three positive equilibrium states (u*, v*), two of the states are trivial and one is non-trivial. Two of the trivial states are (0, 0) which corresponds to total extinction of prey-predator population and (1, 0) which indicates the persistence of prey and extinction of predator population. These two trivial equilibrium states are in fact saddles, which are standard results and can be found in many texts on differential equations.
Thus, let us focus the analysis on the non-trivial equilibrium state (u*, v*). Further calculation shows that another equilibrium state of Eq. 8 is:
| (9) |
Equation 9 is a non-trivial equilibrium state which indicates the co-existence case of prey and predator population. This steady state is very interesting to be analyzed in detail. Linearizing the system of Eq. 8 about the equilibrium Eq. 9 will result in the Jacobian matrix, A:
| (10) |
with each element of A are calculated as follows:
| |
Each element of Jacobian matrix A has been determined for the stability analysis. Now, consider the system of equations below:
| (11) |
which is actually the solutions to Eq. 8. C1, C2, D1 and D2 are constants which corresponds to the eigenvectors of A, i.e.:
| (12) |
where, λ1 and λ2 are eigenvalues of A. The eigenvalues of A are the roots of characteristic equation where:
| (13) |
The characteristic Eq. 12 has two roots, which are:
| (14) |
Based on the values of λ1 and λ2, the equilibrium points can be classified as node, saddle point, spiral or centre. It will also tell us about the stability, whether the equilibrium point is stable or unstable.
Next, the local stability of Eq. 8 is analyzed by choosing δ = 0.4, ε = 2 and μ = 0.8 as an example. So the non-trivial equilibrium point (based on formula (Eq. 9)) is:
and the Jacobian matrix A (based on formula (Eq. 10)) evaluated about the equilibrium is:
The resulting eigenvalues are:
Since these eigenvalues are complex conjugate pair with positive real part, thus it will result in the formation of unstable spiral. This result will be verified later on, when the numerical simulations in two-space dimensions is conducted, so as to study the spatio-temporal dynamics of diffusive prey-predator model (Fig. 1).
| Fig. 1: | Local phase plane of model (8) with parameters, δ = 0.4, ε = 2 and μ = 0.8 |
NUMERICAL SIMULATIONS
Now, let study the rich dynamical structure as a result of prey-predator interactions through numerical simulations. For the sake of brevity, choose α = 1; thus model Eq. 7 becomes:
| (15) |
with:
and 0 < x < Lx, 0 < x < Ly. Lx and Ly are the length and width of the spatial domain. The diffusive prey-predator model (15) is solved numerically using finite-difference method with zero flux boundary condition. And the choice of initial condition is as follow:
| (16) |
is chosen specifically in this form so as to induce a non-trivial spatio-temporal dynamics. This “perturbed” initial condition will also be used in order to see the sensitivity behaviour to initial condition, since it is one of the main characteristics for a system to be chaotic.
The simulation studies are started by using the same parameters as in Local Stability Analysis section i.e., δ = 0.4, ε = 2 and μ = 0.8 and this result in the non-trivial equilibrium point of:
This equilibrium point is actually a spiral, since the eigenvalues is a pair of complex conjugate with positive real part.
The model (15) is solved in two dimensions with a square numerical domain LxxLy where Lx = Ly = 700; using parameters δ = 0.4, ε = 2 and μ = 0.8. The initial condition is chosen as in (16) and β = 1 (so as to exclude the possibility of Turing instability). Refer to Baurmann et al. (2007), Camara and Aziz-Alaoui (2008) and Wang et al. (2010) for the discussion on this matter. Now, consider the two-dimensional approximate prey-predator densities below that are obtained from numerical simulations by increasing the time of experiment.
The evolution of regular spiral pattern as time increases is very interesting to be discussed. Based on Fig. 2a, observe the formation of spiral pattern, which confirm the local stability analysis result that has been obtained. After the spirals form Fig. 2a, these structures grow slightly for a certain time and later their spatial structure becoming more distinct (Fig. 2b). Then, it continues to evolve and as T = 700, the spirals become separated into two circles. Notice that, the destruction of the spirals begins in their centres (Fig. 2c). These circles steadily grow bigger and as T = 1000, the emergence of patchy structure inside these circles is observed (Fig. 2d).
To investigate this evolution further, the time of experiments T is increased in the range of (2000, 4000). Consider the snapshots of spatio-temporal pattern below:
Observe that the irregular patchy structure steadily grow as time increases, refer to Fig. 3a-b. Later, the irregular patches have covered the whole domain as T = 4000 (Fig. 3c) and this is actually consistent with the results reported in Garvie (2007) and Medvinsky et al. (2002). They studied the spatio-temporal dynamics in aquatic community and discovered that the size of these patches has been related to the characteristics length of observed plankton patterns in the ocean. So, the reaction diffusion Eq. 15 exhibit spatio-temporal chaos causing plankton patchiness in aquatic system.
In Fig. 2 and 3, the evolution of spiral patterns with respect to time has been studied. Now for further investigation, the effect of varying parameter δ to the spatio-temporal dynamics of model (15) is studied. By using the same parameters as in Fig. 2 and 3, model (15) is solved in two dimensions with a square numerical domain LxxLy where Lx = Ly = 500. Consider the two-dimensional approximate prey-predator densities below that has been obtained from numerical simulations by varying the values of δ:
| Fig. 2(a-d): | The evolution of spiral pattern at different time (T) (a) 150, (b) 300, (c) 700 and (d) 1000, δ = 0.4, ε = 2, μ = 0.8 and β = 1 and initial condition as in Eq. 16 |
| Fig. 3(a-c): | The evolution of spiral pattern at different time (T) (a) 2000, (b) 3000 and (c) 4000, δ = 0.4, ε = 2, μ = 0.8, β = 1 and initial condition as in Eq. 16 |
Notice that as the values of δ is varied in the range of 0 <δ < 0.5 (i.e., choose δ = 0.45, 0.3 and 0.05), the evolution of the system can be observed, which led to the formation of spiral pattern as in Fig. 4a-c. The diameter of spiral pattern increases as the values of δ is decreased.
Note that, it is very interesting to observe the dynamics of model (15) when δ = 0.05. The spiral structures emerge in Fig. 4c but with small defects in their patterns. The zone in which defects occur appear to oscillate and these oscillations seem to continue in the stationary structures. Refer to Sherratt et al. (2007) for the discussion on this matter.
| Fig. 4(a-c): | Two-dimensional numerical solutions of prey-predator system (a) δ = 0.45, u* = 3/10, v* = 21/40, (b) δ = 0.30, u* = 1/5, v* = 2/5 and (c) δ = 0.05, u* = 1/30, v* = 29/360 using finite difference discretization at T = 150 and Δt = 1/384, ε = 2, μ = 0.8, β = 1 and initial condition as in Eq. 16 |
Apart from that, Garvie (2007) also conducted the same experiment as in Fig. 4c. But, he studied the dynamics of system (15) in one-dimensional case and discovered the occurrence of chaos covering almost the entire domain when δ = 0.05. Based on Bhattacharyay (2001), for two-dimensional case, spiral structures which are generally found to form around some defects are known as spiral defect chaos.
Another interesting phenomenon to be studied is the occurrence of spatio-temporal dynamics of model (15) through Turing instability. Consider the snapshots of spatio-temporal pattern below.
Based on Fig. 5, as the values of diffusion constant β increases, the formation of spiral pattern is observed. The spiral structures get larger, or in other words, its diameter increases as the values of β is varied in the range of 1≤β≤15 (Fig. 5a-c).
| Fig. 5(a-e): | Two-dimensional numerical solutions of prey-predator system (a) β= 1, (b) β= 10, (c) β=15, (d) β= 50 and (d) β= 100 using finite difference discretization at T = 500, Δt = 1/384 and (u* = 6/35, v*, 116/245, δ= 0.4, ε= 2, µ = 0.8 and initial condition as in Eq. 16 |
Refer to Dubey et al. (2009) for the discussion on this matter.
But when the value of β is changed to 50, observe that the spiral pattern loses its stability (Fig. 5d). In order to verify this situation, the value of β is increased to 100 (Fig. 5e). Consequently, it is observed that the dynamics become unstable and spatio-temporal chaos develops.
The formation of chaotic spiral patterns is actually caused by Turing instability. As discussed in the Mathematical Model part, defining:
then Turing instability will occur when D>>1. Therefore, Turing instability will occur with a large choice of values of parameter β i.e., 50, 100 and etc.
The occurrence of such instability is very interesting to be observed and studied. This is because if the constant D is large enough, prey growth rate will reach negative values and prey population will be driven by predators to a very low level in those regions. In other words, where the prey density is at their maximal value, diffusion will lower the prey density at that point. Conversely, where the prey density is at their minimal value, diffusion will increase the prey density at that point (Alonso et al., 2002).
CONCLUSION
Numerous spatio-temporal dynamical structures of diffusive prey-predator model viz., spiral waves, patchy structures, spiral defect chaos and spatio-temporal chaos due to Turing instability are studied. These amazingly beautiful structures emerge due to the proper selection of model parameters and suitable initial condition factors.
Hence by changing such factors, a number of future works and open problems remain. First, it is very interesting to investigate the formation of spatio-temporal dynamics with other choices of parameter. Next, other different forms of “disturbed” initial conditions can be considered so as to induce the non-trivial spatio-temporal patterns.
Another aspect to be studied is about the choice of boundary condition. In this study, the Neumann boundary condition as in Eq. 3 is chosen, but it is really interesting to investigate the interactions of prey-predator with other type of boundary conditions, namely, Dirichlet and Robin. The appearance of even more complicated phenomena is expected using other type of boundary conditions.
All these are very important to be observed and studied in order to understand the ecological interactions between prey and predator.
ACKNOWLEDGMENTS
Research supported by USM Fundamental Research Grant Scheme (FRGS).