- Email: [email protected]
Spatial organization in a plant model
Applied Mathematics and Computation 216 (2010) 983–988
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Spatial organization in a plant model Zhao-Hua Liu, Ai-Ling Wang, Biao Wang, Yong-Jiang Liu * Key Laboratory for AMT Shanxi, North University of China, Taiyuan 030051, People’s Republic of China
a r t i c l e
i n f o
Keywords: Plant Diffusion Spatial organization
a b s t r a c t In this paper, a spatial plant–wrack model is presented. We find such model has spotted pattern dynamics by using both mathematical analysis and numerical simulations. The obtained results may point that the distribution of the plant in space is self-organized. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction Vegetation pattern, which is a subfield of spatial ecology [3,13,1], has been extensively studied by arid land ecologists [14]. The underlying forces in arid lands are water scarcity, plant competition over water resources, and redistribution of water by runoff. The view of the pattern formation phenomenon involving symmetry breaking, which is supported by mathematical models that identify vegetation patterns with instabilities of uniform vegetation states is fairly new [6,5,9,10]. The models account for field observations of various vegetation patterns including bands on hill slopes and spots in plain areas [4]. In [8], Koppel et al presented an experimental investigation of regular spatial patterning in Carex stricta and used this ecosystem as a model to investigate alternate theoretical mechanisms that drive regular spatial patterning in general. Carex stricta, the tussock sedge, is a species with widespread distribution in freshwater marshes of North America. Tussock formation in wetland plants is a fairly common occurrence and appears to be a stress avoidance strategy for plants that are able to escape waterlogging and associated low oxygen levels and soil toxicity by elevating their rooting substrate. In this paper, we investigate a model given by Koppel et al. [8], and investigate its potential for pattern formation analytically. The main objective is to determine how patterning varies with model parameters. The paper is organized as follows. In Section 2, we give the model.And we analyse the spatial model in Section 3. With respect to these parameters, we locate the Turing domain within the parameter space. In Section 3, by performing a series of simulations, we illustrate the emergence of Turing patterns. Finally, some conclusions are given.
2. Main model The model, which describes the interaction of the plant and wrack is as follows [8]:
dP ¼ Pð1 PÞFðPÞ sP GðP; WÞ; dt dW ¼ sP bW; dt
* Corresponding author. E-mail address: [email protected] (Y.-J. Liu). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.01.116
ð1aÞ ð1bÞ
984
Z.-H. Liu et al. / Applied Mathematics and Computation 216 (2010) 983–988
where P is plant biomass, W is wrack biomass, F(P) is a function describing the positive effect of plant biomass on its own growth, s is the specific rate of plant senescence, G(P, W) is a function describing the inhibiting effect of wrack on plant growth as a function of plant and wrack biomass, b is the decay rate of wrack. In this model, we consider the implications of facilitation of plant growth by the formation of a root mound, modeled as a positive effect of vegetation on its own growth, which leads to that:
dP P ¼ Pð1 PÞ sP aPW; dt PþK dW ¼ sP bW; dt
ð2aÞ ð2bÞ
where the term P=ðP þ KÞ models the increase of growth due to facilitation by the root mound, assuming that the size of the mound increases linearly with the amount of plant biomass. Hence, it assumes that without a root mound, growth is close to 0 and the per capita growth increases to 1 as plant growth increases. The quantity K represents the plant biomass where the facilitation term is half maximal. Next, we intend to introduce the spatial model. Fick’s law holds for the reason that the dispersal of individuals can be taken as random, which results in flux terms given as
@P ¼ D1 r2 P; @t
@W ¼ D2 r2 W; @t
ð3Þ
where r2 ¼ @ 2 =@X 2 þ @ 2 =@Y 2 is the usual Laplacian operator in two-dimensional space. Here, X and Y represent space. The diffusion coefficients are denoted as D1 and D1 , respectively. Combining the spatial term, we have the following model:
@P P ¼ Pð1 PÞ sP aPW þ D1 r2 P; dt PþK @W ¼ sP bW þ D2 r2 W: dt
ð4aÞ ð4bÞ
3. Bifurcation analysis We first consider a spatially homogeneous system and find the steady state as follows: (i) E1 ¼ ð0; 0Þ, which is corresponding to extinction of the plant and wrack; (ii) positive equilibrium point E1 ¼ ðP 1 ; W 1 Þ, which is corresponding to coexistence of plant and wrack and 1
P ¼
b asK sb
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 b 2basK 2sb þ a2 s2 K 2 2as2 Kb þ s2 b 4sKb 2ðb þ asÞ
ð5Þ
and
W1 ¼
sP 1 : b
ð6Þ
(iii) positive equilibrium point E ¼ ðP ; W Þ, which is corresponding to coexistence of plant and wrack and
P ¼
b asK sb þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 b 2basK 2sb þ a2 s2 K 2 2as2 Kb þ s2 b 4sKb 2ðb þ asÞ
ð7Þ
and
W ¼
sP : b
ð8Þ 1
By linear stability analysis, we found that E is an unstable saddle. To consider pattern formation for system (4), we need to look for the dispersion relation. Following [12], we will give the linear stability of system (4) following the standard route [2,7]. We make the following substitute:
P ¼ P þ Pð~ r; tÞ
ð9Þ
and
W ¼ W þ Wð~ r; tÞ;
ð10Þ
into the Eq. (4) and assume jPj; jWj 1. Here, ~ r ¼ X (Y) or ~ r ¼ ðX; YÞ, which is corresponding to the one- or two-dimension space. Then, in the linear approximation, we have
@P ¼ a11 P þ a12 W þ D1 r2 P; @t @W ¼ a21 P þ a22 W þ D2 r2 W: @t
ð11aÞ ð11bÞ
Z.-H. Liu et al. / Applied Mathematics and Computation 216 (2010) 983–988
985
The initial conditions are assumed as that
Pjt¼0 ¼ f ð~ rÞ
ð12Þ
Wjt¼0 ¼ gð~ rÞ;
ð13Þ
and
where the functions f ð~ rÞ and gð~ rÞ decay rapidly for ~ r ! 1. Following the standard approach, let us now perform a Laplace transformation of the linearized equations over the two independent variables ~ r and t. For ~ r we use the so-called two-sided version of the transformation. The relations for the forward and backward transforms are
Pkq ¼
Z
1
ekt dt
Z
0
þ1
Pð~ r; tÞeq~r d~ r
ð14Þ
1
and
Pð~ r; tÞ ¼
1 4p2
Z
bþi1
ekt ds
Z
i1
Pkq eq~r dq;
ð15Þ
i1
bi1
where k and q are complex variables. And s is the Laplace transform variable, q is the Fourier transform variable. That is to say that, q ¼ ik or (ik, il) corresponding to one- and two-dimensional space, and the wave numbers k and l are real number. In formula (10) for the backward transformation, the integration contour in the q-plane is the imaginary axis. In the s-plane the contour is parallel to the imaginary axis and located to the right of all singularities of the integrand. After this transformation, the kinetic equations read
ðk a11 D1 q2 ÞP kq a12 W kq ¼ HðqÞ
ð16Þ
ðk a22 D2 q2 ÞW kq a21 Pkq ¼ ZðqÞ;
ð17Þ
and
where H(q) and Z(q) are the transforms of f ð~ rÞ and gð~ rÞ. To reveal the presence of an instability and disclose its character, it is sufficient to consider one variable. The temporal growth of the perturbations can now be found by inverting the Laplace transforms, which follows directly after factorizing the denominator. By solving the linear Eqs. (12) and (13) we find P kq and then use the backward transformation (11) to obtain the following formal solution:
Pð~ r; tÞ ¼
1 4p2
Z
bþi1
ekt ds
Z
bi1
i1
i1
ðk a22 d2 q2 ÞHðqÞ þ a12 ZðqÞ q~r e dq: Dðk; qÞ
ð18Þ
Then we obtain the denominator
Dðk; qÞ ¼ ðk a22 d2 q2 Þðk a11 d1 q2 Þ a12 a21 :
ð19Þ
Hence the dispersion equation is that
ðk a22 d2 q2 Þðk a11 d1 q2 Þ a12 a21 ¼ 0:
ð20Þ
A symmetry breaking occurs when a homogeneous steady state solution of the system (4) is linearly stable to perturbations in the absence of the diffusion and advection terms but linearly unstable to small spatial perturbations in the presence of diffusion and advection. The condition for instability is that (16) has a root s with positive real part. We find the roots from that
k¼
aðq2 Þ þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ½aðq2 Þ2 4bðq2 Þ 2
;
ð21Þ
where
aðq2 Þ ¼ ðd1 þ d2 Þq2 þ ða11 þ a22 Þ; 2
4
2
bðq Þ ¼ d1 d2 q þ ða11 þ a22 Þq þ a11 a22 a12 a21 :
ð22aÞ ð22bÞ
The condition for a spatial mode q (in one- or two-dimensional space) to be unstable and thus grow into a pattern is that ReðkÞ > 0. To well see the dispersion relation, we set that s ¼ 0:15; b ¼ 0:1; a ¼ 1; D1 ¼ 0:01; D2 ¼ 1 and vary K. Fig. 1 depicts the range of the values of q for some constant values, and the small perturbation may bring about an instability with time. From Fig. 1, we can see that the spatial pattern can occur due to the positive real parts of k.
986
Z.-H. Liu et al. / Applied Mathematics and Computation 216 (2010) 983–988
Fig. 1. An illustration of the dispersion relation from the Eq. (17) and the parameter values are in the text. Green line: K ¼ 0:14; Red line: K ¼ 0:16. (For interpretation of colour in figures, the reader is referred to the Web version of this article.)
Fig. 2. Snapshots of contour pictures of the time evolution of the wrack at different s ¼ 0:15; b ¼ 0:1; a ¼ 1; D1 ¼ 0:01; D2 ¼ 1 and K ¼ 0:14. (A) t ¼ 0; (B) t ¼ 100; (C) t ¼ 350; (D) t ¼ 550.
times.
Parameter
values
are
used
as:
4. Main results The dynamics behavior of the spatial model can not be studied by using analytical methods To solve differential equations using computers, one has to discretize the space and time of the problem; that is to say, transform it from an infinite-dimensional (continuous) to a finite-dimensional (discrete) form. In practice the continuous problem defined by the reaction-dif-
Z.-H. Liu et al. / Applied Mathematics and Computation 216 (2010) 983–988
987
fusion system in two-dimensional space is solved in a discrete domain with M N lattice sites. The spacing between the lattice points is defined by the lattice constant Dh. In the discrete system the Laplacian describing diffusion is calculated using finite differences, i.e., the derivatives are approximated by differences over Dh. For Dh ! 0 the differences approach the derivatives. The time evolution is also discrete, i.e., the time goes in steps of Dt. The time evolution can be solved by using the Euler method, which means approximating the value of the concentration at the next time step on the basis of the change rate of the concentration at the previous time step [11]. In the following, we will perform a series of numerical simulations of the spatially extended model (4) in two-dimensional spaces, and the qualitative results are shown in figures. All our numerical simulations employ zero-flux boundary conditions with a system size of 50 50 space units. The model (4) is solved numerically in two-dimensional space using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration with a time step size of Dt ¼ 0:05 and space step size (lattice constant) Dh ¼ 1. The scale of the space and time is average for the Euler method. The initial density distribution corresponds to random perturbations around the stationary state ðP ; W Þ in model (4) with a variance significantly lower than the amplitude of the final patterns, which seems to be more general from the biological point of view. In the simulations, we have found that the distributions of plant and wrack are always of the same type. As a result, we can restrict our analysis of pattern formation to one distribution (in this paper, we show the distribution of wrack, for instance). In Fig. 2, we show the evolution of the spatial pattern of wrack at t ¼ 0; 100; 350 and 550, with random perturbation of the steady state P and W . The random initial distribution leads to the formation of a strongly irregular transient pattern in the domain. At the final (cf. Fig. 2(C)–(D)), regular spots prevail over the whole domain and the dynamics of the system does not undergo any further changes. The temporal development of the spatially averaged values of the wrack and plant density shows that in the first intervals of simulations these values change fast the variance increases (cf. Fig. 3(a)–(b)). One can see that at t ¼ 250, the variance reaches a constant values and the mean increases slower. At this time, the pattern formation has almost completely evolved. And the system has reached its steady state at t > 350.
Mean value of the wrack density
0.334 0.332 0.33 0.328 0.326 0.324 0.322 0.32 0.318 0
100
200
300 Time
400
500
Mean value of the plant dentisy
0.222
0.22
0.218
0.216
0.214
0.212 0
100
200
300 Time
400
500
Fig. 3. The temporal development of the spatial average and variance of wrack and plant density plotted against time. (a) The time series of the mean value of the wrack; (b) The time series of the mean value of the plant. Parameter values are used as the same in Fig. 3.
988
Z.-H. Liu et al. / Applied Mathematics and Computation 216 (2010) 983–988
5. Discussion and conclusion In this paper, a spatial plant–wrack model with spatial diffusion and is investigated by mathematical analysis and numerical simulations. We found that the model exhibits spotted pattern, which may indicate that the plant and wrack are selforganized in the space. A question naturally arising is whether the number of the grid, used in numerical simulations, has great effect on the spatial pattern of the plant and wrack. To address this issue, we performed extensive computer simulations for different parameter values of the number of the grid (in total, more than 10 parameter sets were examined). We found that the behaviors are very similar (for the sake of brevity, we do not provide separated figures). Acknowledgments Supported by the PhD Research Fund of the Ministry of Education of China under Grant No 20050204JJ, and Scientific and Technological Project of the Shanxi Province (20080209ZX). References [1] M.R. Aguiar, O.E. Sala, Plant structure, dynamics and implications for the functioning of arid ecosystems, Trends Ecol. Evol. 14 (1999) 273–277. [2] P. Andresen, M. Bache, E. Mosekilde, G. Dewel, P. Borckmanns, Stationary space-periodic structures with equal diffusion coefficients, Phys. Rev. E 60 (1999) 297–301. [3] J. Bascompte, R.V. Sole, Modeling Spatiotemporal Dynamics in Ecology, Springer-Verlag, Berlin, 1998. [4] J. von Hardenberg, E. Meron, M. Shachak, Y. Zarmi, Diversity of vegetation patterns and desertification, Phys. Rev. Lett. 87 (2001) 198101. [5] R. HilleRisLambers, M. Rietkerk, F. van den Bosch, H.H.T. Prins, H. de Kroon, Vegetation pattern formation in semi-arid grazing systems, Ecology 82 (2001) 50–61. [6] C.A. Klausmeier, Regular and irregular patterns in semiarid vegetation, Science 284 (1999) 1826–1828. [7] S.P. Kuznetsov, E. Mosekilde, G. Dewel, P. Borckmans, Absolute and convective instabilities in a one-dimensional brusselator flow model, J. Chem. Phys. 106 (1997) 7609–7616. [8] J. van de Koppel, C.M. Crain, Scale-dependent inhibition drives regular tussock spacing in a freshwater marsh, Amer. Nat. 168 (2006) E136–E147. [9] R. Lefever, O. Lejeune, On the origin of tiger bush, Bull. Math. Biol. 59 (1997) 263–294. [10] R. Lefever, O. Lejeune, P. Couteron, in: P.K. Maini, H. Othmer (Eds.), Mathematical Models for Biological Pattern Formation, IMA, Springer, New York, 2000. [11] H. Malchow, S.V. Petrovskii, E. Venturino, Spatiotemporal Patterns in Ecology and Epidemiology: Theory, Models, and Simulations, Chapman and Hall, CRC, London, 2008. [12] J.D. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications, Springer, New York, 2003. [13] D. Tilman, P. Kareiva, Spatial Ecology, Princeton University Press, Princeton, 1997. edited book. [14] C. Valentin, J.M. dHerbes, J. Poesen, Soil and water components of banded vegetation patterns, Catena 37 (1999) 1–24.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Spatial organization in a plant model Zhao-Hua Liu, Ai-Ling Wang, Biao Wang, Yong-Jiang Liu * Key Laboratory for AMT Shanxi, North University of China, Taiyuan 030051, People’s Republic of China
a r t i c l e
i n f o
Keywords: Plant Diffusion Spatial organization
a b s t r a c t In this paper, a spatial plant–wrack model is presented. We find such model has spotted pattern dynamics by using both mathematical analysis and numerical simulations. The obtained results may point that the distribution of the plant in space is self-organized. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction Vegetation pattern, which is a subfield of spatial ecology [3,13,1], has been extensively studied by arid land ecologists [14]. The underlying forces in arid lands are water scarcity, plant competition over water resources, and redistribution of water by runoff. The view of the pattern formation phenomenon involving symmetry breaking, which is supported by mathematical models that identify vegetation patterns with instabilities of uniform vegetation states is fairly new [6,5,9,10]. The models account for field observations of various vegetation patterns including bands on hill slopes and spots in plain areas [4]. In [8], Koppel et al presented an experimental investigation of regular spatial patterning in Carex stricta and used this ecosystem as a model to investigate alternate theoretical mechanisms that drive regular spatial patterning in general. Carex stricta, the tussock sedge, is a species with widespread distribution in freshwater marshes of North America. Tussock formation in wetland plants is a fairly common occurrence and appears to be a stress avoidance strategy for plants that are able to escape waterlogging and associated low oxygen levels and soil toxicity by elevating their rooting substrate. In this paper, we investigate a model given by Koppel et al. [8], and investigate its potential for pattern formation analytically. The main objective is to determine how patterning varies with model parameters. The paper is organized as follows. In Section 2, we give the model.And we analyse the spatial model in Section 3. With respect to these parameters, we locate the Turing domain within the parameter space. In Section 3, by performing a series of simulations, we illustrate the emergence of Turing patterns. Finally, some conclusions are given.
2. Main model The model, which describes the interaction of the plant and wrack is as follows [8]:
dP ¼ Pð1 PÞFðPÞ sP GðP; WÞ; dt dW ¼ sP bW; dt
* Corresponding author. E-mail address: [email protected] (Y.-J. Liu). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.01.116
ð1aÞ ð1bÞ
984
Z.-H. Liu et al. / Applied Mathematics and Computation 216 (2010) 983–988
where P is plant biomass, W is wrack biomass, F(P) is a function describing the positive effect of plant biomass on its own growth, s is the specific rate of plant senescence, G(P, W) is a function describing the inhibiting effect of wrack on plant growth as a function of plant and wrack biomass, b is the decay rate of wrack. In this model, we consider the implications of facilitation of plant growth by the formation of a root mound, modeled as a positive effect of vegetation on its own growth, which leads to that:
dP P ¼ Pð1 PÞ sP aPW; dt PþK dW ¼ sP bW; dt
ð2aÞ ð2bÞ
where the term P=ðP þ KÞ models the increase of growth due to facilitation by the root mound, assuming that the size of the mound increases linearly with the amount of plant biomass. Hence, it assumes that without a root mound, growth is close to 0 and the per capita growth increases to 1 as plant growth increases. The quantity K represents the plant biomass where the facilitation term is half maximal. Next, we intend to introduce the spatial model. Fick’s law holds for the reason that the dispersal of individuals can be taken as random, which results in flux terms given as
@P ¼ D1 r2 P; @t
@W ¼ D2 r2 W; @t
ð3Þ
where r2 ¼ @ 2 =@X 2 þ @ 2 =@Y 2 is the usual Laplacian operator in two-dimensional space. Here, X and Y represent space. The diffusion coefficients are denoted as D1 and D1 , respectively. Combining the spatial term, we have the following model:
@P P ¼ Pð1 PÞ sP aPW þ D1 r2 P; dt PþK @W ¼ sP bW þ D2 r2 W: dt
ð4aÞ ð4bÞ
3. Bifurcation analysis We first consider a spatially homogeneous system and find the steady state as follows: (i) E1 ¼ ð0; 0Þ, which is corresponding to extinction of the plant and wrack; (ii) positive equilibrium point E1 ¼ ðP 1 ; W 1 Þ, which is corresponding to coexistence of plant and wrack and 1
P ¼
b asK sb
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 b 2basK 2sb þ a2 s2 K 2 2as2 Kb þ s2 b 4sKb 2ðb þ asÞ
ð5Þ
and
W1 ¼
sP 1 : b
ð6Þ
(iii) positive equilibrium point E ¼ ðP ; W Þ, which is corresponding to coexistence of plant and wrack and
P ¼
b asK sb þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 b 2basK 2sb þ a2 s2 K 2 2as2 Kb þ s2 b 4sKb 2ðb þ asÞ
ð7Þ
and
W ¼
sP : b
ð8Þ 1
By linear stability analysis, we found that E is an unstable saddle. To consider pattern formation for system (4), we need to look for the dispersion relation. Following [12], we will give the linear stability of system (4) following the standard route [2,7]. We make the following substitute:
P ¼ P þ Pð~ r; tÞ
ð9Þ
and
W ¼ W þ Wð~ r; tÞ;
ð10Þ
into the Eq. (4) and assume jPj; jWj 1. Here, ~ r ¼ X (Y) or ~ r ¼ ðX; YÞ, which is corresponding to the one- or two-dimension space. Then, in the linear approximation, we have
@P ¼ a11 P þ a12 W þ D1 r2 P; @t @W ¼ a21 P þ a22 W þ D2 r2 W: @t
ð11aÞ ð11bÞ
Z.-H. Liu et al. / Applied Mathematics and Computation 216 (2010) 983–988
985
The initial conditions are assumed as that
Pjt¼0 ¼ f ð~ rÞ
ð12Þ
Wjt¼0 ¼ gð~ rÞ;
ð13Þ
and
where the functions f ð~ rÞ and gð~ rÞ decay rapidly for ~ r ! 1. Following the standard approach, let us now perform a Laplace transformation of the linearized equations over the two independent variables ~ r and t. For ~ r we use the so-called two-sided version of the transformation. The relations for the forward and backward transforms are
Pkq ¼
Z
1
ekt dt
Z
0
þ1
Pð~ r; tÞeq~r d~ r
ð14Þ
1
and
Pð~ r; tÞ ¼
1 4p2
Z
bþi1
ekt ds
Z
i1
Pkq eq~r dq;
ð15Þ
i1
bi1
where k and q are complex variables. And s is the Laplace transform variable, q is the Fourier transform variable. That is to say that, q ¼ ik or (ik, il) corresponding to one- and two-dimensional space, and the wave numbers k and l are real number. In formula (10) for the backward transformation, the integration contour in the q-plane is the imaginary axis. In the s-plane the contour is parallel to the imaginary axis and located to the right of all singularities of the integrand. After this transformation, the kinetic equations read
ðk a11 D1 q2 ÞP kq a12 W kq ¼ HðqÞ
ð16Þ
ðk a22 D2 q2 ÞW kq a21 Pkq ¼ ZðqÞ;
ð17Þ
and
where H(q) and Z(q) are the transforms of f ð~ rÞ and gð~ rÞ. To reveal the presence of an instability and disclose its character, it is sufficient to consider one variable. The temporal growth of the perturbations can now be found by inverting the Laplace transforms, which follows directly after factorizing the denominator. By solving the linear Eqs. (12) and (13) we find P kq and then use the backward transformation (11) to obtain the following formal solution:
Pð~ r; tÞ ¼
1 4p2
Z
bþi1
ekt ds
Z
bi1
i1
i1
ðk a22 d2 q2 ÞHðqÞ þ a12 ZðqÞ q~r e dq: Dðk; qÞ
ð18Þ
Then we obtain the denominator
Dðk; qÞ ¼ ðk a22 d2 q2 Þðk a11 d1 q2 Þ a12 a21 :
ð19Þ
Hence the dispersion equation is that
ðk a22 d2 q2 Þðk a11 d1 q2 Þ a12 a21 ¼ 0:
ð20Þ
A symmetry breaking occurs when a homogeneous steady state solution of the system (4) is linearly stable to perturbations in the absence of the diffusion and advection terms but linearly unstable to small spatial perturbations in the presence of diffusion and advection. The condition for instability is that (16) has a root s with positive real part. We find the roots from that
k¼
aðq2 Þ þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ½aðq2 Þ2 4bðq2 Þ 2
;
ð21Þ
where
aðq2 Þ ¼ ðd1 þ d2 Þq2 þ ða11 þ a22 Þ; 2
4
2
bðq Þ ¼ d1 d2 q þ ða11 þ a22 Þq þ a11 a22 a12 a21 :
ð22aÞ ð22bÞ
The condition for a spatial mode q (in one- or two-dimensional space) to be unstable and thus grow into a pattern is that ReðkÞ > 0. To well see the dispersion relation, we set that s ¼ 0:15; b ¼ 0:1; a ¼ 1; D1 ¼ 0:01; D2 ¼ 1 and vary K. Fig. 1 depicts the range of the values of q for some constant values, and the small perturbation may bring about an instability with time. From Fig. 1, we can see that the spatial pattern can occur due to the positive real parts of k.
986
Z.-H. Liu et al. / Applied Mathematics and Computation 216 (2010) 983–988
Fig. 1. An illustration of the dispersion relation from the Eq. (17) and the parameter values are in the text. Green line: K ¼ 0:14; Red line: K ¼ 0:16. (For interpretation of colour in figures, the reader is referred to the Web version of this article.)
Fig. 2. Snapshots of contour pictures of the time evolution of the wrack at different s ¼ 0:15; b ¼ 0:1; a ¼ 1; D1 ¼ 0:01; D2 ¼ 1 and K ¼ 0:14. (A) t ¼ 0; (B) t ¼ 100; (C) t ¼ 350; (D) t ¼ 550.
times.
Parameter
values
are
used
as:
4. Main results The dynamics behavior of the spatial model can not be studied by using analytical methods To solve differential equations using computers, one has to discretize the space and time of the problem; that is to say, transform it from an infinite-dimensional (continuous) to a finite-dimensional (discrete) form. In practice the continuous problem defined by the reaction-dif-
Z.-H. Liu et al. / Applied Mathematics and Computation 216 (2010) 983–988
987
fusion system in two-dimensional space is solved in a discrete domain with M N lattice sites. The spacing between the lattice points is defined by the lattice constant Dh. In the discrete system the Laplacian describing diffusion is calculated using finite differences, i.e., the derivatives are approximated by differences over Dh. For Dh ! 0 the differences approach the derivatives. The time evolution is also discrete, i.e., the time goes in steps of Dt. The time evolution can be solved by using the Euler method, which means approximating the value of the concentration at the next time step on the basis of the change rate of the concentration at the previous time step [11]. In the following, we will perform a series of numerical simulations of the spatially extended model (4) in two-dimensional spaces, and the qualitative results are shown in figures. All our numerical simulations employ zero-flux boundary conditions with a system size of 50 50 space units. The model (4) is solved numerically in two-dimensional space using a finite difference approximation for the spatial derivatives and an explicit Euler method for the time integration with a time step size of Dt ¼ 0:05 and space step size (lattice constant) Dh ¼ 1. The scale of the space and time is average for the Euler method. The initial density distribution corresponds to random perturbations around the stationary state ðP ; W Þ in model (4) with a variance significantly lower than the amplitude of the final patterns, which seems to be more general from the biological point of view. In the simulations, we have found that the distributions of plant and wrack are always of the same type. As a result, we can restrict our analysis of pattern formation to one distribution (in this paper, we show the distribution of wrack, for instance). In Fig. 2, we show the evolution of the spatial pattern of wrack at t ¼ 0; 100; 350 and 550, with random perturbation of the steady state P and W . The random initial distribution leads to the formation of a strongly irregular transient pattern in the domain. At the final (cf. Fig. 2(C)–(D)), regular spots prevail over the whole domain and the dynamics of the system does not undergo any further changes. The temporal development of the spatially averaged values of the wrack and plant density shows that in the first intervals of simulations these values change fast the variance increases (cf. Fig. 3(a)–(b)). One can see that at t ¼ 250, the variance reaches a constant values and the mean increases slower. At this time, the pattern formation has almost completely evolved. And the system has reached its steady state at t > 350.
Mean value of the wrack density
0.334 0.332 0.33 0.328 0.326 0.324 0.322 0.32 0.318 0
100
200
300 Time
400
500
Mean value of the plant dentisy
0.222
0.22
0.218
0.216
0.214
0.212 0
100
200
300 Time
400
500
Fig. 3. The temporal development of the spatial average and variance of wrack and plant density plotted against time. (a) The time series of the mean value of the wrack; (b) The time series of the mean value of the plant. Parameter values are used as the same in Fig. 3.
988
Z.-H. Liu et al. / Applied Mathematics and Computation 216 (2010) 983–988
5. Discussion and conclusion In this paper, a spatial plant–wrack model with spatial diffusion and is investigated by mathematical analysis and numerical simulations. We found that the model exhibits spotted pattern, which may indicate that the plant and wrack are selforganized in the space. A question naturally arising is whether the number of the grid, used in numerical simulations, has great effect on the spatial pattern of the plant and wrack. To address this issue, we performed extensive computer simulations for different parameter values of the number of the grid (in total, more than 10 parameter sets were examined). We found that the behaviors are very similar (for the sake of brevity, we do not provide separated figures). Acknowledgments Supported by the PhD Research Fund of the Ministry of Education of China under Grant No 20050204JJ, and Scientific and Technological Project of the Shanxi Province (20080209ZX). References [1] M.R. Aguiar, O.E. Sala, Plant structure, dynamics and implications for the functioning of arid ecosystems, Trends Ecol. Evol. 14 (1999) 273–277. [2] P. Andresen, M. Bache, E. Mosekilde, G. Dewel, P. Borckmanns, Stationary space-periodic structures with equal diffusion coefficients, Phys. Rev. E 60 (1999) 297–301. [3] J. Bascompte, R.V. Sole, Modeling Spatiotemporal Dynamics in Ecology, Springer-Verlag, Berlin, 1998. [4] J. von Hardenberg, E. Meron, M. Shachak, Y. Zarmi, Diversity of vegetation patterns and desertification, Phys. Rev. Lett. 87 (2001) 198101. [5] R. HilleRisLambers, M. Rietkerk, F. van den Bosch, H.H.T. Prins, H. de Kroon, Vegetation pattern formation in semi-arid grazing systems, Ecology 82 (2001) 50–61. [6] C.A. Klausmeier, Regular and irregular patterns in semiarid vegetation, Science 284 (1999) 1826–1828. [7] S.P. Kuznetsov, E. Mosekilde, G. Dewel, P. Borckmans, Absolute and convective instabilities in a one-dimensional brusselator flow model, J. Chem. Phys. 106 (1997) 7609–7616. [8] J. van de Koppel, C.M. Crain, Scale-dependent inhibition drives regular tussock spacing in a freshwater marsh, Amer. Nat. 168 (2006) E136–E147. [9] R. Lefever, O. Lejeune, On the origin of tiger bush, Bull. Math. Biol. 59 (1997) 263–294. [10] R. Lefever, O. Lejeune, P. Couteron, in: P.K. Maini, H. Othmer (Eds.), Mathematical Models for Biological Pattern Formation, IMA, Springer, New York, 2000. [11] H. Malchow, S.V. Petrovskii, E. Venturino, Spatiotemporal Patterns in Ecology and Epidemiology: Theory, Models, and Simulations, Chapman and Hall, CRC, London, 2008. [12] J.D. Murray, Mathematical Biology II: Spatial Models and Biomedical Applications, Springer, New York, 2003. [13] D. Tilman, P. Kareiva, Spatial Ecology, Princeton University Press, Princeton, 1997. edited book. [14] C. Valentin, J.M. dHerbes, J. Poesen, Soil and water components of banded vegetation patterns, Catena 37 (1999) 1–24.