- Email: [email protected]
Wave of chaos in a diffusive system: Generating realistic patterns of patchiness in plankton–fish dynamics
Chaos, Solitons and Fractals 40 (2009) 262–276 www.elsevier.com/locate/chaos
Wave of chaos in a diffusive system: Generating realistic patterns of patchiness in plankton–fish dynamics Ranjit Kumar Upadhyay b
a,*
, Nitu Kumari a, Vikas Rai
b
a Department of Applied Mathematics, Indian School of Mines University, Dhanbad 826 004, Jharkhand, India Department of Applied Mathematics, HMR Institute of Technology and Management, GT Karnal Road, Hamidpur, Delhi 110 03, India
Accepted 23 July 2007
Abstract We show that wave of chaos (WOC) can generate two-dimensional time-independent spatial patterns which can be a potential candidate for understanding planktonic patchiness observed in marine environments. These spatio-temporal patterns were obtained in computer simulations of a minimal model of phytoplankton–zooplankton dynamics driven by forces of diffusion. We also attempt to figure out the average lifetimes of these non-linear non-equilibrium patterns. These spatial patterns serve as a realistic model for patchiness found in aquatic systems (e.g., marine and oceanic). Additionally, spatio-temporal chaos produced by bi-directional WOCs is robust to changes in key parameters of the system; e.g., intra-specific competition among individuals of phytoplankton and the rate of fish predation. The ideas contained in the present paper may find applications in diverse fields of human endeavor. 2007 Elsevier Ltd. All rights reserved.
1. Introduction There is growing interest in the chaotic dynamics of ecological systems and also a need to study the problem of ecological dynamics in space, as well as time [1]. One of the most exciting and well-debated issues of modern ecology is the observation of plankton-patchiness [2]. Plankton acts as the base of the food chain in aquatic systems and plays an important role in the ecology of the ocean. The study of mechanisms which result in patchiness of the plankton populations can have significant impact on the fisheries policy [3]. Over the past few decades, many researches have been initiated to find the causes of such plankton patchiness. Researchers have demonstrated that spatio-temporally varying plankton dynamics can arise from a number of different sources [4]. Abraham [5] found that the change in the environment due to turbulent advection causes characteristic spatial patterns of the phytoplankton and zooplankton populations. Turing’s idea [6] that differential-diffusion can act on a reacting system to produce time-independent spatial patterns in a homogeneous environment was explored by Segel and Jackson [7] and by Levin and Segel [8] in an ecological context in the first part of seventies. The Levin and Segel [8] showed that this could serve as a mechanism for origin of planktonic patchiness in marine ecosystems. Turing spatial patterns have been observed in computer simulations of interaction–diffusion system by many authors. These structures *
Corresponding author. Tel.: +91 326 2200817; fax: +91 326 2210028. E-mail addresses: [email protected] (R.K. Upadhyay), [email protected] (N. Kumari).
0960-0779/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2007.07.078
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
263
of instability are generated under conditions of differential-diffusion for biological species like phytoplankton and zooplankton. This is not the case when dispersal of species is due to turbulent mixing. The more important observation in this context is that spatial patterns generated due to Turing instability are regular and stationary. Medvinsky et al. [9] examined plankton spatial patterns caused by diffusive interaction between nearby habitats some of which were fish free in a reaction–diffusion model with passive movement of individuals of two interacting species caused by turbulent mixing. Significant studies on plankton patchiness using conceptual minimal models have been discussed in the literature [10–13]. Scheffer [14] suggested a minimal phytoplankton–zooplankton model in which planktonic fronts propagation, predator–prey limit cycle oscillations and planktonic Turing patches were observed by Malchow [15,16]. In the same model system without fish predation, Pascual [17] observed diffusion-induced chaos with equal diffusivity constants. Recently, Tikhonov et al. [34] has investigated a hybrid model of spatio-temporal continuous phytoplankton–zooplankton with discrete agent planktivorous fish dynamics for inherent capacities of fish school chaotic movement. Non-stationary irregular patterns corresponding to spatio-temporal chaos were first observed in an interaction–diffusion system by Petrovskii and Malchow [18]. Medvinsky et al. [19] provide a lucid description of these non-Turing spatial patterns and conclude that these spatial patterns could serve as model for those commonly observed in marine systems. The chaos in this system is generated by propagating interfaces (of regions displaying regular and chaotic dynamics) in different directions. These authors termed this mechanism of spatio-temporal chaos generation as ‘‘wave of chaos’’. In this paper, we study an interaction–diffusion system, which is different from the model considered by Petrovskii and Malchow [18]. In the model system, the per capita growth rate of phytoplankton is a function of the nutrient level in the system. We have studied both one-dimensional (1D) and two-dimensional (2D) reaction–diffusion model of nutrient–phytoplankton–zooplankton–fish interaction. The chaotic and regular dynamics of the minimal model is studied in a patchy environment. The basic model was first introduced by Scheffer [14]. These conceptual minimal models are key to our understanding of Non-Turing time-independent spatial patterns.
2. The model system We consider a four-component nutrient–phytoplankton–zooplankton–fish model where at any location (X, Y) and time T, the phytoplankton P(X, Y, T) and zooplankton Z(X, Y, T) populations satisfy the reaction–diffusion equations oP aN tP P bP 2 Z þ DP DP ; ¼ oT ðH N þ NÞ ðH P þ P Þ
ð1aÞ
oZ etP FZ 2 þ DZ DZ; Z dZ 2 ¼ oT ðH P þ P Þ ðH Z þ Z 2 Þ
ð1bÞ
where DP and DZ are the diffusion coefficients of phytoplankton and zooplankton respectively. N is the nutrient level of the system. a is the maximum per capita growth rate of prey population. HN is the phytoplankton density at which specific growth rate becomes half its saturation value. b is the intensity of competition among individuals of phytoplankton. t is the rate at which phytoplankton is grazed and it follows Holling type II functional response. F is the predation rate of zooplankton by fish population, which follows Holling type III functional response. d is the mortality rate of predator while e is the conversion coefficient from individuals of phytoplankton into individuals of zooplankton. We assume the system is defined on a bounded domain (habitat) and is augmented with appropriate initial conditions and the zero-flux boundary conditions. The particular choice of boundary conditions reflects the assumption that the individual species cannot leave the domain. This model was originally formulated as a system of ordinary differential equations by Scheffer [14] and since then it has been studied by many authors [19–24]. Interactions involved in the minimal model system (1a) and (1b) are pictorially represented in Fig. 1. The interaction part of the model system investigated is depicted by solid arrows. The interaction involving dotted arrows indicates that either of the two effects involved, positive or negative is considered in the model. For example the dotted arrow between fish and zooplankton indicates that only the negative effect that fish have on the zooplankton of the system is considered in our model. These minimal models exhibit a wide range of ecologically
Fig. 1. Interactions incorporated in the model system.
264
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
relevant behavior such as spiral waves [19], target waves [25], diffusion-induced instability [22], and chaos [17,19]. Malchow et al. [26] presented a historical overview of modeling plankton dynamics and pattern formation mechanisms. We make the following assumptions for simplicity: r¼
aN ; HN þ N
r ¼ x; b
hðkP Þ ¼
P ; HP þ P
gðmZ 2 Þ ¼
ðH 2Z
Z2 : þ Z2Þ
Now we define dimensionless phytoplankton density, zooplankton density, spatial co-ordinates and time by 1=2 t P r ; xi ¼ X i ; t ¼ rT u¼ ; v¼Z x rx DP and rescaling the parameters via a ¼ kx;
b¼
et ; r
d c¼ ; r
d¼
mx ; t
D¼
DZ ; DP
leads to the following dimensionless system: ou ¼ uð1 uÞ vhðgÞ þ Du; ot ov ¼ bvhðgÞ cv Fgðf2 Þ þ DDv; ot
ð2aÞ ð2bÞ
with hðgÞ ¼
g ðH P þ gÞ
ðg ¼ auÞ;
gðf2 Þ ¼
f2 ðH 2z þ f2 Þ
ðn2 ¼ dv2 Þ;
ð3Þ
where u(x, y, t) and v(x, y, t) are phytoplankton and zooplankton densities in the patch respectively. x is the carrying capacity of an individual patch, b is the parameter measuring the ratio of product of conversion coefficient with grazing rate to the product of intensity of competition among individuals of phytoplankton with carrying capacity. c is the per capita predator’s death rate.
3. Methodology To investigate the spatio-temporal dynamics of the model system, we have numerically solved Eqs. (2a) and (2b) using semi-implicit (in time) finite difference method. We have called this method semi-implicit because right hand side of the method involves approximations at nth and (n 1)th time level. The steps Dx and Dt of the numerical grid are chosen sufficiently small so that the results do not depend upon the value of the steps. This method finally results to a sparse, banded linear system of algebraic equations. The linear system obtained, is solved by using LU factorization method for one-dimensional case and by GMRES algorithm for two-dimensional case. The graph between space and prey–predator density is obtained to study the spatial dynamics of the model system. Also temporal dynamics is studied by observing the effect of time on space vs. density plot of prey and predator populations. We have also used in our study the phase space representation and two-dimensional (2D) parameter scans. 2D scans are studied to identify the parameter regimes in which different dynamics (stable focus, stable limit cycle and chaos) exist. The parameters chosen for 2D scan diagrams are the parameters which control the dynamics of the model system. For two-dimensional case, we have obtained snapshots of prey and predator densities at different times and tried to study the spatio-temporal dynamics of the interaction–diffusion model.
4. Analytical results Before proceeding to the study of spatio-temporal pattern formation, it is reasonable to first consider the local dynamics of the system. Therefore in this section, we present some details of the analytical results obtained on the model system (1a) and (1b) without diffusion, which helps to make appropriate choice of parameter values for computer simulation.
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
265
4.1. Kolmogorov analysis For the model system to be biologically meaningful, the interaction part of the system dP aN tP P bP 2 Z; ¼ dT ðH N þ N Þ ðH P þ P Þ
ð4aÞ
dZ etP FZ 2 Z dZ 2 ; ¼ dT ðH P þ P Þ ðH Z þ Z 2 Þ
ð4bÞ
should be a Kolmogorov system. An application of Kolmogorov theorem for system (4a) and (4b) gives the following conditions H N aðet dÞ 1 ðiÞ et > d and ðiiÞ < N bdH P with ðaÞ Z 2 < H 2Z ; ðbÞ
etH P P 2
ðH P þ P Þ
>
FZ ðH 2Z Z 2 Þ ðH 2Z þ Z 2 Þ2
;
ðcÞ tZ H P þ bP ðH P þ P Þ2 > 0; where (P*, Z*) is the non-trivial equilibrium point. 4.2. Linear stability analysis An application of stability analysis to the Kolmogorov system (4a) and (4b), implies that the system possesses three stationary states which are given as under: (i) The total extinction point E0 = (0, 0) always exists and it is a saddle point. (ii) The equilibrium point E1 ¼ ðP ; 0Þ with P ¼ bðHaN , where only the extinction of the predator is seen always exists N þN Þ and is locally asymptotically stable. (iii) The non-trivial equilibrium point E* = (P*, Z*). From Eq. (4a), we obtain Z ¼ 1t ðr bP ÞðH P þ P Þ where r ¼ ðH NaNþN Þ. Substituting this Z* in Eq. (4b), we get the quantic polynomial A1P5 + A2P4 + A3P3 + A4P2 + A5P + A6 = 0. where A1 ¼ b2 ðet dÞ; A2 ¼ 2bðr bH P Þðet dÞ þ b2 H P d; A3 A4 A5 A6
¼ ðet dÞðr2 þ b2 H 2P 4rbH P Þ þ 2bdH P ðr bH P Þ; ¼ 2rH P ðr bH P Þðet dÞ dH P ðr2 þ b2 H 2P þ 4rbH P Þ; ¼ ðet dÞðr2 H 2P þ t2 H 2Z Þ 2rdH 2P ðr bH P Þ; ¼ dH P ðr2 H 2P þ t2 H 2Z Þ:
P* is a real positive root of the polynomial. Since all the parameters are positive and using the results of Kolmogorov analysis, we have A1 > 0 and A6 < 0. For r < bHP and 2b(r bHP)(et d) > b2HPd, we have A2 < 0 and A5 > 0. Therefore irrespective of the sign of A3 and A4, there are two positive real roots of the above quantic polynomial. P* is one of the two real positive roots of the polynomial. We now discuss the local behavior of system of ODEs (4a) and (4b) around each of the three equilibrium points. The variational matrix at any point E(P, Z) is given by 3 2 aN tZH P tP 2bP 7 6 ðH N þ N Þ ðH P þ P Þ ðH P þ P Þ2 7 6 V ðP ; ZÞ ¼ 6 7: 4 2FZH 2Z 5 etH P Z etP d 2 2 2 ðH P þ P Þ ðH P þ P Þ ðH z þ Z 2 Þ From the above general variational matrix, the variational matrix about the equilibrium point E0 is given by " # aN 0 ðH N þN Þ : V0 ¼ 0 d
266
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
It has eigenvalues l01 ¼ H NaNþN and l02 = d which shows that E0 is a saddle point with an unstable manifold along X-direction and a stable manifold along Y-direction. The variational matrix about the point E1 is 3 2 aN taN bH ðH þ NÞ þ aN 7 6 P N V 1 ¼ 4 ðH N þ N Þ 5: etaN d 0 bH P ðH N þ N Þ þ aN The eigenvalues are l11 ¼ H NaNþN and l12 ¼ bH P ðHetaN d. Therefore the equilibrium point E1 is locally asymptotiN þN ÞþaN d < 0. cally stable provided bH P ðHetaN N þN ÞþaN The variational matrix for non-trivial point E* = (P*, Z*) is 3 2 tP Z tP bP 6 ðH þ P Þ2 ðH P þ P Þ 7 7 6 P V ¼6 7: 4 FZ ðZ 2 H 2Z Þ 5 etH P Z ðH P þ P Þ2
ðZ 2 þ H 2Z Þ2
Theorem . The unique non-zero positive equilibrium point E* is locally asymptotically stable provided the following condition is satisfied. tZ P 2
ðH P þ P Þ
<
FZ ðH 2Z Z 2 Þ ðH 2Z þ Z 2 Þ2
þ bP :
Proof . The characteristic equation of the variational matrix V* is given by k2 þ Ak þ B ¼ 0; where A¼ B ¼ FZ
tZ P
þ 2
FZ ðZ 2 H 2Z Þ
! bP ;
ðH 2Z þ Z 2 Þ2 ! tZ P ðZ 2 H 2Z Þ et2 H P P Z bP þ : 2 ðH P þ P Þ ðH 2Z þ Z 2 Þ2 ðH P þ P Þ3
ðH P þ P Þ
Applying Routh–Hurwitz criteria, E* is locally asymptotically stable provided A > 0 and B > 0 which gives tZ P ðH P þP Þ2
<
FZ ðH 2Z Z 2 Þ ðH 2Z þZ 2 Þ2
þ bP , where (P*, Z*) is the non-trivial equilibrium point.
h
5. Simulation results To better understand the spatial dynamics of the proposed diffusion model, we turn to its numerical simulation. The computer simulation of the model system (2a) and (2b) was performed using MATLAB 7.0. 5.1. The spatio-temporal dynamics in one-dimensional case In this subsection, we focus on the one-dimensional dynamics of the system (2a) and (2b) with u(x, t), v(x, t) and D = o2/ox2. Eqs. (2a) and (2b) takes the following form: ou vu o2 u þ 2; ¼ uð1 uÞ ot ðH P þ uÞ ox
ð5aÞ
ov vu v2 o2 v cv F þD 2: ¼b 2 2 ot ðH P þ uÞ ox ðH Z þ v Þ
ð5bÞ
In Eqs. (5a) and (5b), 0 < x < Lx, 0 < y < Ly. Model system (5a) and (5b) describes the dynamics of an aquatic system in a horizontal layer, vertical distribution of species inside the layer is assumed to be homogeneous. In our simulation, we
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
267
have considered a square domain of dimension Lx = Ly = 400. To be more realistic from the biological point of view, we consider a complex and non-monotic form of initial condition which determines the initial spatial distribution of the species in the real community as under uðx; 0Þ ¼ u þ eðx x1 Þðx x2 Þ; vðx; 0Þ ¼ v ;
ð6Þ
where (u*, v*) is the non-trivial state for the co-existence of prey and predator. e is the parameter which affects the system dynamics. At the domain boundary, we consider the zero flux boundary conditions: ou ou ov ov ¼ ¼ ¼ ¼ 0: ð7Þ oxð0;tÞ oxðL;tÞ oxð0;tÞ oxðL;tÞ This type of boundary condition is specifically used for modeling the dynamics of spatially bounded aquatic ecosystems [27]. We call a part of space a patch if it has a distinct kind of local spatial dynamics than adjacent parts. WOC propagates in both forward as well as in backward directions replacing a patch with local regular dynamics (characterized by a stable limit cycle in the phase plane of the system) with chaotic dynamics. These patches are shown in Fig. 2. These patches [9] are formed as a result of the action of diffusion on the interacting biological species (phytoplankton and zooplankton in the present case). With initial condition (6), when computer simulation is performed, a complex spatial dynamics is observed at t = 1000 as shown in Fig. 4a where a smooth regular patch separates chaotic patches at x1 = 1200 and x2 = 2800 with other parameters as in Fig. 3. We found a set of parameter values HP = 0.3, b = 2, c = 0.8, F = 1 · 104, HZ = 2.5, e = 1 · 108 where chaotic dynamics is observed with an appearance of jagged pattern at t = 600 as shown in Fig. 3c. Slightly changing the parameter values of b and F, keeping the remaining parameters same as above stable focus and limit cycle are obtained as in Figs. 3a and 3b respectively. At t = 600, numerical simulation for following set of parameter values HP = 0.3, b = 1.1, c = 0.8, F = 0.2, HZ = 2.5, u* = 0.6908, v* = 0.1566 gives a stable focus as shown in Fig. 3a. While taking b = 2, F = 1 · 103 and keeping remaining parameters same as above we get a limit cycle solution (cf. Fig. 3b). We obtained chaotic behavior for the parameter values HP = 0.3, b = 2, c = 0.8, F = 1 · 104, HZ = 2.5 and the corresponding (u*, v*) = (0.2, 0.4). Then after fixing these parameters, we tried to study the effect of time on the chaotic dynamics of the model system. We observed from our simulation results that the chaotic dynamics persists for a long period of time. The spatial distributions vary gradually in time after the local temporal behavior of the dynamical variables u and v follow a stable limit cycle of the homogeneous system (5a) and (5b) (cf. Fig. 4). The one-dimensional analysis shows the formation of strongly irregular transient pattern inside a subdomain. These jagged chaotic patterns and regular smooth patterns occupy the space alternatively. From Fig. 4, we observed that this jagged pattern representing chaotic behavior of the system grows steadily with time. The size of the domain occupied by the irregular chaotic patterns slowly grows with time in both directions equally, displacing the regular pattern and hence occupying the whole region. This phenomenon has been termed as ‘‘wave of chaos’’ (WOC) [18]. This phenomenon persists for the long period of time and after the irregular pattern occupies the whole domain, the dynamics of the system does not undergo any further changes. For the set of parameter values of Fig. 3c, the temporal behavior of prey and predator abundances are also chaotic as shown in Fig. 5. The phase plane is obtained at fixed points x1 = 1200 and x2 = 2800 of the domain. By changing the value of the parameter b in the above set of parameters, a limit cycle is obtained at b = 1.9 shown in Fig. 6 and a stable focus is observed at b = 1.2.
Fig. 2. Patchy structure of the system. Individual patches have distinct local dynamics in the beginning. Patches with regular spatial dynamics are constantly replaced by ones with chaotic dynamics. WOC generates chaos in the whole domain.
268
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276 0.9
1.4 prey predator
0.8
1.2
0.7
1
0.6
prey,predator
prey,predator
prey predator
0.5 0.4
0.8
0.6
0.4 0.3 0.2
0.2 0.1
0
0.2
0.4
0.6
0.8
1 space
1.2
1.4
1.6
1.8
2
0
0
0.2
0.4
0.6
0.8
4
x 10
1 space
1.2
1.4
1.6
1.8
2 4
x 10
1.4 prey predator
1.2
prey,predator
1
0.8
0.6
0.4
0.2
0
0
0.2
0.4
0.6
0.8
1 space
1.2
1.4
1.6
1.8
2 4
x 10
Fig. 3. Spatial distribution of the population showing the onset of chaotic phase after gradually displacing the regular phase at t = 600 of model (5a) and (5b) for parameter values HP = 0.3, c = 0.8, HZ = 2.5, Dt = 101, h = 1. (a) b = 1.1, F = 0.2, (b) b = 2, F = 0.001, (c) b = 2, F = 0.0001.
2D scan results on the parameter space (c, b) and (F, b) has been presented to observe simultaneous effect of these parameters on the dynamics of the model system (5a) and (5b) as shown in Fig. 7. From Fig. 7, it is found that the chaotic dynamics is observed at discrete points for b = 2 and in the range 0 6 F 6 0.06 (Fig. 7b). In the parameter space (c, b) chaos is observed along the line b = c. We performed extensive numerical computations to test whether spatio-temporal chaos generated through the action of WOC is robust with respect to changes in system parameters. It is evident from these figures that the chaotic dynamics is robust to changes for the range 0 6 F 6 0.06 at b = 2 in parameters space (F, b). 5.2. The spatio-temporal dynamics in two-dimensional case In this subsection, we consider dynamics of the two-dimensional model system (2a) and (2b) with u = u(x, y, t), v = v(x, y, t) and D = o2/ox2 + o2/oy2. Eqs. (2a) and (2b) takes the following form: ou vu o2 u o2 u þ 2þ 2; ¼ uð1 uÞ ot ðH P þ uÞ ox oy
2 ov vu v2 o v o2 v cv F ; þ þ D ¼b ot ðH P þ uÞ ox2 oy 2 ðH 2Z þ v2 Þ
ð8aÞ ð8bÞ
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276 1.4
1.4 prey predator
1.2
prey,predator
prey,predator
1
0.8
0.6
0.8
0.6
0.4
0.4
0.2
0.2
0
1000
2000
3000 4000 space
5000
6000
0
7000
0
prey predator
1.2
2000
3000 4000 space
5000
6000
7000
prey predator
1.2
1
prey,predator
1
0.8
0.6
0.8
0.6
0.4
0.4
0.2
0.2
0
1000
1.4
1.4
prey,predator
prey predator
1.2
1
0
269
0
1000
2000
3000 4000 space
5000
6000
0
7000
0
1000
2000
3000
4000 space
5000
6000
7000
1.4 prey predator
1.2
prey,predator
1 0.8 0.6 0.4 0.2 0
0
1000
2000
3000
4000 space
5000
6000
7000
Fig. 4. Simulation showing the long-term persistence of chaotic nature as the regular phase is gradually displaced by the chaotic phase (a) t = 1000 (b) t = 2000 (c) t = 5000 (d) t = 10,000 (e) t = 15,000.
The spatio-temporal dynamics of a chaotic system largely depends upon the choice of initial conditions [19]. We present the results of computer simulations with two different initial spatial distributions. First simulation result of two-dimensional case, the initial distribution of the species considered is a two-dimensional generalization of initial condition given in Eq. (6) (note that the initial conditions are deliberately chosen so that it is asymmetrical in order to make more influence of the corners of the domain). It is of the following form: uðx; y; 0Þ ¼ u e1 ðx 0:1y 225Þðx 0:1y 675Þ; vðx; y; 0Þ ¼ v e2 ðx 450Þ e3 ðy 150Þ:
ð9Þ
270
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276 1.4
1.2
Predator
1
0.8
0.6
0.4
0.2
0
0
0.1
0.2
0.3
0.4 Prey
0.5
0.6
0.7
0.8
Fig. 5. Phase plane of the system for the spatial range defined by fixed points x1 = 1200, x2 = 2800 inside the domain occupied by irregular spatio-temporal oscillation. Parameters are given in the text. 1 0.9 0.8 0.7
Predator
0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4 Prey
0.5
0.6
0.7
0.8
Fig. 6. Phase plane of the system showing limit cycle behavior with b = 1.9. Other parameters are same as in Fig. 4.
3.5
2.9
3
2.6 2.3
2.5
2
2
1.7
b 1.5
b 1.4
1
1.1
0.5
0.8
0
0.5
0
0.5
1
c
1.5
2
2.5
0
0.02
0.04
0.06
0.08
0.1
F
Fig. 7. 2D scan diagrams for the model system (5a) and (5b) in the parameter space (c, b) and (F, b).
The numerical results obtained using semi-implicit numerical techniques of two-dimensional phytoplankton–zooplankton system are presented in Fig. 8. The parameter values at which intermittent chaos was obtained in 1D analysis and initial condition (9) with zero-flux boundary conditions (as explained in 1D analysis) in two dimension are used to perform simulations of two-dimensional spatio-temporal system. The parameter values taken are, HP = 0.3, b = 2, c = 0.8, HZ = 2.5, F = 1 · 104, e1 = 2 · 107, e2 = 3 · 105, e3 = 1.2 · 104 Dt = 1/3, h = 1, d = 1.
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
271
Phytoplankton population density is plotted in first column and zooplankton population density in the second column of Fig. 8 at t = 0, 150, 200, 300, 600 and 1000. Fig. 8 clearly depicts the variations in the population behavior of phytoplankton and zooplankton at the early stages i.e. at t = 0, 150, 200. But at later stages, i.e. at t = 300, 600, 1000
Fig. 8. Spatial distribution of prey densities (phytoplankton) [first column] and predator densities (zooplankton) [second column] of model system (8a) and (8b) with initial condition given in Eq. (9) for (a) t = 0, (b) t = 150, (c) t = 200, (d) t = 300, (e) t = 600, (f) t = 1000. Parameter values are given in the text.
272
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
Fig 8. (continued)
both prey and predator densities exhibit almost similar spatial population dynamics. The regular spiral spatial pattern is followed by the irregular patchy patterns. Absence of spiral patterns at t = 0 justifies that the appearance of spirals is not induced by initial condition. These spiral patterns first grow till a particular stage after which its destruction begins from their centers. Finally an irregular pattern prevails over the entire square domain at t = 600 which persists with time as seen at t = 1000.
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
273
Fig. 9. Spatial distribution of prey densities (phytoplankton) [first column] and predator densities (zooplankton) [second column] using initial condition given in Eq. (10) for (a) t = 0, (b) t = 150, (c) t = 200, (d) t = 300, (e) t = 600, (f) t = 1000. Parameter values are given in the text.
Although choice of the initial condition drastically affects the dynamics of the system but the spatial dynamics observed in Fig. 8 is not induced by initial condition. To justify this fact, we consider another initial condition for the second set of simulation, which is
274
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
Fig 9. (continued)
uðx; y; 0Þ ¼ u e1 ðx 180Þðx 720Þ e2 ðy 90Þðy 210Þ; vðx; y; 0Þ ¼ v e3 ðx 450Þ e4 ðy 135Þ and the parameters are the same as those used to generate Fig. 8.
ð10Þ
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
275
Similar scenario can be seen in Fig. 9. An irregular patchy pattern is finally obtained preceded by a somewhat less regular spiral pattern. We observed that the dynamics of the system is regular only for a specified period of time after which an irregular dynamics prevails. Also the chaotic nature of the minimal diffusion model has been studied and it is observed that spatio-temporal chaotic dynamics of this interaction–diffusion system is self-organized.
6. Discussions and conclusion These spatial patterns can be considered as time-independent as after certain amount of time is elapsed, there were not visible changes in spatial patterns (cf. Figs. 8e and 8f, 9e and 9f). The spatio-temporal chaos generated by propagation of WOC is robust to changes in key system parameters such as intensity of competition between individuals of phytoplankton species and changes in the per capita rate of predation by fishes. Two-dimensional spatial patterns of phytoplankton–zooplankton dynamics are self-organized and, therefore, can be considered to provide a theoretical framework to understand patchiness in marine environments. Recently, Medvinsky and collaborators [28] have proposed a new mechanism for the origin of spatio-temporal chaos in a heterogeneous environment wherein nearby habitats possess different kind of local dynamics. It remains to be seen which of these mechanisms is more common in marine and other aquatic systems. The discovery of spatio-temporal chaos [29] in early nineties has potential to lead to many noble applications in Engineering [30,31] and Physiology [32,33]. There exist many mechanisms, which cause spatio-temporal chaos in a given system. It is imperative to assess their potential to solve problems in Physics, Biology and Engineering. Here we attempt to do the same in a model ecological system.
Acknowledgement This research work is supported by the All India Council for Technical Education which has offered the award of AICTE National Doctoral Fellowship to Nitu Kumari (F.No.:1-10/FD/NDF-PG/(29)/2006-07).
References [1] Rai V, Schaffer WM. Chaos in ecology. Chaos, Solitons & Fractals 2001;12:197–203. [2] Franks PJS. Spatial patterns in dense algal blooms. Limnol Oceanol 1997;42:1297–305. [3] Legendre L. The significance of microalgal blooms for fisheries and for the export of particulate organic carbon in Oceans. J Plank Res 1990;12:681–99. [4] Hillary RM, Bees MA. Plankton lattices and the role of chaos in plankton patchiness. Phys Rev E 2004;69:031913. [5] Abraham ER. The generation of plankton patchiness by turbulent stirring. Nature 1998;391:577–80. [6] Turing AM. On the chemical basis of morphogenesis. Philos Trans R Soc London Ser B 1952;237:37–72. [7] Segel LA, Jackson JL. Dissipative structures: an explanation and an ecological example. J Theor Biol 1972;37:545–59. [8] Levin SA, Segel LA. Hypothesis for origin of planktonic patchiness. Nature 1976;259:659. [9] Medvinsky AB, Tikhonova IA, Aliev RR, Li BL, Lin ZS, Malchow H. Patchy environment as a factor of complex plankton dynamics. Phys Rev E 2001;64:021915–7. [10] Fasham MJR. The statistical and mathematical analysis of plankton patchiness. Oceanogr Mar Biol Ann Rev 1978;16:43–79. [11] Franks PJS, Wroblewski JS, Fileri GR. Behavior of simple plankton model with food-level acclimation by herbivores. Mar Biol 1986;91:121–9. [12] Steel JH, Hunderson EW. A simple model for plankton patchiness. J Plankton Res 1992;14:1397–403. [13] Truscott JE, Brindley J. Equilibria, stability and excitability in a general class of plankton population models. Philos Trans R Soc London Ser A 1994;347:703–18. [14] Scheffer M. Fish and nutrients interplay determines algal biomass: a minimal model. OIKOS 1991;62:271–82. [15] Malchow H. Spatio-temporal pattern formation in non-linear nonequilibrium plankton dynamics. Proc R Soc London Ser B 1993;251:103–9. [16] Malchow H. Nonequilibrium structures in plankton dynamics. Ecol Model 1994;75–76:123–34. [17] Pascual M. Diffusion induced chaos in a spatial predator–prey system. Proc R Soc London Ser B 1993;103:251–7. [18] Petrovskii SV, Malchow H. Wave of chaos: new mechanism of pattern formation in spatio-temporal population dynamics. Theor Populat Biol 2001;59:157–74. [19] Medvinsky A, Petrovskii S, Tikhonova I, Malchow H, Li Bai-Lian. Spatiotemporal complexity of plankton and fish dynamics. SIAM Rev 2002;44:311–70. [20] Malchow H, Petrovskii SV, Medvinsky AB. Numerical study of plankton–fish dynamics in a spatially structured and noisy environment. Ecol Model 2002;149:247–55.
276
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
[21] Malchow H. Motional instabilities in prey–predator systems. J Theor Biol 2000;204:639–47. [22] Malchow H, Hilker F, Petrovskii S. Noise and productivity dependence of spatiotemporal pattern formation in a prey–predator system. Discr Cont Dyn – B 2004;4:705–11. [23] Malchow H, Radtke B, Kallache M, Medvinsky A, Tikhonov D, Petrovskii S. Spatio-temporal pattern formation in coupled models of plankton dynamics and fish school motion. Non-linear Anal – Real World Appl 2000;1:53–67. [24] Tikhonova I, Li B, Malchow H, Medvinsky A. The impact of the phytoplankton growth rule on spatial and temporal dynamics of plankton communities in a heterogeneous environment. Biofizika 2003;48:891–9. [25] Sherratt J, Eagan B, Lewis M. Oscillations and chaos behind predator–prey invasion: mathematical artifact or ecological reality? Philos Trans R Soc London B 1997;352:21–38. [26] Malchow H, Petrovskii S, Medvinsky AB. Pattern formation in models of plankton dynamics: a synthesis. Oceanol Acta 2001;24:479–87. [27] Scheffer M. Population and community biology series 22. Ecology of shallow lakes. London: Chapman & Hall; 1998. [28] Medvinsky AB. Personal communication, 2006. [29] Cross MC, Hohenberg PC. Spatio-temporal chaos. Science 1994;263:1569–70. [30] Hu G, Xiao J, Yang J, Xie F, Qu Z. Synchronization of spatio-temporal chaos and its applications. Phys Rev E 1997;56(3):2738–46. [31] Gang H, Zhilin Q. Controlling spatiotemporal chaos in coupled map lattices. Phys Rev Lett 1994;72:68–71. [32] Zhang H, Patel N. Spiral Wave breakdown in an excitable medium model of cardiac tissue. Chaos, Solitons & Fractals 1995;5(3– ):635–9. [33] Bramann Y, John F, Lindner JF, Ditto WL. Taming spatio-temporal chaos with disorder. Nature 1995;378:465–7. [34] Tikhonov DA, Enderlein J, Malchow H, Medvinsky AB. Chaos and fractals in fish school motion. Chaos, Solitons & Fractals 2001;12:277–88.
Wave of chaos in a diffusive system: Generating realistic patterns of patchiness in plankton–fish dynamics Ranjit Kumar Upadhyay b
a,*
, Nitu Kumari a, Vikas Rai
b
a Department of Applied Mathematics, Indian School of Mines University, Dhanbad 826 004, Jharkhand, India Department of Applied Mathematics, HMR Institute of Technology and Management, GT Karnal Road, Hamidpur, Delhi 110 03, India
Accepted 23 July 2007
Abstract We show that wave of chaos (WOC) can generate two-dimensional time-independent spatial patterns which can be a potential candidate for understanding planktonic patchiness observed in marine environments. These spatio-temporal patterns were obtained in computer simulations of a minimal model of phytoplankton–zooplankton dynamics driven by forces of diffusion. We also attempt to figure out the average lifetimes of these non-linear non-equilibrium patterns. These spatial patterns serve as a realistic model for patchiness found in aquatic systems (e.g., marine and oceanic). Additionally, spatio-temporal chaos produced by bi-directional WOCs is robust to changes in key parameters of the system; e.g., intra-specific competition among individuals of phytoplankton and the rate of fish predation. The ideas contained in the present paper may find applications in diverse fields of human endeavor. 2007 Elsevier Ltd. All rights reserved.
1. Introduction There is growing interest in the chaotic dynamics of ecological systems and also a need to study the problem of ecological dynamics in space, as well as time [1]. One of the most exciting and well-debated issues of modern ecology is the observation of plankton-patchiness [2]. Plankton acts as the base of the food chain in aquatic systems and plays an important role in the ecology of the ocean. The study of mechanisms which result in patchiness of the plankton populations can have significant impact on the fisheries policy [3]. Over the past few decades, many researches have been initiated to find the causes of such plankton patchiness. Researchers have demonstrated that spatio-temporally varying plankton dynamics can arise from a number of different sources [4]. Abraham [5] found that the change in the environment due to turbulent advection causes characteristic spatial patterns of the phytoplankton and zooplankton populations. Turing’s idea [6] that differential-diffusion can act on a reacting system to produce time-independent spatial patterns in a homogeneous environment was explored by Segel and Jackson [7] and by Levin and Segel [8] in an ecological context in the first part of seventies. The Levin and Segel [8] showed that this could serve as a mechanism for origin of planktonic patchiness in marine ecosystems. Turing spatial patterns have been observed in computer simulations of interaction–diffusion system by many authors. These structures *
Corresponding author. Tel.: +91 326 2200817; fax: +91 326 2210028. E-mail addresses: [email protected] (R.K. Upadhyay), [email protected] (N. Kumari).
0960-0779/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2007.07.078
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
263
of instability are generated under conditions of differential-diffusion for biological species like phytoplankton and zooplankton. This is not the case when dispersal of species is due to turbulent mixing. The more important observation in this context is that spatial patterns generated due to Turing instability are regular and stationary. Medvinsky et al. [9] examined plankton spatial patterns caused by diffusive interaction between nearby habitats some of which were fish free in a reaction–diffusion model with passive movement of individuals of two interacting species caused by turbulent mixing. Significant studies on plankton patchiness using conceptual minimal models have been discussed in the literature [10–13]. Scheffer [14] suggested a minimal phytoplankton–zooplankton model in which planktonic fronts propagation, predator–prey limit cycle oscillations and planktonic Turing patches were observed by Malchow [15,16]. In the same model system without fish predation, Pascual [17] observed diffusion-induced chaos with equal diffusivity constants. Recently, Tikhonov et al. [34] has investigated a hybrid model of spatio-temporal continuous phytoplankton–zooplankton with discrete agent planktivorous fish dynamics for inherent capacities of fish school chaotic movement. Non-stationary irregular patterns corresponding to spatio-temporal chaos were first observed in an interaction–diffusion system by Petrovskii and Malchow [18]. Medvinsky et al. [19] provide a lucid description of these non-Turing spatial patterns and conclude that these spatial patterns could serve as model for those commonly observed in marine systems. The chaos in this system is generated by propagating interfaces (of regions displaying regular and chaotic dynamics) in different directions. These authors termed this mechanism of spatio-temporal chaos generation as ‘‘wave of chaos’’. In this paper, we study an interaction–diffusion system, which is different from the model considered by Petrovskii and Malchow [18]. In the model system, the per capita growth rate of phytoplankton is a function of the nutrient level in the system. We have studied both one-dimensional (1D) and two-dimensional (2D) reaction–diffusion model of nutrient–phytoplankton–zooplankton–fish interaction. The chaotic and regular dynamics of the minimal model is studied in a patchy environment. The basic model was first introduced by Scheffer [14]. These conceptual minimal models are key to our understanding of Non-Turing time-independent spatial patterns.
2. The model system We consider a four-component nutrient–phytoplankton–zooplankton–fish model where at any location (X, Y) and time T, the phytoplankton P(X, Y, T) and zooplankton Z(X, Y, T) populations satisfy the reaction–diffusion equations oP aN tP P bP 2 Z þ DP DP ; ¼ oT ðH N þ NÞ ðH P þ P Þ
ð1aÞ
oZ etP FZ 2 þ DZ DZ; Z dZ 2 ¼ oT ðH P þ P Þ ðH Z þ Z 2 Þ
ð1bÞ
where DP and DZ are the diffusion coefficients of phytoplankton and zooplankton respectively. N is the nutrient level of the system. a is the maximum per capita growth rate of prey population. HN is the phytoplankton density at which specific growth rate becomes half its saturation value. b is the intensity of competition among individuals of phytoplankton. t is the rate at which phytoplankton is grazed and it follows Holling type II functional response. F is the predation rate of zooplankton by fish population, which follows Holling type III functional response. d is the mortality rate of predator while e is the conversion coefficient from individuals of phytoplankton into individuals of zooplankton. We assume the system is defined on a bounded domain (habitat) and is augmented with appropriate initial conditions and the zero-flux boundary conditions. The particular choice of boundary conditions reflects the assumption that the individual species cannot leave the domain. This model was originally formulated as a system of ordinary differential equations by Scheffer [14] and since then it has been studied by many authors [19–24]. Interactions involved in the minimal model system (1a) and (1b) are pictorially represented in Fig. 1. The interaction part of the model system investigated is depicted by solid arrows. The interaction involving dotted arrows indicates that either of the two effects involved, positive or negative is considered in the model. For example the dotted arrow between fish and zooplankton indicates that only the negative effect that fish have on the zooplankton of the system is considered in our model. These minimal models exhibit a wide range of ecologically
Fig. 1. Interactions incorporated in the model system.
264
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
relevant behavior such as spiral waves [19], target waves [25], diffusion-induced instability [22], and chaos [17,19]. Malchow et al. [26] presented a historical overview of modeling plankton dynamics and pattern formation mechanisms. We make the following assumptions for simplicity: r¼
aN ; HN þ N
r ¼ x; b
hðkP Þ ¼
P ; HP þ P
gðmZ 2 Þ ¼
ðH 2Z
Z2 : þ Z2Þ
Now we define dimensionless phytoplankton density, zooplankton density, spatial co-ordinates and time by 1=2 t P r ; xi ¼ X i ; t ¼ rT u¼ ; v¼Z x rx DP and rescaling the parameters via a ¼ kx;
b¼
et ; r
d c¼ ; r
d¼
mx ; t
D¼
DZ ; DP
leads to the following dimensionless system: ou ¼ uð1 uÞ vhðgÞ þ Du; ot ov ¼ bvhðgÞ cv Fgðf2 Þ þ DDv; ot
ð2aÞ ð2bÞ
with hðgÞ ¼
g ðH P þ gÞ
ðg ¼ auÞ;
gðf2 Þ ¼
f2 ðH 2z þ f2 Þ
ðn2 ¼ dv2 Þ;
ð3Þ
where u(x, y, t) and v(x, y, t) are phytoplankton and zooplankton densities in the patch respectively. x is the carrying capacity of an individual patch, b is the parameter measuring the ratio of product of conversion coefficient with grazing rate to the product of intensity of competition among individuals of phytoplankton with carrying capacity. c is the per capita predator’s death rate.
3. Methodology To investigate the spatio-temporal dynamics of the model system, we have numerically solved Eqs. (2a) and (2b) using semi-implicit (in time) finite difference method. We have called this method semi-implicit because right hand side of the method involves approximations at nth and (n 1)th time level. The steps Dx and Dt of the numerical grid are chosen sufficiently small so that the results do not depend upon the value of the steps. This method finally results to a sparse, banded linear system of algebraic equations. The linear system obtained, is solved by using LU factorization method for one-dimensional case and by GMRES algorithm for two-dimensional case. The graph between space and prey–predator density is obtained to study the spatial dynamics of the model system. Also temporal dynamics is studied by observing the effect of time on space vs. density plot of prey and predator populations. We have also used in our study the phase space representation and two-dimensional (2D) parameter scans. 2D scans are studied to identify the parameter regimes in which different dynamics (stable focus, stable limit cycle and chaos) exist. The parameters chosen for 2D scan diagrams are the parameters which control the dynamics of the model system. For two-dimensional case, we have obtained snapshots of prey and predator densities at different times and tried to study the spatio-temporal dynamics of the interaction–diffusion model.
4. Analytical results Before proceeding to the study of spatio-temporal pattern formation, it is reasonable to first consider the local dynamics of the system. Therefore in this section, we present some details of the analytical results obtained on the model system (1a) and (1b) without diffusion, which helps to make appropriate choice of parameter values for computer simulation.
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
265
4.1. Kolmogorov analysis For the model system to be biologically meaningful, the interaction part of the system dP aN tP P bP 2 Z; ¼ dT ðH N þ N Þ ðH P þ P Þ
ð4aÞ
dZ etP FZ 2 Z dZ 2 ; ¼ dT ðH P þ P Þ ðH Z þ Z 2 Þ
ð4bÞ
should be a Kolmogorov system. An application of Kolmogorov theorem for system (4a) and (4b) gives the following conditions H N aðet dÞ 1 ðiÞ et > d and ðiiÞ < N bdH P with ðaÞ Z 2 < H 2Z ; ðbÞ
etH P P 2
ðH P þ P Þ
>
FZ ðH 2Z Z 2 Þ ðH 2Z þ Z 2 Þ2
;
ðcÞ tZ H P þ bP ðH P þ P Þ2 > 0; where (P*, Z*) is the non-trivial equilibrium point. 4.2. Linear stability analysis An application of stability analysis to the Kolmogorov system (4a) and (4b), implies that the system possesses three stationary states which are given as under: (i) The total extinction point E0 = (0, 0) always exists and it is a saddle point. (ii) The equilibrium point E1 ¼ ðP ; 0Þ with P ¼ bðHaN , where only the extinction of the predator is seen always exists N þN Þ and is locally asymptotically stable. (iii) The non-trivial equilibrium point E* = (P*, Z*). From Eq. (4a), we obtain Z ¼ 1t ðr bP ÞðH P þ P Þ where r ¼ ðH NaNþN Þ. Substituting this Z* in Eq. (4b), we get the quantic polynomial A1P5 + A2P4 + A3P3 + A4P2 + A5P + A6 = 0. where A1 ¼ b2 ðet dÞ; A2 ¼ 2bðr bH P Þðet dÞ þ b2 H P d; A3 A4 A5 A6
¼ ðet dÞðr2 þ b2 H 2P 4rbH P Þ þ 2bdH P ðr bH P Þ; ¼ 2rH P ðr bH P Þðet dÞ dH P ðr2 þ b2 H 2P þ 4rbH P Þ; ¼ ðet dÞðr2 H 2P þ t2 H 2Z Þ 2rdH 2P ðr bH P Þ; ¼ dH P ðr2 H 2P þ t2 H 2Z Þ:
P* is a real positive root of the polynomial. Since all the parameters are positive and using the results of Kolmogorov analysis, we have A1 > 0 and A6 < 0. For r < bHP and 2b(r bHP)(et d) > b2HPd, we have A2 < 0 and A5 > 0. Therefore irrespective of the sign of A3 and A4, there are two positive real roots of the above quantic polynomial. P* is one of the two real positive roots of the polynomial. We now discuss the local behavior of system of ODEs (4a) and (4b) around each of the three equilibrium points. The variational matrix at any point E(P, Z) is given by 3 2 aN tZH P tP 2bP 7 6 ðH N þ N Þ ðH P þ P Þ ðH P þ P Þ2 7 6 V ðP ; ZÞ ¼ 6 7: 4 2FZH 2Z 5 etH P Z etP d 2 2 2 ðH P þ P Þ ðH P þ P Þ ðH z þ Z 2 Þ From the above general variational matrix, the variational matrix about the equilibrium point E0 is given by " # aN 0 ðH N þN Þ : V0 ¼ 0 d
266
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
It has eigenvalues l01 ¼ H NaNþN and l02 = d which shows that E0 is a saddle point with an unstable manifold along X-direction and a stable manifold along Y-direction. The variational matrix about the point E1 is 3 2 aN taN bH ðH þ NÞ þ aN 7 6 P N V 1 ¼ 4 ðH N þ N Þ 5: etaN d 0 bH P ðH N þ N Þ þ aN The eigenvalues are l11 ¼ H NaNþN and l12 ¼ bH P ðHetaN d. Therefore the equilibrium point E1 is locally asymptotiN þN ÞþaN d < 0. cally stable provided bH P ðHetaN N þN ÞþaN The variational matrix for non-trivial point E* = (P*, Z*) is 3 2 tP Z tP bP 6 ðH þ P Þ2 ðH P þ P Þ 7 7 6 P V ¼6 7: 4 FZ ðZ 2 H 2Z Þ 5 etH P Z ðH P þ P Þ2
ðZ 2 þ H 2Z Þ2
Theorem . The unique non-zero positive equilibrium point E* is locally asymptotically stable provided the following condition is satisfied. tZ P 2
ðH P þ P Þ
<
FZ ðH 2Z Z 2 Þ ðH 2Z þ Z 2 Þ2
þ bP :
Proof . The characteristic equation of the variational matrix V* is given by k2 þ Ak þ B ¼ 0; where A¼ B ¼ FZ
tZ P
þ 2
FZ ðZ 2 H 2Z Þ
! bP ;
ðH 2Z þ Z 2 Þ2 ! tZ P ðZ 2 H 2Z Þ et2 H P P Z bP þ : 2 ðH P þ P Þ ðH 2Z þ Z 2 Þ2 ðH P þ P Þ3
ðH P þ P Þ
Applying Routh–Hurwitz criteria, E* is locally asymptotically stable provided A > 0 and B > 0 which gives tZ P ðH P þP Þ2
<
FZ ðH 2Z Z 2 Þ ðH 2Z þZ 2 Þ2
þ bP , where (P*, Z*) is the non-trivial equilibrium point.
h
5. Simulation results To better understand the spatial dynamics of the proposed diffusion model, we turn to its numerical simulation. The computer simulation of the model system (2a) and (2b) was performed using MATLAB 7.0. 5.1. The spatio-temporal dynamics in one-dimensional case In this subsection, we focus on the one-dimensional dynamics of the system (2a) and (2b) with u(x, t), v(x, t) and D = o2/ox2. Eqs. (2a) and (2b) takes the following form: ou vu o2 u þ 2; ¼ uð1 uÞ ot ðH P þ uÞ ox
ð5aÞ
ov vu v2 o2 v cv F þD 2: ¼b 2 2 ot ðH P þ uÞ ox ðH Z þ v Þ
ð5bÞ
In Eqs. (5a) and (5b), 0 < x < Lx, 0 < y < Ly. Model system (5a) and (5b) describes the dynamics of an aquatic system in a horizontal layer, vertical distribution of species inside the layer is assumed to be homogeneous. In our simulation, we
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
267
have considered a square domain of dimension Lx = Ly = 400. To be more realistic from the biological point of view, we consider a complex and non-monotic form of initial condition which determines the initial spatial distribution of the species in the real community as under uðx; 0Þ ¼ u þ eðx x1 Þðx x2 Þ; vðx; 0Þ ¼ v ;
ð6Þ
where (u*, v*) is the non-trivial state for the co-existence of prey and predator. e is the parameter which affects the system dynamics. At the domain boundary, we consider the zero flux boundary conditions: ou ou ov ov ¼ ¼ ¼ ¼ 0: ð7Þ oxð0;tÞ oxðL;tÞ oxð0;tÞ oxðL;tÞ This type of boundary condition is specifically used for modeling the dynamics of spatially bounded aquatic ecosystems [27]. We call a part of space a patch if it has a distinct kind of local spatial dynamics than adjacent parts. WOC propagates in both forward as well as in backward directions replacing a patch with local regular dynamics (characterized by a stable limit cycle in the phase plane of the system) with chaotic dynamics. These patches are shown in Fig. 2. These patches [9] are formed as a result of the action of diffusion on the interacting biological species (phytoplankton and zooplankton in the present case). With initial condition (6), when computer simulation is performed, a complex spatial dynamics is observed at t = 1000 as shown in Fig. 4a where a smooth regular patch separates chaotic patches at x1 = 1200 and x2 = 2800 with other parameters as in Fig. 3. We found a set of parameter values HP = 0.3, b = 2, c = 0.8, F = 1 · 104, HZ = 2.5, e = 1 · 108 where chaotic dynamics is observed with an appearance of jagged pattern at t = 600 as shown in Fig. 3c. Slightly changing the parameter values of b and F, keeping the remaining parameters same as above stable focus and limit cycle are obtained as in Figs. 3a and 3b respectively. At t = 600, numerical simulation for following set of parameter values HP = 0.3, b = 1.1, c = 0.8, F = 0.2, HZ = 2.5, u* = 0.6908, v* = 0.1566 gives a stable focus as shown in Fig. 3a. While taking b = 2, F = 1 · 103 and keeping remaining parameters same as above we get a limit cycle solution (cf. Fig. 3b). We obtained chaotic behavior for the parameter values HP = 0.3, b = 2, c = 0.8, F = 1 · 104, HZ = 2.5 and the corresponding (u*, v*) = (0.2, 0.4). Then after fixing these parameters, we tried to study the effect of time on the chaotic dynamics of the model system. We observed from our simulation results that the chaotic dynamics persists for a long period of time. The spatial distributions vary gradually in time after the local temporal behavior of the dynamical variables u and v follow a stable limit cycle of the homogeneous system (5a) and (5b) (cf. Fig. 4). The one-dimensional analysis shows the formation of strongly irregular transient pattern inside a subdomain. These jagged chaotic patterns and regular smooth patterns occupy the space alternatively. From Fig. 4, we observed that this jagged pattern representing chaotic behavior of the system grows steadily with time. The size of the domain occupied by the irregular chaotic patterns slowly grows with time in both directions equally, displacing the regular pattern and hence occupying the whole region. This phenomenon has been termed as ‘‘wave of chaos’’ (WOC) [18]. This phenomenon persists for the long period of time and after the irregular pattern occupies the whole domain, the dynamics of the system does not undergo any further changes. For the set of parameter values of Fig. 3c, the temporal behavior of prey and predator abundances are also chaotic as shown in Fig. 5. The phase plane is obtained at fixed points x1 = 1200 and x2 = 2800 of the domain. By changing the value of the parameter b in the above set of parameters, a limit cycle is obtained at b = 1.9 shown in Fig. 6 and a stable focus is observed at b = 1.2.
Fig. 2. Patchy structure of the system. Individual patches have distinct local dynamics in the beginning. Patches with regular spatial dynamics are constantly replaced by ones with chaotic dynamics. WOC generates chaos in the whole domain.
268
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276 0.9
1.4 prey predator
0.8
1.2
0.7
1
0.6
prey,predator
prey,predator
prey predator
0.5 0.4
0.8
0.6
0.4 0.3 0.2
0.2 0.1
0
0.2
0.4
0.6
0.8
1 space
1.2
1.4
1.6
1.8
2
0
0
0.2
0.4
0.6
0.8
4
x 10
1 space
1.2
1.4
1.6
1.8
2 4
x 10
1.4 prey predator
1.2
prey,predator
1
0.8
0.6
0.4
0.2
0
0
0.2
0.4
0.6
0.8
1 space
1.2
1.4
1.6
1.8
2 4
x 10
Fig. 3. Spatial distribution of the population showing the onset of chaotic phase after gradually displacing the regular phase at t = 600 of model (5a) and (5b) for parameter values HP = 0.3, c = 0.8, HZ = 2.5, Dt = 101, h = 1. (a) b = 1.1, F = 0.2, (b) b = 2, F = 0.001, (c) b = 2, F = 0.0001.
2D scan results on the parameter space (c, b) and (F, b) has been presented to observe simultaneous effect of these parameters on the dynamics of the model system (5a) and (5b) as shown in Fig. 7. From Fig. 7, it is found that the chaotic dynamics is observed at discrete points for b = 2 and in the range 0 6 F 6 0.06 (Fig. 7b). In the parameter space (c, b) chaos is observed along the line b = c. We performed extensive numerical computations to test whether spatio-temporal chaos generated through the action of WOC is robust with respect to changes in system parameters. It is evident from these figures that the chaotic dynamics is robust to changes for the range 0 6 F 6 0.06 at b = 2 in parameters space (F, b). 5.2. The spatio-temporal dynamics in two-dimensional case In this subsection, we consider dynamics of the two-dimensional model system (2a) and (2b) with u = u(x, y, t), v = v(x, y, t) and D = o2/ox2 + o2/oy2. Eqs. (2a) and (2b) takes the following form: ou vu o2 u o2 u þ 2þ 2; ¼ uð1 uÞ ot ðH P þ uÞ ox oy
2 ov vu v2 o v o2 v cv F ; þ þ D ¼b ot ðH P þ uÞ ox2 oy 2 ðH 2Z þ v2 Þ
ð8aÞ ð8bÞ
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276 1.4
1.4 prey predator
1.2
prey,predator
prey,predator
1
0.8
0.6
0.8
0.6
0.4
0.4
0.2
0.2
0
1000
2000
3000 4000 space
5000
6000
0
7000
0
prey predator
1.2
2000
3000 4000 space
5000
6000
7000
prey predator
1.2
1
prey,predator
1
0.8
0.6
0.8
0.6
0.4
0.4
0.2
0.2
0
1000
1.4
1.4
prey,predator
prey predator
1.2
1
0
269
0
1000
2000
3000 4000 space
5000
6000
0
7000
0
1000
2000
3000
4000 space
5000
6000
7000
1.4 prey predator
1.2
prey,predator
1 0.8 0.6 0.4 0.2 0
0
1000
2000
3000
4000 space
5000
6000
7000
Fig. 4. Simulation showing the long-term persistence of chaotic nature as the regular phase is gradually displaced by the chaotic phase (a) t = 1000 (b) t = 2000 (c) t = 5000 (d) t = 10,000 (e) t = 15,000.
The spatio-temporal dynamics of a chaotic system largely depends upon the choice of initial conditions [19]. We present the results of computer simulations with two different initial spatial distributions. First simulation result of two-dimensional case, the initial distribution of the species considered is a two-dimensional generalization of initial condition given in Eq. (6) (note that the initial conditions are deliberately chosen so that it is asymmetrical in order to make more influence of the corners of the domain). It is of the following form: uðx; y; 0Þ ¼ u e1 ðx 0:1y 225Þðx 0:1y 675Þ; vðx; y; 0Þ ¼ v e2 ðx 450Þ e3 ðy 150Þ:
ð9Þ
270
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276 1.4
1.2
Predator
1
0.8
0.6
0.4
0.2
0
0
0.1
0.2
0.3
0.4 Prey
0.5
0.6
0.7
0.8
Fig. 5. Phase plane of the system for the spatial range defined by fixed points x1 = 1200, x2 = 2800 inside the domain occupied by irregular spatio-temporal oscillation. Parameters are given in the text. 1 0.9 0.8 0.7
Predator
0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.1
0.2
0.3
0.4 Prey
0.5
0.6
0.7
0.8
Fig. 6. Phase plane of the system showing limit cycle behavior with b = 1.9. Other parameters are same as in Fig. 4.
3.5
2.9
3
2.6 2.3
2.5
2
2
1.7
b 1.5
b 1.4
1
1.1
0.5
0.8
0
0.5
0
0.5
1
c
1.5
2
2.5
0
0.02
0.04
0.06
0.08
0.1
F
Fig. 7. 2D scan diagrams for the model system (5a) and (5b) in the parameter space (c, b) and (F, b).
The numerical results obtained using semi-implicit numerical techniques of two-dimensional phytoplankton–zooplankton system are presented in Fig. 8. The parameter values at which intermittent chaos was obtained in 1D analysis and initial condition (9) with zero-flux boundary conditions (as explained in 1D analysis) in two dimension are used to perform simulations of two-dimensional spatio-temporal system. The parameter values taken are, HP = 0.3, b = 2, c = 0.8, HZ = 2.5, F = 1 · 104, e1 = 2 · 107, e2 = 3 · 105, e3 = 1.2 · 104 Dt = 1/3, h = 1, d = 1.
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
271
Phytoplankton population density is plotted in first column and zooplankton population density in the second column of Fig. 8 at t = 0, 150, 200, 300, 600 and 1000. Fig. 8 clearly depicts the variations in the population behavior of phytoplankton and zooplankton at the early stages i.e. at t = 0, 150, 200. But at later stages, i.e. at t = 300, 600, 1000
Fig. 8. Spatial distribution of prey densities (phytoplankton) [first column] and predator densities (zooplankton) [second column] of model system (8a) and (8b) with initial condition given in Eq. (9) for (a) t = 0, (b) t = 150, (c) t = 200, (d) t = 300, (e) t = 600, (f) t = 1000. Parameter values are given in the text.
272
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
Fig 8. (continued)
both prey and predator densities exhibit almost similar spatial population dynamics. The regular spiral spatial pattern is followed by the irregular patchy patterns. Absence of spiral patterns at t = 0 justifies that the appearance of spirals is not induced by initial condition. These spiral patterns first grow till a particular stage after which its destruction begins from their centers. Finally an irregular pattern prevails over the entire square domain at t = 600 which persists with time as seen at t = 1000.
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
273
Fig. 9. Spatial distribution of prey densities (phytoplankton) [first column] and predator densities (zooplankton) [second column] using initial condition given in Eq. (10) for (a) t = 0, (b) t = 150, (c) t = 200, (d) t = 300, (e) t = 600, (f) t = 1000. Parameter values are given in the text.
Although choice of the initial condition drastically affects the dynamics of the system but the spatial dynamics observed in Fig. 8 is not induced by initial condition. To justify this fact, we consider another initial condition for the second set of simulation, which is
274
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
Fig 9. (continued)
uðx; y; 0Þ ¼ u e1 ðx 180Þðx 720Þ e2 ðy 90Þðy 210Þ; vðx; y; 0Þ ¼ v e3 ðx 450Þ e4 ðy 135Þ and the parameters are the same as those used to generate Fig. 8.
ð10Þ
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
275
Similar scenario can be seen in Fig. 9. An irregular patchy pattern is finally obtained preceded by a somewhat less regular spiral pattern. We observed that the dynamics of the system is regular only for a specified period of time after which an irregular dynamics prevails. Also the chaotic nature of the minimal diffusion model has been studied and it is observed that spatio-temporal chaotic dynamics of this interaction–diffusion system is self-organized.
6. Discussions and conclusion These spatial patterns can be considered as time-independent as after certain amount of time is elapsed, there were not visible changes in spatial patterns (cf. Figs. 8e and 8f, 9e and 9f). The spatio-temporal chaos generated by propagation of WOC is robust to changes in key system parameters such as intensity of competition between individuals of phytoplankton species and changes in the per capita rate of predation by fishes. Two-dimensional spatial patterns of phytoplankton–zooplankton dynamics are self-organized and, therefore, can be considered to provide a theoretical framework to understand patchiness in marine environments. Recently, Medvinsky and collaborators [28] have proposed a new mechanism for the origin of spatio-temporal chaos in a heterogeneous environment wherein nearby habitats possess different kind of local dynamics. It remains to be seen which of these mechanisms is more common in marine and other aquatic systems. The discovery of spatio-temporal chaos [29] in early nineties has potential to lead to many noble applications in Engineering [30,31] and Physiology [32,33]. There exist many mechanisms, which cause spatio-temporal chaos in a given system. It is imperative to assess their potential to solve problems in Physics, Biology and Engineering. Here we attempt to do the same in a model ecological system.
Acknowledgement This research work is supported by the All India Council for Technical Education which has offered the award of AICTE National Doctoral Fellowship to Nitu Kumari (F.No.:1-10/FD/NDF-PG/(29)/2006-07).
References [1] Rai V, Schaffer WM. Chaos in ecology. Chaos, Solitons & Fractals 2001;12:197–203. [2] Franks PJS. Spatial patterns in dense algal blooms. Limnol Oceanol 1997;42:1297–305. [3] Legendre L. The significance of microalgal blooms for fisheries and for the export of particulate organic carbon in Oceans. J Plank Res 1990;12:681–99. [4] Hillary RM, Bees MA. Plankton lattices and the role of chaos in plankton patchiness. Phys Rev E 2004;69:031913. [5] Abraham ER. The generation of plankton patchiness by turbulent stirring. Nature 1998;391:577–80. [6] Turing AM. On the chemical basis of morphogenesis. Philos Trans R Soc London Ser B 1952;237:37–72. [7] Segel LA, Jackson JL. Dissipative structures: an explanation and an ecological example. J Theor Biol 1972;37:545–59. [8] Levin SA, Segel LA. Hypothesis for origin of planktonic patchiness. Nature 1976;259:659. [9] Medvinsky AB, Tikhonova IA, Aliev RR, Li BL, Lin ZS, Malchow H. Patchy environment as a factor of complex plankton dynamics. Phys Rev E 2001;64:021915–7. [10] Fasham MJR. The statistical and mathematical analysis of plankton patchiness. Oceanogr Mar Biol Ann Rev 1978;16:43–79. [11] Franks PJS, Wroblewski JS, Fileri GR. Behavior of simple plankton model with food-level acclimation by herbivores. Mar Biol 1986;91:121–9. [12] Steel JH, Hunderson EW. A simple model for plankton patchiness. J Plankton Res 1992;14:1397–403. [13] Truscott JE, Brindley J. Equilibria, stability and excitability in a general class of plankton population models. Philos Trans R Soc London Ser A 1994;347:703–18. [14] Scheffer M. Fish and nutrients interplay determines algal biomass: a minimal model. OIKOS 1991;62:271–82. [15] Malchow H. Spatio-temporal pattern formation in non-linear nonequilibrium plankton dynamics. Proc R Soc London Ser B 1993;251:103–9. [16] Malchow H. Nonequilibrium structures in plankton dynamics. Ecol Model 1994;75–76:123–34. [17] Pascual M. Diffusion induced chaos in a spatial predator–prey system. Proc R Soc London Ser B 1993;103:251–7. [18] Petrovskii SV, Malchow H. Wave of chaos: new mechanism of pattern formation in spatio-temporal population dynamics. Theor Populat Biol 2001;59:157–74. [19] Medvinsky A, Petrovskii S, Tikhonova I, Malchow H, Li Bai-Lian. Spatiotemporal complexity of plankton and fish dynamics. SIAM Rev 2002;44:311–70. [20] Malchow H, Petrovskii SV, Medvinsky AB. Numerical study of plankton–fish dynamics in a spatially structured and noisy environment. Ecol Model 2002;149:247–55.
276
R.K. Upadhyay et al. / Chaos, Solitons and Fractals 40 (2009) 262–276
[21] Malchow H. Motional instabilities in prey–predator systems. J Theor Biol 2000;204:639–47. [22] Malchow H, Hilker F, Petrovskii S. Noise and productivity dependence of spatiotemporal pattern formation in a prey–predator system. Discr Cont Dyn – B 2004;4:705–11. [23] Malchow H, Radtke B, Kallache M, Medvinsky A, Tikhonov D, Petrovskii S. Spatio-temporal pattern formation in coupled models of plankton dynamics and fish school motion. Non-linear Anal – Real World Appl 2000;1:53–67. [24] Tikhonova I, Li B, Malchow H, Medvinsky A. The impact of the phytoplankton growth rule on spatial and temporal dynamics of plankton communities in a heterogeneous environment. Biofizika 2003;48:891–9. [25] Sherratt J, Eagan B, Lewis M. Oscillations and chaos behind predator–prey invasion: mathematical artifact or ecological reality? Philos Trans R Soc London B 1997;352:21–38. [26] Malchow H, Petrovskii S, Medvinsky AB. Pattern formation in models of plankton dynamics: a synthesis. Oceanol Acta 2001;24:479–87. [27] Scheffer M. Population and community biology series 22. Ecology of shallow lakes. London: Chapman & Hall; 1998. [28] Medvinsky AB. Personal communication, 2006. [29] Cross MC, Hohenberg PC. Spatio-temporal chaos. Science 1994;263:1569–70. [30] Hu G, Xiao J, Yang J, Xie F, Qu Z. Synchronization of spatio-temporal chaos and its applications. Phys Rev E 1997;56(3):2738–46. [31] Gang H, Zhilin Q. Controlling spatiotemporal chaos in coupled map lattices. Phys Rev Lett 1994;72:68–71. [32] Zhang H, Patel N. Spiral Wave breakdown in an excitable medium model of cardiac tissue. Chaos, Solitons & Fractals 1995;5(3– ):635–9. [33] Bramann Y, John F, Lindner JF, Ditto WL. Taming spatio-temporal chaos with disorder. Nature 1995;378:465–7. [34] Tikhonov DA, Enderlein J, Malchow H, Medvinsky AB. Chaos and fractals in fish school motion. Chaos, Solitons & Fractals 2001;12:277–88.