- Email: [email protected]
Turing patterns in a simple model of a nutrient–microorganism system in the sediment
Ecological Complexity 1 (2004) 77–94
Turing patterns in a simple model of a nutrient–microorganism system in the sediment Martin Baurmann∗ , Ulrike Feudel ICBM, Carl von Ossietzky Universität, PF 2503, 26111 Oldenburg, Germany Received 2 September 2003; received in revised form 19 December 2003; accepted 11 January 2004
Abstract Sediments are characterized by heterogeneous distributions of nutrients and microorganisms which emerge as a result of the interaction between chemical and biological processes with physical transport. We study in a simplified model the dynamics of one population of microorganisms and its nutrient, taking into account that the considered bacteria possess an active as well as an inactive state. Furthermore, the nutrients are transported by bioirrigation. It is shown that under certain conditions Turing patterns can occur which yield heterogeneous vertical spatial distributions of species. Furthermore, this model exhibits several stable coexisting spatial profiles, so that it depends crucially on the initial condition which of the distributions will be realized. This phenomenon of multistability can still be observed when spatial profiles are externally imposed by considering a depth-dependent bioirrigation. © 2004 Elsevier B.V. All rights reserved. Keywords: Reaction–diffusion models; Turing instability; Pattern formation; Sediment
1. Introduction The study of pattern formation in spatio-temporal systems has continued to be a central issue to modelers since the fundamental paper of Turing (1952). He analyzed the spontaneous emergence of inhomogeneous distributions of chemical substances in reaction–diffusion systems (RDS) and their possible importance for morphogenetic processes. Based on his seminal paper about diffusion instability a broad theory of pattern formation in RDS has been developed (Nicolis and Prigogine, 1977; Nicolis, 1995; ∗
Corresponding author. Fax: +49-441-798-3404. E-mail address: [email protected] (M. Baurmann).
Murray, 1993; Petrovskii et al., 2003; Satnoianu et al., 2000). A necessary condition for a diffusion instability to occur is the difference in diffusion constants of two nonlinearly reacting chemicals (activator and inhibitor). Since this is impossible to realize in solutions it took more than 30 years before Turing’s theoretical findings have been verified in experiment (Castets et al., 1990). In the last decade, the analysis of pattern formation has been expanded by including advection as a transport process (Rovinsky and Menzinger, 1994; Kuznetsov et al., 1997; Satnoianu et al., 1998, 2001). In particular, it has been shown, that the instability of a homogeneous distribution of chemical substances can be achieved by the interaction of nonlinear, not necessarily autocatalytic, reaction of three
1476-945X/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ecocom.2004.01.001
78
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
chemical species with advective transport (Rovinsky and Menzinger, 1994). Recently, it has been shown by Henry and Wearne (2002) that anomalous diffusion of the activator gives rise to Turing-patterns even if the diffusion coefficients of activator and inhibitor are equal. Pattern formation has also been a matter of research in biology and ecology. In particular, the aggregation behavior of microorganisms was studied. This paper focuses on possible mechanisms for the emergence of spatial patterns, i.e. inhomogeneous distributions of bacteria and nutrients in marine sediments. We show that the interaction between growth processes of bacteria, fed upon nutrients, and diffusion can lead to spatial pattern formation in the sediment. Our approach yields a qualitative understanding of the appearance of inhomogeneous spatial distributions of bacteria and nutrients using a very simplified model which reflects some fundamental properties of bacteria populations. Although this work has been inspired by experimental studies which show that chemical substances like sulfate or DOC are inhomogeneously distributed in depth (Mudryk et al., 2000; Madani et al., 2003), our goal is not to reproduce certain experiments quantitatively, but to explore principles of emerging patterns in marine sediments. The model we develop in this paper describes a type of interaction, that is an integral part of the complex ensemble of interactions found in benthic systems (see VanCappellen and Wang, 1996; Hunter et al., 1997). In contrast to these complex models only one bacteria species and its nutrient is considered in our model system. As a fundamental aspect of the model we divide the population of bacteria in an active and a dormant part. Additionally, we assume that the activation of dormant bacteria is mediated by the exchange of certain signal molecules. This assumption is motivated by several findings reported in litterateur: Kitsunezaki (1997) showed that the surfaces of bacteria colonies can propagate with different velocities yielding a fingering of the moving fronts. This special behavior is the result of the interaction of nonlinear diffusion and an activation mechanism of the bacteria. Similarly, Tsimring et al. (1995) divided the considered bacteria population in a motile and nonmotile part. In their model chemotaxis is the basic mechanism that induces the formation of spatial patterns. Camazine et al. (2001) showed that in slime
molds the chemotaxis, induced by a chemoattractant (cyclic adenosine monophosphate, cAMP), can lead to pattern-formation. The same substance was found by Bruns et al. (2002, 2003) to effect a stimulation of the growth of dormant cells in special planktonic bacteria in fresh and sea water and the authors state, that signal molecules (like cAMP) may be of relevance for the growth of a broader spectrum of marine bacteria. The paper is organized as follows: in Section 2, we introduce a simple model considering a bacteria species in a sediment that feeds upon one chemical substrate. To identify some basic properties of this model, we restrict ourselves first to a local model, in which only the reaction part of the system is taken into account, while diffusion is neglected. We find two equilibrium states and discuss their stability properties with respect to the variation of several parameters. These equilibrium states correspond to homogeneous distributions of species. While the equilibrium characterized by the extinction of bacteria is always stable, the equilibrium corresponding to positive population densities can lose its stability either in a Hopf-bifurcation (hb) or in a turning point (tp). Section 3 is devoted to the formation of spatial structures on a vertical 1D domain. We show that under certain conditions Turing structures appear which correspond to an inhomogeneous distribution of nutrients and microorganisms. These inhomogeneities emerge due to a diffusion instability without any spatial external forcing. This phenomenon presumes that the nutrient is much more mobile than the bacteria. Additionally, a certain supply of nutrient and a limit of the bacteria’s mortality must be guaranteed. Finally, we consider a spatial external forcing by incorporating boundary conditions or a bioirrigation function which decreases exponentially with depth. The aim of this part of our study is to check, whether pattern formation due to diffusion instabilities is robust against the influence of an external force. Although the stationary profiles differ from those found for the model with no spatial forcing, the general situation does not change too much: again, we find the coexistence of different spatially heterogeneous profiles in depth indicating that external forces cannot prevent pattern formation in large parameter intervals.
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
2. The minimal sediment model 2.1. The local model and the transport processes This paper focuses on a simplified model of the interaction of microorganisms and a chemical substance in the sediment. The microorganisms in the sediment build up their biomass by degrading special chemical species, using them as a nutrient. Consequently both, the population density of microorganisms, more specific of bacteria, and the distribution of chemicals in the sediment, is influenced by this interaction. Many different microorganisms and chemical substances are involved in the biogeochemical processes in the sediment (see e.g. VanCappellen and Wang, 1996; Hunter et al., 1997). Instead of examining the whole multitude of different pathways, we concentrate on a simple predator–prey model that we deploy as a prototype that exhibits some of the typical behavior,
79
found in the complete system. Thus, we set up a minimal sediment model (MS-model) in which we restrain to observe only one chemical species (concentration: N (mol/m3 )) acting as a nutrient and a population of one kind of bacteria (density: B (kg biomass/m3 )) that feeds on the nutrient. In the MS-model, five different processes are taken into account for the interaction between bacteria and nutrient (see Fig. 1): we consider two sources of input of nutrient into the system. On one hand, we include the diffusive flux through the sediment’s surface into the system, on the other hand, chemicals can be transported in lower regions of the sediment by bioirrigation leading to a nonlocal transport term. The interaction between bacteria and the chemical species consists of a consumption of nutrient and the corresponding increase of biomass of the bacteria due to the degradation of nutrient. As a last factor we include the mortality of microorganisms.
Diffusion
Boundary Flux
Bioirrigation
Seawater
local Processes Sediment
bioirrigation
Nutrient
consumption
Fig. 1. Minimal sediment model.
growth
mortality
Bacteria active inactive part part
(de)activation
80
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
The effects of bioirrigation are modeled by a diffusive input of Fickian type. It is sufficient to parameˆ − N), terize this diffusive term by a certain flux Φ(N ˆ where N is the nutrient concentration in the sea water. Thus, the nonlocal transport is considered as a diffusive connection between the sediment in each depth with the sea water above the sediment. In the model, this process is realized as a lateral boundary condition and, hence, a local process at each depth level. When examining the spatial model we will replace the fixed parameter Φ by a function Φ(z), where z is the depth below the sediment surface. By this means a decreasing effect of bioirrigation with increasing depth can be incorporated into the model. The flux of nutrient through the sediment surface is another boundary condition that has to be applied to the model’s uppermost point. It will be considered in a spatial model, whereas we will neglect it in the local formulation. To improve the consistency of the model with natural systems, a certain part of the bacteria population will be considered as inactive. Only the active part of the bacteria Ba feeds on the nutrient leading to a certain growth rate and dies with a mortality m, while the inactive part of bacteria may die but will not effect the nutrient consumption. The amount of active bacteria (Ba ) is not bound by a rigid factor but is determined by a (dimensionless) function r(B), which depends only on the biomass of microorganisms. The actual form of r is difficult to specify. There exist no well established models for the activation-mechanism of bacteria in the sediment. Nevertheless, we can state, that the biomass of active bacteria cannot exceed the total population density. So, r(B) is restricted to: 0 ≤ r(B) ≤ 1
(1)
Additionally, we expect an increase of r if B increases, presuming some activation mechanism induced by the organisms of B. At least for certain ranges of biomass B we have to guarantee: ∂r >0 (2) ∂B By that approach we model an activation mechanism that agrees with the results found by Bruns et al. (2002, 2003) in experiments with fresh and sea water bacteria. Realizing (1) and (2), we employ the function: B r(B) = (3) B+L
That means we assume this function to be of Monod-type. This type of function is one of the simplest mathematical functions possessing a saturation characteristics. In principle, any other function for which (1) and (2) holds is possible to choose and would give the same results. Applying this formula, the system is governed by the following equations: dB = αBa N − mB dτ dN ˆ − N) − βBa N = Φ(N dτ B Ba = r(B)B = B B+L
(4)
where α is the rate of conversion into biomass, β is the capturing rate, m is the mortality rate and L is the ˆ parameterize the exhalf saturation constant. Φ and N ˆ is the concentration of change due to bioirrigation: N N in the external pool of sea water and Φ parameterizes the hydraulic conductivity between sediment and sea water. All listed parameters have positive values. To reduce the number of model parameters and to make the Eq. (4) dimensionless we rescale them. Using m Y, α Φ L = K, β
N=
m Yˆ , α 1 τ= t Φ
ˆ = N
B=
Φ X, β
our resulting system of equations reads: m X dX = XY − X dt Φ X+K X dY = Yˆ − Y − XY dt X+K
(5) (6)
This is the dimensionless three parameter formulation (m/Φ can be regarded as one independent parameter) of the local MS-model, which shall be examined in the following. Though the consumption of nutrient and growth of bacteria’s population are modeled by a Lotka–Volterra approach, both are highly non-linear with respect to X, because the feeding term depends on the biomass of active bacteria Ba only. Therefore, the MS-model can be considered as a predator–prey system with a non linear predator-dependent response function.
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
2.2. Dynamical properties of the local model The analysis of the dynamical properties of the local model is a necessary assumption to understand the possible states of the homogeneous system without local transport. The system is in equilibrium when the variables do not change in time. It has then reached a stationary state. The locations and characteristics of these states are of high interest for us, because they give answers, how the system will behave in the long time run. Particularly we are interested in the question for which values of the parameters the population of bacteria will be able to survive. Analyzing the model equations, we find a trivial stationary state at [X0∗ = 0, Y0∗ = Yˆ ], at which the bacteria are extinct and the nutrient has reached the concentration of the sea water (Yˆ ). There might exist two other stationary states in which nutrient and microorganisms coexist at ∗ , Y ∗ ), where (X1,2 1,2 ∗ X1,2 = 21 (Yˆ − 1 ± (Yˆ − 1)2 − 4K) (7) and Yi∗ =
Xi∗ + K Xi∗
i = 1, 2
(8)
Regarding (7), the local model will show these additional stationary states, if (Yˆ − 1)2 − 4K ≥ 0
(9)
Inequality (9) indicates the ranges where different stationary states can coexist for a given set of parameters: √ K ≤ 41 (Yˆ − 1)2 or Yˆ ≥ 1 + 2 K The parameters m and Φ have no impact on the ranges of coexistence. Since we examine a natural system, it is fundamental, to check, whether stationary states are stable with respect to small perturbations or not. If they are not, even small changes in the nutrient’s concentration or the bacteria’s biomass will destabilize the equilibrium and the system’s variables will converge to a stable stationary state. Due to the inherent fluctuations of natural systems unstable stationary states will only exist in a mathematical context and cannot be observed by measurements or experiments. Nevertheless, they
81
might be of importance for the dynamics of the system as we will discuss later. At this point, we will not reproduce the whole analysis, that can be used to determine a stationary state’s stability; it is described in detail in a large number of publications (e.g. Murray, 1993; Nicolis and Prigogine, 1977; Arrowsmith and Place, 1992). The stability of a stationary state is determined by the signs of the real parts of the eigenvalues computed from the Jacobian matrix at the steady state. If the real parts of all eigenvalues are negative, then the steady state is stable, otherwise not. For the MS-model we obtain two eigenvalues: √ m (T ± T 2 − 4S) (10) λ1,2 = Φ 2(Xi∗ + K) with Φ Yˆ Xi∗ , S = (Yˆ − 1)Xi∗ − 2K (11) m According to the number of negative real parts of the eigenvalues we can classify the stationary state. If both eigenvalues have negative real parts the steady state is stable and called attractor. If at least one real part of the eigenvalues is positive, the steady state is unstable and called a saddle (one positive and one negative real part) or repellor (two positive real parts). Let us first discuss the stability of (X2∗ , Y2∗ ). In this case we obtain S < 0. Hence, for the steady state (X2∗ , Y2∗ ) both eigenvalues of the Jacobian are real and have different signs, so that (X2∗ , Y2∗ ) is a saddle. Contrary, at (X1∗ , Y1∗ ) we have S > 0. The real parts of both eigenvalues have the same sign, which is just the negative of the sign of T . Consequently we can distinguish: T =K−
if X1∗ >
mK Φ Yˆ
[X1∗ , Y1∗ ] is an attractor
if X1∗ <
mK Φ Yˆ
[X1∗ , Y1∗ ] is a repellor
We omit substituting the merely complex algebraic expression for X1∗ . Instead we sketch the two qualitatively different stability properties instead, in order to illustrate all possible transitions of stationary states. For that purpose, we will vary the parameter Yˆ , corresponding to the concentration of nutrients in the sea water, whereas for the other parameters we choose two different sets of constants listed in Table 1.
82
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
Table 1 Parameter sets for the MS-model
Parameter set 1 Parameter set 2
X*
K
m
Φ
0.8 8.0
0.5 8.0
1.0 1.0
14 12 10 8 6 4
These sets differ in the values of K and m. m is the mortality rate, so that for parameter set 2 the microorganisms are more likely to die. The parameter K affects the activation ratio of the bacteria population. If the value of K is relatively small (parameter set 1), a slight increase of bacteria’s biomass causes an almost total activation of microorganisms, whereas a higher value of K (parameter set 2) results in a slower activation. Generally, the ratio of active bacteria is (at same amounts of biomass) higher, if K is smaller. Fig. 2 illustrates the situation for the two different parameter sets depending on the concentration of nutrients in the sea water Yˆ . We can distinguish two kinds of qualitative behavior corresponding to different environmental conditions. Let us first discuss the stationary states and their stability. For both parameter sets the trivial steady state corresponding to an extinct bacteria population (X0∗ , Y0∗ ) is stable in the whole Yˆ -interval under consideration. But for a part of this interval there exist the two other ∗ steady-states X1,2 which are created in a turning √ point at Yˆ tp = 1 + K. X1∗ and X2∗ exhibit different stability properties. For parameter set 1 (Fig. 2(a)) the turning point corresponds to a saddle-node bifurcation leading to the emergence of a stable attracting equilibrium state (X1∗ , Y1∗ ) and an unstable saddle (X2∗ , Y2∗ ). Thus, the system possesses two coexisting stable equilibria beyond the turning point (Yˆ > Yˆ tp ), i.e. the system is bistable. In this case it depends crucially on the initial conditions which of the two stable states (X0∗ , Y0∗ ) or (X1∗ , Y1∗ ), will be realized by the system in the long-term limit. Starting with a nonzero bacteria population we can change the environmental conditions in such a way that we decrease the nutrient concentration in the sea water. When reaching the turning point, we observe a sudden collapse of the bacteria population due to the loss of stability. As a result the bacteria population becomes extinct.
2 0 (a) 0
2
4
6
8
^ 10 12 14 Y
2
4
6
8
^ 10 12 14 Y
X* 14 12 10 8 6 4 2 0 (b) 0
Fig. 2. Stability of stationary states for two different parameterizations and variable value Yˆ (axis of abscissae). Stable states are marked by thick solid lines, unstable by thin dashed lines. The box symbolizes a turning point, the star a Hopf-bifurcation. From the Hopf-bifurcation point a steep dotted line emerges, signifing a set of periodical points on a unstable orbit.
For parameter set 2 we observe bistability only for ∗ exist. Stability a part of the Yˆ -interval, where X1,2 of a nonzero bacteria population occurs only for high values of the nutrient concentration in the sea water Yˆ . Decreasing Yˆ leads again to a collapse of the bacteria population, but in a different dynamical scenario. The stationary state loses stability in a Hopf-bifurcation instead of a saddle-node bifurcation. The Hopf-bifurcation leads to the emergence of periodic behavior (dotted line in Fig. 2b). However, these oscillations cannot be observed in a natural system because they are unstable (subcritical Hopf-bifurcation). A Hopf-bifurcation occurs if a pair of complex conjugate eigenvalues of the Jacobian crosses the imaginary axis. Presuming m 2Yˆ > Φ Yˆ − 1
(12)
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
83
Fig. 3. The parameter space of the MS-model. Parameter sets corresponding to a turning point situation are symbolized by the bright surface, those connected with a Hopf-bifurcation lie on the dark surface. The thick black line is the line where both surfaces touch each other. It is the set of Takens Bogdanov points.
this can be obtained by setting Yˆ Φ Yˆ Φ Yˆ − 1 − K = Khb := m m
(13)
Fig. 3 illustrates the different kinds of dynamical behavior in the three dimensional parameter space spanned by (K, m/Φ, Yˆ ) The set of turning points and Hopf-bifurcations are the surfaces Ktp and Khb , respectively. These surfaces divide the parameter space into three subspaces. Parameter sets corresponding to points which lie above the Ktp -surface designate pa∗ , Y ∗ ) do not exist. If (12) rameter sets where (X1,2 1,2 holds and Khb < K < Ktp (space between the surfaces) the non-trivial stationary state exists and is a repellor. Only the points below both surfaces stand for parameter sets that generate bistable dynamics, because besides of the trivial equilibrium (0, Yˆ ) also (X1∗ , Y1∗ ) is an attractor. Since we have to fulfil (12) the Khb -surface is bounded. At its boundaries it touches the Ktp -surface and is not continued across. Both surfaces touch in a line, that consists of so called Taken Bogdanov points. They symbolize parameter sets, showing eigenvalues λ1 = λ2 = 0 at (X1∗ , Y1∗ ) = (X2∗ , Y2∗ ) = √ √ (2 K, 2 + K).
Fig. 3 shows that there are two fundamental ways, how the steady state of nonzero bacteria’s biomass may loose its stability, when environmental parameters change. The first one occurs by crossing a turning point (tp), the second one coincides with traversing of a subcritical Hopf-bifurcation. Both diagrams in Fig. 2 correspond to different lines in Fig. 3. Let us discuss the distinct ways of destabilization also from another point of view. As already mentioned above bistability means that the final state approached by the system in the long-term limit depends on the initial condition. To study the robustness of the system with respect to larger perturbations, it is useful to analyze the basins of attraction of the different coexisting stable states. The basin of attraction is the set of initial conditions which converges to one stable state of the system. In our example we have two basins of attraction for the stable nonzero state (X1∗ , Y1∗ ) and the extinct state (X0∗ , Y0∗ ), respectively. Fig. 4 shows the two stable stationary states and their basins of attraction for values of Yˆ close to the critical values for the loss of stability (Yˆ ≈ 2.78 for tp and Yˆ ≈ 9.1 for hb). In the hb-type system, the basin boundary is made up by the unstable periodic orbit emerging at the Hopf-bifurcation, while in the tp-type
84
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
Fig. 4. Steady states and their basins of attraction; series A: Parameter set 1 and Yˆ = [3.0, 2.9, 2.8] (tp-type); series B: Parameter set 2 and Yˆ = [9.5, 9.25, 9.15] (hb-type). Stable equilibira are marked by a +, unstable ones by a 䊊.
the basin boundary is connected with the unstable equilibrium. In the latter case, the role of the unstable equilibrium becomes evident. Even though it is not observable in nature, its location is important for the dynamics of the system. The unstable steady-state
(saddle point) “sits on” the boundary of the basin of attraction, its stable manifolds make up the basin boundaries. Both diagrams illustrate that the states are moving closer to their basin boundaries, while Yˆ decreases.
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
The difference between the two examined ways is, that the size of hb-type’s basin converges to zero, while the tp-type’s basin collapses abruptly. According to that mathematical finding, natural systems representing the hb-type are more likely to lose stability, than their tp-type counterparts, since a fluctuation to any direction will cause the bacteria to become extinct, if it is only large enough. As a consequence hb-type-systems with parameters close to the critical bounds can be destabilized by reasonably small perturbation even in the case, when they are still stable.
3. The spatial minimal sedimentmodel 3.1. Setup of the spatial model In the sediment, the local processes analyzed in Section 2 are influenced by transport processes like diffusion leading to spatial distributions of chemicals and bacteria populations. Therefore, we now consider a spatially extended system in order to include the effects of diffusion. In particular, we are interested in the changes of the stability properties of homogeneous species distributions due to diffusion. We apply the model on a vertical one-dimensional domain with space coordinate ζ (later with dimensionless coordinate z). The diffusion process takes place along the vertical direction and is modeled by Fickian law. At the lower and upper boundaries we have to establish boundary conditions: we assume that we have no fluxes over the lower boundary and accept a certain flux of nutrient at the sediments surface, which forms the upper boundary. As another step towards a realistic model we will take into account, that effects of bioirrigation will decrease with increasing depth. To incorporate this phenomenon, we exchange the fixed parameter Φ by a monotonously decreasing function ζ Φ(ζ) = Φ exp −a (14) ζl where a is a non-negative factor parameterizing the exponential function and ζl is the length of the domain. Inducing these changes in the original equations (formulated with respect to B and N) and scaling (z =
85
ζ/ζl )leads to: m ∂2 X XY ∂X = X − 1 + DX 2 ∂t Φ X+K ∂z ∂Y X ∂2 Y = exp(−az)(Yˆ − Y) − XY + DY 2 ∂t X+K ∂z (15) and ∂Y = Φsurf (Yˆ − Y) ∂z z=0
∂Y =0 ∂z z=1
(16)
with the conductivity across the surface Φsurf and the dimensionless diffusion coefficients: DB DN DX = and DY = (17) 2 Φζl Φζl2 To simplify the differential equations we will, excluding the case DX = 0, implement δ with δ = DY /DX . This substitution is helpful, because the ratio of the diffusion coefficients is of more interest than their absolute values. √ In case of constant bioirrigation a scaling by z = DX ξ leads to: ∂X m XY ∂2 X = X −1 + 2 ∂t Φ X+K ∂ξ (18) ∂Y X ∂2 Y = Yˆ − Y − XY + δ 2 ∂t X+K ∂ξ 3.2. Diffusion-induced patterns in the spatial model In this section, we want to examine the effects of diffusion only, so we set Φsurf = 0 and a = 0. In that case the conditions of nutrients to flow into the system are equal in all points of the domain. Since diffusive fluxes are driven by gradients, they only affect the solution, if gradients occur. As a consequence, all homogeneous profiles referring to the local model’s solutions will also solve (18). Thus, one solution of (18) can always be found: these are the trivial homogeneous functions (X = 0, Y = Yˆ ) for all values z. As we have shown above an additional non-trivial, homogeneous solution may exist in particular intervals of parameters. Besides these two homogeneous solutions no other profile without concentration gradients exists. Again we look at the stability of the homogeneous depth profiles of the concentration with respect to perturbations. We will show in the following, that the
86
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
non-trivial homogeneous stationary state may loose its stability, if diffusion is taken into account as an additional process. Instead non-homogeneous concentration profiles will arise. At first sight this possibility seems to be absurd: diffusion is known as a process, that will extinguish spatial gradients. However, it has been shown by Turing (1952) that nonlinear chemical reactions in spatially extended systems can lead to an inhomogeneous distribution if diffusion is taken into account. These Turing instabilities are excited by perturbations of specific spatial wavelengths and will only appear within certain parameter spaces (Alonso et al., 2002). Our aim is to identify these Turing-spaces with respect to parameters and wavelengths and to discuss associated system properties. With respect to the spatial MS-model we study the question: under which conditions is it possible, that homogeneous profiles (X(z) = X1∗ , Y(z) = Y1∗ ) become unstable, presuming (X1∗ , Y1∗ ) is an attractor in the local model? To answer this question we analyze the dynamics of the spatio-temporal system (15) with respect to a perturbation modeled with a wave ansatz (compare Murray, 1993). We consider perturbation functions x = ζx exp(λt − x), i2πκz = ζy exp(λt − x) and investigate their time evolution in the linearized model. Taking into account the chosen boundary conditions, we can restrict ourselves to the following class of functions: x(z, t) = cx cos (2πκz) exp(λt)
1/2 j11 DY +j22 DX
− (j11 DY +j22 DX )2 −4DX DY det(J) 1 κˆ − = 2π 2DX DY (21) and
1/2 j11 DY + j22 DX
+ (j11 DY +j22 DX )2 −4DX DY det(J) 1 κˆ + = 2π 2DX DY (22) where J denotes the Jacobian matrix of the local model at (X1∗ , Y1∗ ) and its elements ˙ ∂X j11 = ∂X X∗ ,Y ∗ 1
and j22
1
∂Y˙ = . ∂Y X∗ ,Y ∗ 1
1
To avoid the complex terms of κˆ − and κˆ + , we first check that j22 < 0 and gain a weaker formulation of (20) by the estimate: Re(λ(κ)) > 0 ⇒ κ− < κ < κ+
y(z, t) = cy cos (2πκz) exp(λt)
with
Linear stability analysis yields intervals for the wavelength κ, within which the homogeneous solution is unstable with respect to perturbations with this wavelength. Again the eigenvalues of the corresponding Jacobian which depend now on the wavelength κ and the diffusion constants DX and DY have to be positive. For reaction–diffusion systems like that of the MS-model, we have (compare Petrovskii et al., 2003): Re(λ(κ)) < 0 if η := (j11 DY + j22 DX )2 − 4DX DY det(J) < 0 (19) If η > 0 it follows: Re(λ(κ)) > 0 ⇔ κˆ − < κ < κˆ +
with
(20)
1 κ− = 2π
(23)
det(J) , DY j11
1 κ+ = 2π
j11 DX
(24)
Since the considered equilibrium is an attractor in the local formulation, we have det(J) > 0. Hence, we need j11 > 0 to fulfil (24). This means that the increase of the bacteria’s biomass has to be autocatalytic. We obtain this special property in the MS-model because the activation of bacteria is more pronounced for large populations, since we assume that bacteria are capable to activate dormant cells by excreting signal molecules. κ− < κ < κ+ is a necessary condition to obtain Turing instabilities. By this inequality we can outline a space in the κ −DY -plane, where Turing instabilities
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
87
may be obtained. Outside this space, no such instabilities occur. In particular, we have no Turing instabilities, if κ− > κ+ . Substituting (24) we obtain: det(J) (25) δ < 2 ⇒ no Turing instabilities j11
Condition (26) can be formulated stricter: by derivating Re(λ(κ)) with respect to κ, we find a maximum at 1 DX j22 + DY j11 κe = (27) 2π 2DX DY
So for the MS-model we can state:
Because Re(λ(κ)) ≤ Re(λ(κe )) it is obvious that there is a Turing space only if Re(λ(κe )) > 0 demanding as a necessary condition:
Φ(X1∗ + K)[(Yˆ − 1)X1∗ − 2K)] δ< mK2 ⇒ no Turing instabilities
(26)
Eq. (26) shows, that the increase of Yˆ , Φ and DX leads to stabilization. On the other hand, increasing K, m and DY destabilizes the system. Altogether we can state that the occurrence of Turing instabilities is enhanced by a low inflow of nutrient into the system and a high loss of bacteria due to mortality. 12 κ
+
10
κ
8
κe
6 4 2κ 0 −5
−4
(a) 10
14 12
10
−3
10
Dx
κ
10 κ
8 6 4κ + 2 0 −2
(b) 10
κe
10
−1
0
10
DY
Fig. 5. κ− , κ+ and κe for different values of DX and DY , respectively. Outside the area limited by the thick lines, no Turing instabilities can be met.
κ− < κe < κ+
(28)
Since we can show that the inequality κe < κ+ holds on the entire parameter space, the essential condition is κ− < κe (see Fig. 5) leading to: δ<
det(J) − j12 j21 ⇒ no Turing instabilities 2 j11
(29)
3.3. Simulations on a discretized domain The analytical examination of the sediment model led to the important result, that an attracting homogeneous state of nonzero bacteria population may loose its stability, if diffusion is taken into account. Beyond the diffusion instability inhomogeneous distributions can occur which cannot be computed analytically. To illustrate the emergence of inhomogeneous distributions of bacteria and nutrient, we consider 41 boxes. We use the parameter set 1 and set Yˆ = 8, DX = 2 × 10−4 and δ = 5000. The initial values are the homogeneous profiles that are perturbed by a cosine function with five periods over the domain length (pattern type 0). Fig. 6 shows, that for parameter set 1 the initial spatial frequency of the profile’s perturbation is changing quickly in time and finally a pattern will form, consisting of two cosine cycles (pattern type 1). The transient dynamics is rather short, so that already after time t = 60 the final distribution exists. To illustrate how a change of bacteria’s mobility will affect the pattern formation process, we will change now the value of δ. For a merely small interval of δ ∈ [0, 5000] a homogeneous profile will appear in the longtime-run (see Fig. 7). This spatially constant profile will loose its stability when δ increases. An abrupt transition to two maxima occurs after crossing the Turing instability threshold. Choosing higher values 5000 < δ < 6000,
88
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
Fig. 6. Long-term behavior of a spatial cosine oscillation in the MS-model (parameter set 1, Yˆ = 8, DX = 2 × 10−4 , δ = 5000).
the amplitudes of the profile are growing rapidly until they reach a range in which the increase is moderate 6000 < δ < 13,000 (pattern type 1). At a value of δ = 13,000, the qualitative behavior of the MS-model’s long-term-solution changes abruptly again: instead of the spatial oscillation with two cosines, a two-modal profile occurs (pattern type 2). Fig. 8 illustrates the abrupt change of the stationary states. Additionally, it is interesting to note, that the small increase of δ causes the profile to show a short spatial interval, where X is extinct. This interval widens, if δ is increased further. The consequence of the appearance of this “dead-zone” in the middle of the domain is that at this point, the domain is separated into two subdomains with no diffusive mass transfer between them. This would mean that life of microorganisms in deeper layers of the sediment may be maintained only by nutrient inflow due to bioirrigation. Further analysis shows that the two different types of inhomogeneous profiles, shown in Fig. 8 coexist
for a broad range of δ-values. Thus, we can conclude, that the transition from pattern type 1 to pattern type 2 at δ = 13,000 in Fig. 8 is not due to a bifurcation but rather to a transition of the initial condition of our simulation from one basin of attraction to another one. To study, which stationary states exist for the MS-model and how they change, when the parameter values are modified, we use the program Content (Kuznetsov and Levitin, 1996) to follow equilibrium points under variation of DY . Content allows the analysis of one-dimensional PDEs at a user defined discretization of the spatial domain. We use parameter set 1 and chose a discretization (1000 boxes) fine enough that our results are stable with respect to the discretization length and the detection of equilibrium points that are numerical artifacts is avoided. Fig. 9 shows a diagram in which the biomass of bacteria X in the uppermost discretization-box is plotted versus DY . It illustrates the emergence of stationary states from the homogeneous profile indicated by the
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
89
Fig. 7. The increase of δ effects (abrupt) changes of the equilibrium states of bacteria population (subplots only differ by the location of viewpoints).
90
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
20
X
15
10
5
0
0
0,2
0,4
0,6
0,8
1
depth Fig. 8. Stable profiles at δ = 13,000 (thin line) and δ = 13,001 (thick line).
horizontal line at X(z = 0) = 6.88, Y(z = 0) = 1.14. At each bifurcation point on that line (marked by squares) a non-homogeneous profile emerges. Since these profiles are characterized by a certain wavenumber (at least in the neighborhood of the bifurcation) we can assign this characteristic wavenumber to the bifurcation points. Knowing that Re(λ(κ)) = 0 at the bifurcation points we can evaluate the location of the bifurcation points using (20) and deriving the critical values DYc (κ). We obtain: κ 1/2 1 3/2 2 5/2 ≥3 DYc () 6.1614 1.7712 1.0370 1.0076 5.0632 <0
Fig. 9. Projection of the MS-model bifurcation diagram on the DY /X[0]-plane. Squares symbolize bifurcations (not all plotted). Most bifurcation branches are labeled due to their main oscillation mode (number) and their orientation at the origin (sign). The branch of homogeneous profiles we marked by 0, and branches with extincted population of bacteria in the middle part by exc. There exist also a pair of branches representing unimodal profiles with bacteria extinction at the boundaries (exb). They superpose the 2- and 3/2- branches displayed very close to the abscissa.
At DY ≈ 1.0, we find two bifurcation points in close succession. They can be assigned to the wavenumbers 2 and 3/2, respectively. At DY ≈ 1.77, profiles of wavenumber 1 emerge, followed by a bifurcation with κ = 5/2 at DY ≈ 5.06. Profiles showing an half cosine oscillation κ = 1/2 emerge at DY ≈ 6.16. Note that solutions on the κ = 1 branch do not exist for DY > 3.6 . . . , since the branch “disappears” in a bifurcation with the κ = 2 branch. In the same way, the branches characterized by κ = 1/2 meet in a bifurcation on the 2+ set at DY ≈ 60 (not shown).
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
All other branches can be followed to arbitrary high values of DY . At [DY ≈ 2.4, X ≈ 6.88], the projection shows a superposition of homogeneous profiles and another curve exc. This second curve represents profiles with zones of extinction (as shown in Fig. 7). There are no equilibrium points on this curve showing homogeneous distributions of microorganisms. At [DY ≈ 2.4, X ≈ 6.88] we find two distinct equilibrium points (one homogeneous, the other not) that are superimposed by projection. A similar branch of inhomogeneous solutions representing profiles with extinction of bacteria at the boundaries and a unimodal distribution of biomass in the middle can also be followed, but can hardly be distinguished from the abscissa in the projection plot. Additionally to the location of equilibrium profiles their stability with respect to perturbations is of interest. Unfortunately, Content does not provide the necessary eigenvalues. Nevertheless, we were able to get some information on stability by employing the program Candys/QA (applying a much coarser spatial discretization) (Feudel and Jansen, 1992). We could derive a similar bifurcation diagram but got additional branches, which we consider as numerical artifacts (due to the coarse discretization). Furthermore, we check the stability for some values of DY and some wavenumbers κ by numerical experiments: we start with the homogeneous profile and perturb it with a small oscillation (wavenumber κpert ). By forward integration we test to which stationary state the perturbed profile converges. Table 2 summarizes the results of our numerical experiments. We see, that the homogeneous profile is loosing stability in favor of profiles with two cosine oscillations (2+ and the phase shifted 2−). Increasing DY further this attractor coexists with stable 3/2+ and 3/2− profiles. At higher values of DY , profiles with a zone of extinction (exc+ and exc-) will become attractive, instead of the 2+ and 2− distributions (cf. Fig. 7). On the 1/2±, 1± or 5/2±
91
Fig. 10. Location of stationary states in a MS-model with nonuniform bioirrigation. The 3/2+ branch (compare Fig. 9) could not be found.
branches, we found no attracting intervals, they seems to be completely unstable. 3.4. Effects of non-uniform bioirrigation So far we discussed a model with uniform bioirrigation acting as a spatially homogeneous forcing. With respect to natural systems this is not realistic. As mentioned above, the input of nutrient—due to bioirrigation—becomes smaller with increasing depth. This induces an inhomogeneous external forcing and thus we expect an imposed external depth-dependence. Let us check whether this externally imposed forcing destroys the pattern formation process. To take the depth dependence into account we have to analyze the Eq. (15) with a > 0. Fig. 10 shows the bifurcation diagram for a model with a = 0.1. So, the effective hydraulic conductivity at the top of the domain is Φ(z = 0) = Φ whereas it is Φ(z = 1) = exp(−0.1)Φ = 0.9Φ at the lower boundary. Of course this is a very moderate decrease
Table 2 Numerical experiments: the homogeneous profiles were perturbed by a (cosine) oscillation of wavenumber κpert and converge to the noted branch in the bifurcation diagram (cf. Fig. 9) DY κpert Converges to branch
1.3 3/2 3/2+
1.5 2 2+
2.0 2 2+
4.0 2 2+
4.0 3/2 3/2+
6.0 2 exc+
6.0 3/2 3/2+
10.0 2 exc+
10.0 3/2 3/2+
10.0 5/2 exc+
92
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
of bioirrigation effects, but it allows us to compare the diagram with Fig. 9. Surveying both diagrams we are able to identify most of the branches (except of a part of the 3/2+ branch, that is missing in Fig. 10). In contrast to the model with uniform forcing the diagram displayed in Fig. 10 shows no bifurcations connecting branches. Pitchfork bifurcations split up and lead to the disconnection of branches. The reason for this splitting is a symmetry breaking, which is introduced by the depth-dependent bioirrigation. Additionally, there is a shift of the branches to the left, so that we can say that coexistence is more easy to achieve (i.e. with a smaller DY ). This observation is consistent with our findings that the decrease of inflow of nutrient leads to a weakened stability of the homogeneous profile. The depth-dependence of bioirrigation applies already a certain spatial structure in form of an external forcing. Surprisingly, the emerging spatial distribution of species beyond the instability do not show only the imposed spatial profile but a much broader variety of different spatial distributions. Thus, depthdependent external forcing do not prevent the emergence of spatial patterns different from the imposed ones. 3.5. Effects of boundary fluxes Very similar to Fig. 10 is the bifurcation diagram we obtain by running the model with uniform bioirrigation but including an upper boundary flux of nutrient (Fig. 11). Following the same strategy as in the last subsection, we limit the effects of the changed boundary conditions very strictly to allow the comparison with the bifurcation plot in the uniform case (Fig. 9). Fig. 11 shows some common features with Fig. 10. In particular, we again obtain symmetry breaking and hence a splitting up of bifurcation points. The shift of branches towards lower values (in comparison to the uniform case) is neglectible, although we would expect a shift to the right due to the stabilizing effects of higher inflow of nutrients. An important change occurs looking at the spatial profiles with a “dead-zone”: the formerly isolated exc+/ exc- branch has now contact to other branches (of type 3/2). So, when changing parameters a “dead zone” may appear without a sudden change of the system’s behavior but by the expansion of a zero point.
Fig. 11. Bifurcation diagram for the spatial MS-model with uniform bioirrigation and a small upper boundary-flux Φsurf = 0.01.
4. Discussion We have studied a minimal sediment model, describing the interaction of a bacteria population and its nutrient. In this model we take into account that bacteria can live in an active state, in which they reproduce consuming nutrient, as well as in a passive state of dormancy. Furthermore, we assume that the activation of bacteria is enhanced by the size of population. This assumption is supported by experimental studies which show a stimulation of the growth of certain marine bacteria by signal molecules. The local dynamics of this model, which can be regarded as a two-species predator–prey system (with input of prey, predation and the mortality of the predator) is characterized by the coexistence of two equilibria over a large range of parameters. One equilibrium is a state, in which the bacteria population is extinct. This (trivial) equilibrium is stable in the whole parameter range considered. Assuming that the activation of dormant bacteria is fast and the input of nutrient is high, there is another stable equilibrium, characterized by a positive population density of bacteria. Varying the parameters of the system this steady state can loose its stability by crossing a Hopf-bifurcation or a turning point. In the first case, the basin of attraction of the equilibrium is contracting to the equilibrium point, while in the second case it collapses abruptly. Thus, one can conclude that the system is more sensitive to perturbations in parameter regions in which
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
the Hopf-bifurcation is close. The reason is that the Hopf-bifurcation is subcritical leading to an unstable limit cycle whose location determines the boundary of the basin of attraction for the nontrivial equilibrium state. The main focus of this paper is the study of pattern formation processes due to the interaction of nonlinear growth processes of the bacteria and diffusive transport. We have shown that in our simple minimal sediment model Turing patterns can emerge. Neglecting diffusion and nonconstant forcing the only spatial profiles, that are stationary and stable, are those, that correspond to the equilibria, found for the local model. When diffusion comes into play, the non-trivial equilibrium can loose its stability in favor of a heterogeneous distribution of species due to a diffusion instability. We have derived the conditions for the instability to occur: the nutrient has to be much more mobile than the bacteria which is usually fulfilled in the sediment due to the stickiness of bacteria to the sediment matrix. Additionally, a minimal supply of nutrient must be guaranteed and the mortality of bacteria must be relatively low. A further analysis of the model proved, that several stable heterogeneous profiles can coexist. In that case it depends on the perturbation that is added to the constant profiles to which the model converges. We depicted the location of equilibria in the state-space in a bifurcation diagram (see Fig. 9). Further analyses of the spatial model yield the result that the minimal sediment model exhibits multistability, i.e. it possesses several stable coexisting heterogeneous profiles. In that case it depends crucially on the initial condition to which heterogeneous distribution of nutrient and bacteria the model converges. Thus, there are different spatial profiles possible under the same environmental conditions. Finally, we have studied the influence of a nonconstant spatial forcing on the pattern formation process. This forcing describes decreasing bioirrigation with increasing depth and/or boundary fluxes. Both situations are more realistic in natural environments than constant forcings. Surprisingly, the externally imposed spatial structures did not destroy the pattern formation process. We again find multiple coexisting heterogeneous profiles different from the one induced by the non-constant forcing. Qualitative changes in
93
the bifurcation diagram are only observed when the non-uniformity of the forcing is pronounced. It is important to note that the non-constant forcing leads to a symmetry-breaking which results in a split up of bifurcation points. Furthermore, the original homogeneous profile is replaced by an inhomogeneous one. According to our studies the appearance of Turing patterns can be regarded as one possible mechanism of the formation of heterogeneous spatial distributions of nutrients and microorganisms in the sediment. But there are many other processes which may influence the emergence of spatial patterns in the sediment and which are not taken into account here. Such processes include advection, chemotaxis and the competition between different degradation pathways. Future work needs to take these processes into account.
Acknowledgements We thank H. Cypionka, W. Ebenhöh and T. Gross for discussions. This work was supported by the DFG research group BioGeoChemistry of Tidal Flats.
References Alonso, D., Bartumeus, F., Catalan, J., 2002. Mutual interference between predators can give rise to Turing spatial patterns. Ecology 83 (1), 28–34. Arrowsmith, D., Place, C., 1992. Dynamical Systems—Differential Equations, Maps and Chaotic Behavior. Chapman & Hall/CRC Press, London. Bruns, U., Cypionka, H., Overmann, J., 2002. Cyclic AMP and acryl homoserine lactones increase the cultivation efficiency of heterotrophic bacteria from the central baltic sea. Appl. Environ. Microbiol. 68 (8), 3978–3987. Bruns, A., Nübel, U., Cypionka, H., Overmann, J., 2003. Effect of signal compounds and incubation conditions on the culturability of freshwater bacterioplankton. Appl. Environ. Microbiol. 69 (4), 1980–1989. Camazine, S., Deneubourg, J.-L., Franks, N., Sneyd, J., Theraulaz, G., Bonabeau, E., 2001. Self-Organization in Biological Systems. Princeton University Press, Princeton. Castets, V., Dulos, E., Boissonade, J., De Kepper, P., 1990. Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern. Phys. Rev. Lett. 64 (24), 2953–2956. Feudel, U., Jansen, W., 1992. CANDYS/QA—a software system for the qualitative analysis of nonlinear dynamical systems. Int. J. Bifurc. Chaos 2, 773–794.
94
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
Henry, B., Wearne, S., 2002. Existence of Turing instabilities in a two-species fractional reaction–diffusion system. SIAM J. Appl. Math. 62 (3), 870–887. Hunter, K., Wang, Y., VanCappellen, P., 1997. Kinetic modeling of microbially-driven redox chemistry of subsurface environments: coupling transport, microbial metabolism and geochemistry. J. Hydrol. 209, 53–80. Kitsunezaki, S., 1997. Interface dynamics for bacterial colony formation. J. Phys. Soc. Jpn. 66, 1544–1550. Kuznetsov, Y., Levitin, V., 1996. Content: a multiplatform environment for analyzing dynamical systems. Dynamical Systems Laboratory, Centrum voor Wiskunde en Informatica, Amsterdam. Kuznetsov, S., Mosekilde, E., Dewel, G., Borckmans, P., 1997. Absolute and convective instabilities in a one dimensional Brusselator flow modell. J. Chem. Phys. 106 (18), 7609–7616. Madani, S., Meysman, F., Middelburg, J., 2003. Biogeochemical modeling of sediments from the Santa Barbara Basin (California). In: Rullkötter, J. (Ed.), BioGeoChemistry of Tidal Flats, vol. 12. Forschungszentrum Terramare, Wilhelmshaven, pp. 91–93. Mudryk, Z., Podgórska, B., Ameryk, A., Bolałek, J., 2000. The occurrence and activity of sulfate-reducing bacteria in the bottom sediments of the gulf of Gda´nsk. Oceanologia 42 (1), 105–117. Murray, J., 1993. Mathematical Biology, second ed. Springer, Berlin. Nicolis, G., 1995. Introduction to Nonlinear Science. Cambridge University Press, Cambridge.
Nicolis, G., Prigogine, I., 1977. Self-Organization in Nonequilibrium Systems—From Dissipative Structures to Order through Fluctuations. Wiley, New York. Petrovskii, S., Li, B.-L., Malchow, H., 2003. Quantification of the spatial aspect of chaotic dynamics in biological and chemical systems. Bull. Math. Biol. 65 (3), 425– 446. Rovinsky, A., Menzinger, M., 1994. Differential flow instability in dynamical systems without an unstable (activator) subsystem. Phys. Rev. Lett. 72 (13), 2017–2020. Satnoianu, R., Merkin, J., Scott, S., 1998. Spatio-temporal structures in a differential flow reactor with cubic auto catalator kinetics. Physica D 124, 345–367. Satnoianu, R., Menzinger, M., Maini, P., 2000. Turing instabilities in general systems. Math. Biol. 41 (6), 493–512. Satnoianu, R., Maini, P., Menzinger, M., 2001. Parameter space analysis, pattern sensitivity and model comparison for Turing and stationary flow-distributed waves (FDS). Physica D 160, 79–102. Tsimring, L., Levine, H., Aranson, I., Ben-Jacob, E., Cohen, I., Shochet, O., Reynolds, W., 1995. Aggregation patterns in stressed bacteria. Phys. Rev. Lett. 75 (9), 1859– 1862. Turing, A., 1952. The chemical basis of morphogenesis. Philos. Trans. R. Soc. Lond. B 237 (B 641), 37–72. VanCappellen, P., Wang, Y., 1996. Cycling of iron and manganese in surface sediments: a general theory for the coupled transport and reaction of carbon oxygen, nitrogen, sulfur, iron and manganese. Am. J. Sci. 296, 197–243.
Turing patterns in a simple model of a nutrient–microorganism system in the sediment Martin Baurmann∗ , Ulrike Feudel ICBM, Carl von Ossietzky Universität, PF 2503, 26111 Oldenburg, Germany Received 2 September 2003; received in revised form 19 December 2003; accepted 11 January 2004
Abstract Sediments are characterized by heterogeneous distributions of nutrients and microorganisms which emerge as a result of the interaction between chemical and biological processes with physical transport. We study in a simplified model the dynamics of one population of microorganisms and its nutrient, taking into account that the considered bacteria possess an active as well as an inactive state. Furthermore, the nutrients are transported by bioirrigation. It is shown that under certain conditions Turing patterns can occur which yield heterogeneous vertical spatial distributions of species. Furthermore, this model exhibits several stable coexisting spatial profiles, so that it depends crucially on the initial condition which of the distributions will be realized. This phenomenon of multistability can still be observed when spatial profiles are externally imposed by considering a depth-dependent bioirrigation. © 2004 Elsevier B.V. All rights reserved. Keywords: Reaction–diffusion models; Turing instability; Pattern formation; Sediment
1. Introduction The study of pattern formation in spatio-temporal systems has continued to be a central issue to modelers since the fundamental paper of Turing (1952). He analyzed the spontaneous emergence of inhomogeneous distributions of chemical substances in reaction–diffusion systems (RDS) and their possible importance for morphogenetic processes. Based on his seminal paper about diffusion instability a broad theory of pattern formation in RDS has been developed (Nicolis and Prigogine, 1977; Nicolis, 1995; ∗
Corresponding author. Fax: +49-441-798-3404. E-mail address: [email protected] (M. Baurmann).
Murray, 1993; Petrovskii et al., 2003; Satnoianu et al., 2000). A necessary condition for a diffusion instability to occur is the difference in diffusion constants of two nonlinearly reacting chemicals (activator and inhibitor). Since this is impossible to realize in solutions it took more than 30 years before Turing’s theoretical findings have been verified in experiment (Castets et al., 1990). In the last decade, the analysis of pattern formation has been expanded by including advection as a transport process (Rovinsky and Menzinger, 1994; Kuznetsov et al., 1997; Satnoianu et al., 1998, 2001). In particular, it has been shown, that the instability of a homogeneous distribution of chemical substances can be achieved by the interaction of nonlinear, not necessarily autocatalytic, reaction of three
1476-945X/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ecocom.2004.01.001
78
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
chemical species with advective transport (Rovinsky and Menzinger, 1994). Recently, it has been shown by Henry and Wearne (2002) that anomalous diffusion of the activator gives rise to Turing-patterns even if the diffusion coefficients of activator and inhibitor are equal. Pattern formation has also been a matter of research in biology and ecology. In particular, the aggregation behavior of microorganisms was studied. This paper focuses on possible mechanisms for the emergence of spatial patterns, i.e. inhomogeneous distributions of bacteria and nutrients in marine sediments. We show that the interaction between growth processes of bacteria, fed upon nutrients, and diffusion can lead to spatial pattern formation in the sediment. Our approach yields a qualitative understanding of the appearance of inhomogeneous spatial distributions of bacteria and nutrients using a very simplified model which reflects some fundamental properties of bacteria populations. Although this work has been inspired by experimental studies which show that chemical substances like sulfate or DOC are inhomogeneously distributed in depth (Mudryk et al., 2000; Madani et al., 2003), our goal is not to reproduce certain experiments quantitatively, but to explore principles of emerging patterns in marine sediments. The model we develop in this paper describes a type of interaction, that is an integral part of the complex ensemble of interactions found in benthic systems (see VanCappellen and Wang, 1996; Hunter et al., 1997). In contrast to these complex models only one bacteria species and its nutrient is considered in our model system. As a fundamental aspect of the model we divide the population of bacteria in an active and a dormant part. Additionally, we assume that the activation of dormant bacteria is mediated by the exchange of certain signal molecules. This assumption is motivated by several findings reported in litterateur: Kitsunezaki (1997) showed that the surfaces of bacteria colonies can propagate with different velocities yielding a fingering of the moving fronts. This special behavior is the result of the interaction of nonlinear diffusion and an activation mechanism of the bacteria. Similarly, Tsimring et al. (1995) divided the considered bacteria population in a motile and nonmotile part. In their model chemotaxis is the basic mechanism that induces the formation of spatial patterns. Camazine et al. (2001) showed that in slime
molds the chemotaxis, induced by a chemoattractant (cyclic adenosine monophosphate, cAMP), can lead to pattern-formation. The same substance was found by Bruns et al. (2002, 2003) to effect a stimulation of the growth of dormant cells in special planktonic bacteria in fresh and sea water and the authors state, that signal molecules (like cAMP) may be of relevance for the growth of a broader spectrum of marine bacteria. The paper is organized as follows: in Section 2, we introduce a simple model considering a bacteria species in a sediment that feeds upon one chemical substrate. To identify some basic properties of this model, we restrict ourselves first to a local model, in which only the reaction part of the system is taken into account, while diffusion is neglected. We find two equilibrium states and discuss their stability properties with respect to the variation of several parameters. These equilibrium states correspond to homogeneous distributions of species. While the equilibrium characterized by the extinction of bacteria is always stable, the equilibrium corresponding to positive population densities can lose its stability either in a Hopf-bifurcation (hb) or in a turning point (tp). Section 3 is devoted to the formation of spatial structures on a vertical 1D domain. We show that under certain conditions Turing structures appear which correspond to an inhomogeneous distribution of nutrients and microorganisms. These inhomogeneities emerge due to a diffusion instability without any spatial external forcing. This phenomenon presumes that the nutrient is much more mobile than the bacteria. Additionally, a certain supply of nutrient and a limit of the bacteria’s mortality must be guaranteed. Finally, we consider a spatial external forcing by incorporating boundary conditions or a bioirrigation function which decreases exponentially with depth. The aim of this part of our study is to check, whether pattern formation due to diffusion instabilities is robust against the influence of an external force. Although the stationary profiles differ from those found for the model with no spatial forcing, the general situation does not change too much: again, we find the coexistence of different spatially heterogeneous profiles in depth indicating that external forces cannot prevent pattern formation in large parameter intervals.
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
2. The minimal sediment model 2.1. The local model and the transport processes This paper focuses on a simplified model of the interaction of microorganisms and a chemical substance in the sediment. The microorganisms in the sediment build up their biomass by degrading special chemical species, using them as a nutrient. Consequently both, the population density of microorganisms, more specific of bacteria, and the distribution of chemicals in the sediment, is influenced by this interaction. Many different microorganisms and chemical substances are involved in the biogeochemical processes in the sediment (see e.g. VanCappellen and Wang, 1996; Hunter et al., 1997). Instead of examining the whole multitude of different pathways, we concentrate on a simple predator–prey model that we deploy as a prototype that exhibits some of the typical behavior,
79
found in the complete system. Thus, we set up a minimal sediment model (MS-model) in which we restrain to observe only one chemical species (concentration: N (mol/m3 )) acting as a nutrient and a population of one kind of bacteria (density: B (kg biomass/m3 )) that feeds on the nutrient. In the MS-model, five different processes are taken into account for the interaction between bacteria and nutrient (see Fig. 1): we consider two sources of input of nutrient into the system. On one hand, we include the diffusive flux through the sediment’s surface into the system, on the other hand, chemicals can be transported in lower regions of the sediment by bioirrigation leading to a nonlocal transport term. The interaction between bacteria and the chemical species consists of a consumption of nutrient and the corresponding increase of biomass of the bacteria due to the degradation of nutrient. As a last factor we include the mortality of microorganisms.
Diffusion
Boundary Flux
Bioirrigation
Seawater
local Processes Sediment
bioirrigation
Nutrient
consumption
Fig. 1. Minimal sediment model.
growth
mortality
Bacteria active inactive part part
(de)activation
80
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
The effects of bioirrigation are modeled by a diffusive input of Fickian type. It is sufficient to parameˆ − N), terize this diffusive term by a certain flux Φ(N ˆ where N is the nutrient concentration in the sea water. Thus, the nonlocal transport is considered as a diffusive connection between the sediment in each depth with the sea water above the sediment. In the model, this process is realized as a lateral boundary condition and, hence, a local process at each depth level. When examining the spatial model we will replace the fixed parameter Φ by a function Φ(z), where z is the depth below the sediment surface. By this means a decreasing effect of bioirrigation with increasing depth can be incorporated into the model. The flux of nutrient through the sediment surface is another boundary condition that has to be applied to the model’s uppermost point. It will be considered in a spatial model, whereas we will neglect it in the local formulation. To improve the consistency of the model with natural systems, a certain part of the bacteria population will be considered as inactive. Only the active part of the bacteria Ba feeds on the nutrient leading to a certain growth rate and dies with a mortality m, while the inactive part of bacteria may die but will not effect the nutrient consumption. The amount of active bacteria (Ba ) is not bound by a rigid factor but is determined by a (dimensionless) function r(B), which depends only on the biomass of microorganisms. The actual form of r is difficult to specify. There exist no well established models for the activation-mechanism of bacteria in the sediment. Nevertheless, we can state, that the biomass of active bacteria cannot exceed the total population density. So, r(B) is restricted to: 0 ≤ r(B) ≤ 1
(1)
Additionally, we expect an increase of r if B increases, presuming some activation mechanism induced by the organisms of B. At least for certain ranges of biomass B we have to guarantee: ∂r >0 (2) ∂B By that approach we model an activation mechanism that agrees with the results found by Bruns et al. (2002, 2003) in experiments with fresh and sea water bacteria. Realizing (1) and (2), we employ the function: B r(B) = (3) B+L
That means we assume this function to be of Monod-type. This type of function is one of the simplest mathematical functions possessing a saturation characteristics. In principle, any other function for which (1) and (2) holds is possible to choose and would give the same results. Applying this formula, the system is governed by the following equations: dB = αBa N − mB dτ dN ˆ − N) − βBa N = Φ(N dτ B Ba = r(B)B = B B+L
(4)
where α is the rate of conversion into biomass, β is the capturing rate, m is the mortality rate and L is the ˆ parameterize the exhalf saturation constant. Φ and N ˆ is the concentration of change due to bioirrigation: N N in the external pool of sea water and Φ parameterizes the hydraulic conductivity between sediment and sea water. All listed parameters have positive values. To reduce the number of model parameters and to make the Eq. (4) dimensionless we rescale them. Using m Y, α Φ L = K, β
N=
m Yˆ , α 1 τ= t Φ
ˆ = N
B=
Φ X, β
our resulting system of equations reads: m X dX = XY − X dt Φ X+K X dY = Yˆ − Y − XY dt X+K
(5) (6)
This is the dimensionless three parameter formulation (m/Φ can be regarded as one independent parameter) of the local MS-model, which shall be examined in the following. Though the consumption of nutrient and growth of bacteria’s population are modeled by a Lotka–Volterra approach, both are highly non-linear with respect to X, because the feeding term depends on the biomass of active bacteria Ba only. Therefore, the MS-model can be considered as a predator–prey system with a non linear predator-dependent response function.
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
2.2. Dynamical properties of the local model The analysis of the dynamical properties of the local model is a necessary assumption to understand the possible states of the homogeneous system without local transport. The system is in equilibrium when the variables do not change in time. It has then reached a stationary state. The locations and characteristics of these states are of high interest for us, because they give answers, how the system will behave in the long time run. Particularly we are interested in the question for which values of the parameters the population of bacteria will be able to survive. Analyzing the model equations, we find a trivial stationary state at [X0∗ = 0, Y0∗ = Yˆ ], at which the bacteria are extinct and the nutrient has reached the concentration of the sea water (Yˆ ). There might exist two other stationary states in which nutrient and microorganisms coexist at ∗ , Y ∗ ), where (X1,2 1,2 ∗ X1,2 = 21 (Yˆ − 1 ± (Yˆ − 1)2 − 4K) (7) and Yi∗ =
Xi∗ + K Xi∗
i = 1, 2
(8)
Regarding (7), the local model will show these additional stationary states, if (Yˆ − 1)2 − 4K ≥ 0
(9)
Inequality (9) indicates the ranges where different stationary states can coexist for a given set of parameters: √ K ≤ 41 (Yˆ − 1)2 or Yˆ ≥ 1 + 2 K The parameters m and Φ have no impact on the ranges of coexistence. Since we examine a natural system, it is fundamental, to check, whether stationary states are stable with respect to small perturbations or not. If they are not, even small changes in the nutrient’s concentration or the bacteria’s biomass will destabilize the equilibrium and the system’s variables will converge to a stable stationary state. Due to the inherent fluctuations of natural systems unstable stationary states will only exist in a mathematical context and cannot be observed by measurements or experiments. Nevertheless, they
81
might be of importance for the dynamics of the system as we will discuss later. At this point, we will not reproduce the whole analysis, that can be used to determine a stationary state’s stability; it is described in detail in a large number of publications (e.g. Murray, 1993; Nicolis and Prigogine, 1977; Arrowsmith and Place, 1992). The stability of a stationary state is determined by the signs of the real parts of the eigenvalues computed from the Jacobian matrix at the steady state. If the real parts of all eigenvalues are negative, then the steady state is stable, otherwise not. For the MS-model we obtain two eigenvalues: √ m (T ± T 2 − 4S) (10) λ1,2 = Φ 2(Xi∗ + K) with Φ Yˆ Xi∗ , S = (Yˆ − 1)Xi∗ − 2K (11) m According to the number of negative real parts of the eigenvalues we can classify the stationary state. If both eigenvalues have negative real parts the steady state is stable and called attractor. If at least one real part of the eigenvalues is positive, the steady state is unstable and called a saddle (one positive and one negative real part) or repellor (two positive real parts). Let us first discuss the stability of (X2∗ , Y2∗ ). In this case we obtain S < 0. Hence, for the steady state (X2∗ , Y2∗ ) both eigenvalues of the Jacobian are real and have different signs, so that (X2∗ , Y2∗ ) is a saddle. Contrary, at (X1∗ , Y1∗ ) we have S > 0. The real parts of both eigenvalues have the same sign, which is just the negative of the sign of T . Consequently we can distinguish: T =K−
if X1∗ >
mK Φ Yˆ
[X1∗ , Y1∗ ] is an attractor
if X1∗ <
mK Φ Yˆ
[X1∗ , Y1∗ ] is a repellor
We omit substituting the merely complex algebraic expression for X1∗ . Instead we sketch the two qualitatively different stability properties instead, in order to illustrate all possible transitions of stationary states. For that purpose, we will vary the parameter Yˆ , corresponding to the concentration of nutrients in the sea water, whereas for the other parameters we choose two different sets of constants listed in Table 1.
82
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
Table 1 Parameter sets for the MS-model
Parameter set 1 Parameter set 2
X*
K
m
Φ
0.8 8.0
0.5 8.0
1.0 1.0
14 12 10 8 6 4
These sets differ in the values of K and m. m is the mortality rate, so that for parameter set 2 the microorganisms are more likely to die. The parameter K affects the activation ratio of the bacteria population. If the value of K is relatively small (parameter set 1), a slight increase of bacteria’s biomass causes an almost total activation of microorganisms, whereas a higher value of K (parameter set 2) results in a slower activation. Generally, the ratio of active bacteria is (at same amounts of biomass) higher, if K is smaller. Fig. 2 illustrates the situation for the two different parameter sets depending on the concentration of nutrients in the sea water Yˆ . We can distinguish two kinds of qualitative behavior corresponding to different environmental conditions. Let us first discuss the stationary states and their stability. For both parameter sets the trivial steady state corresponding to an extinct bacteria population (X0∗ , Y0∗ ) is stable in the whole Yˆ -interval under consideration. But for a part of this interval there exist the two other ∗ steady-states X1,2 which are created in a turning √ point at Yˆ tp = 1 + K. X1∗ and X2∗ exhibit different stability properties. For parameter set 1 (Fig. 2(a)) the turning point corresponds to a saddle-node bifurcation leading to the emergence of a stable attracting equilibrium state (X1∗ , Y1∗ ) and an unstable saddle (X2∗ , Y2∗ ). Thus, the system possesses two coexisting stable equilibria beyond the turning point (Yˆ > Yˆ tp ), i.e. the system is bistable. In this case it depends crucially on the initial conditions which of the two stable states (X0∗ , Y0∗ ) or (X1∗ , Y1∗ ), will be realized by the system in the long-term limit. Starting with a nonzero bacteria population we can change the environmental conditions in such a way that we decrease the nutrient concentration in the sea water. When reaching the turning point, we observe a sudden collapse of the bacteria population due to the loss of stability. As a result the bacteria population becomes extinct.
2 0 (a) 0
2
4
6
8
^ 10 12 14 Y
2
4
6
8
^ 10 12 14 Y
X* 14 12 10 8 6 4 2 0 (b) 0
Fig. 2. Stability of stationary states for two different parameterizations and variable value Yˆ (axis of abscissae). Stable states are marked by thick solid lines, unstable by thin dashed lines. The box symbolizes a turning point, the star a Hopf-bifurcation. From the Hopf-bifurcation point a steep dotted line emerges, signifing a set of periodical points on a unstable orbit.
For parameter set 2 we observe bistability only for ∗ exist. Stability a part of the Yˆ -interval, where X1,2 of a nonzero bacteria population occurs only for high values of the nutrient concentration in the sea water Yˆ . Decreasing Yˆ leads again to a collapse of the bacteria population, but in a different dynamical scenario. The stationary state loses stability in a Hopf-bifurcation instead of a saddle-node bifurcation. The Hopf-bifurcation leads to the emergence of periodic behavior (dotted line in Fig. 2b). However, these oscillations cannot be observed in a natural system because they are unstable (subcritical Hopf-bifurcation). A Hopf-bifurcation occurs if a pair of complex conjugate eigenvalues of the Jacobian crosses the imaginary axis. Presuming m 2Yˆ > Φ Yˆ − 1
(12)
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
83
Fig. 3. The parameter space of the MS-model. Parameter sets corresponding to a turning point situation are symbolized by the bright surface, those connected with a Hopf-bifurcation lie on the dark surface. The thick black line is the line where both surfaces touch each other. It is the set of Takens Bogdanov points.
this can be obtained by setting Yˆ Φ Yˆ Φ Yˆ − 1 − K = Khb := m m
(13)
Fig. 3 illustrates the different kinds of dynamical behavior in the three dimensional parameter space spanned by (K, m/Φ, Yˆ ) The set of turning points and Hopf-bifurcations are the surfaces Ktp and Khb , respectively. These surfaces divide the parameter space into three subspaces. Parameter sets corresponding to points which lie above the Ktp -surface designate pa∗ , Y ∗ ) do not exist. If (12) rameter sets where (X1,2 1,2 holds and Khb < K < Ktp (space between the surfaces) the non-trivial stationary state exists and is a repellor. Only the points below both surfaces stand for parameter sets that generate bistable dynamics, because besides of the trivial equilibrium (0, Yˆ ) also (X1∗ , Y1∗ ) is an attractor. Since we have to fulfil (12) the Khb -surface is bounded. At its boundaries it touches the Ktp -surface and is not continued across. Both surfaces touch in a line, that consists of so called Taken Bogdanov points. They symbolize parameter sets, showing eigenvalues λ1 = λ2 = 0 at (X1∗ , Y1∗ ) = (X2∗ , Y2∗ ) = √ √ (2 K, 2 + K).
Fig. 3 shows that there are two fundamental ways, how the steady state of nonzero bacteria’s biomass may loose its stability, when environmental parameters change. The first one occurs by crossing a turning point (tp), the second one coincides with traversing of a subcritical Hopf-bifurcation. Both diagrams in Fig. 2 correspond to different lines in Fig. 3. Let us discuss the distinct ways of destabilization also from another point of view. As already mentioned above bistability means that the final state approached by the system in the long-term limit depends on the initial condition. To study the robustness of the system with respect to larger perturbations, it is useful to analyze the basins of attraction of the different coexisting stable states. The basin of attraction is the set of initial conditions which converges to one stable state of the system. In our example we have two basins of attraction for the stable nonzero state (X1∗ , Y1∗ ) and the extinct state (X0∗ , Y0∗ ), respectively. Fig. 4 shows the two stable stationary states and their basins of attraction for values of Yˆ close to the critical values for the loss of stability (Yˆ ≈ 2.78 for tp and Yˆ ≈ 9.1 for hb). In the hb-type system, the basin boundary is made up by the unstable periodic orbit emerging at the Hopf-bifurcation, while in the tp-type
84
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
Fig. 4. Steady states and their basins of attraction; series A: Parameter set 1 and Yˆ = [3.0, 2.9, 2.8] (tp-type); series B: Parameter set 2 and Yˆ = [9.5, 9.25, 9.15] (hb-type). Stable equilibira are marked by a +, unstable ones by a 䊊.
the basin boundary is connected with the unstable equilibrium. In the latter case, the role of the unstable equilibrium becomes evident. Even though it is not observable in nature, its location is important for the dynamics of the system. The unstable steady-state
(saddle point) “sits on” the boundary of the basin of attraction, its stable manifolds make up the basin boundaries. Both diagrams illustrate that the states are moving closer to their basin boundaries, while Yˆ decreases.
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
The difference between the two examined ways is, that the size of hb-type’s basin converges to zero, while the tp-type’s basin collapses abruptly. According to that mathematical finding, natural systems representing the hb-type are more likely to lose stability, than their tp-type counterparts, since a fluctuation to any direction will cause the bacteria to become extinct, if it is only large enough. As a consequence hb-type-systems with parameters close to the critical bounds can be destabilized by reasonably small perturbation even in the case, when they are still stable.
3. The spatial minimal sedimentmodel 3.1. Setup of the spatial model In the sediment, the local processes analyzed in Section 2 are influenced by transport processes like diffusion leading to spatial distributions of chemicals and bacteria populations. Therefore, we now consider a spatially extended system in order to include the effects of diffusion. In particular, we are interested in the changes of the stability properties of homogeneous species distributions due to diffusion. We apply the model on a vertical one-dimensional domain with space coordinate ζ (later with dimensionless coordinate z). The diffusion process takes place along the vertical direction and is modeled by Fickian law. At the lower and upper boundaries we have to establish boundary conditions: we assume that we have no fluxes over the lower boundary and accept a certain flux of nutrient at the sediments surface, which forms the upper boundary. As another step towards a realistic model we will take into account, that effects of bioirrigation will decrease with increasing depth. To incorporate this phenomenon, we exchange the fixed parameter Φ by a monotonously decreasing function ζ Φ(ζ) = Φ exp −a (14) ζl where a is a non-negative factor parameterizing the exponential function and ζl is the length of the domain. Inducing these changes in the original equations (formulated with respect to B and N) and scaling (z =
85
ζ/ζl )leads to: m ∂2 X XY ∂X = X − 1 + DX 2 ∂t Φ X+K ∂z ∂Y X ∂2 Y = exp(−az)(Yˆ − Y) − XY + DY 2 ∂t X+K ∂z (15) and ∂Y = Φsurf (Yˆ − Y) ∂z z=0
∂Y =0 ∂z z=1
(16)
with the conductivity across the surface Φsurf and the dimensionless diffusion coefficients: DB DN DX = and DY = (17) 2 Φζl Φζl2 To simplify the differential equations we will, excluding the case DX = 0, implement δ with δ = DY /DX . This substitution is helpful, because the ratio of the diffusion coefficients is of more interest than their absolute values. √ In case of constant bioirrigation a scaling by z = DX ξ leads to: ∂X m XY ∂2 X = X −1 + 2 ∂t Φ X+K ∂ξ (18) ∂Y X ∂2 Y = Yˆ − Y − XY + δ 2 ∂t X+K ∂ξ 3.2. Diffusion-induced patterns in the spatial model In this section, we want to examine the effects of diffusion only, so we set Φsurf = 0 and a = 0. In that case the conditions of nutrients to flow into the system are equal in all points of the domain. Since diffusive fluxes are driven by gradients, they only affect the solution, if gradients occur. As a consequence, all homogeneous profiles referring to the local model’s solutions will also solve (18). Thus, one solution of (18) can always be found: these are the trivial homogeneous functions (X = 0, Y = Yˆ ) for all values z. As we have shown above an additional non-trivial, homogeneous solution may exist in particular intervals of parameters. Besides these two homogeneous solutions no other profile without concentration gradients exists. Again we look at the stability of the homogeneous depth profiles of the concentration with respect to perturbations. We will show in the following, that the
86
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
non-trivial homogeneous stationary state may loose its stability, if diffusion is taken into account as an additional process. Instead non-homogeneous concentration profiles will arise. At first sight this possibility seems to be absurd: diffusion is known as a process, that will extinguish spatial gradients. However, it has been shown by Turing (1952) that nonlinear chemical reactions in spatially extended systems can lead to an inhomogeneous distribution if diffusion is taken into account. These Turing instabilities are excited by perturbations of specific spatial wavelengths and will only appear within certain parameter spaces (Alonso et al., 2002). Our aim is to identify these Turing-spaces with respect to parameters and wavelengths and to discuss associated system properties. With respect to the spatial MS-model we study the question: under which conditions is it possible, that homogeneous profiles (X(z) = X1∗ , Y(z) = Y1∗ ) become unstable, presuming (X1∗ , Y1∗ ) is an attractor in the local model? To answer this question we analyze the dynamics of the spatio-temporal system (15) with respect to a perturbation modeled with a wave ansatz (compare Murray, 1993). We consider perturbation functions x = ζx exp(λt − x), i2πκz = ζy exp(λt − x) and investigate their time evolution in the linearized model. Taking into account the chosen boundary conditions, we can restrict ourselves to the following class of functions: x(z, t) = cx cos (2πκz) exp(λt)
1/2 j11 DY +j22 DX
− (j11 DY +j22 DX )2 −4DX DY det(J) 1 κˆ − = 2π 2DX DY (21) and
1/2 j11 DY + j22 DX
+ (j11 DY +j22 DX )2 −4DX DY det(J) 1 κˆ + = 2π 2DX DY (22) where J denotes the Jacobian matrix of the local model at (X1∗ , Y1∗ ) and its elements ˙ ∂X j11 = ∂X X∗ ,Y ∗ 1
and j22
1
∂Y˙ = . ∂Y X∗ ,Y ∗ 1
1
To avoid the complex terms of κˆ − and κˆ + , we first check that j22 < 0 and gain a weaker formulation of (20) by the estimate: Re(λ(κ)) > 0 ⇒ κ− < κ < κ+
y(z, t) = cy cos (2πκz) exp(λt)
with
Linear stability analysis yields intervals for the wavelength κ, within which the homogeneous solution is unstable with respect to perturbations with this wavelength. Again the eigenvalues of the corresponding Jacobian which depend now on the wavelength κ and the diffusion constants DX and DY have to be positive. For reaction–diffusion systems like that of the MS-model, we have (compare Petrovskii et al., 2003): Re(λ(κ)) < 0 if η := (j11 DY + j22 DX )2 − 4DX DY det(J) < 0 (19) If η > 0 it follows: Re(λ(κ)) > 0 ⇔ κˆ − < κ < κˆ +
with
(20)
1 κ− = 2π
(23)
det(J) , DY j11
1 κ+ = 2π
j11 DX
(24)
Since the considered equilibrium is an attractor in the local formulation, we have det(J) > 0. Hence, we need j11 > 0 to fulfil (24). This means that the increase of the bacteria’s biomass has to be autocatalytic. We obtain this special property in the MS-model because the activation of bacteria is more pronounced for large populations, since we assume that bacteria are capable to activate dormant cells by excreting signal molecules. κ− < κ < κ+ is a necessary condition to obtain Turing instabilities. By this inequality we can outline a space in the κ −DY -plane, where Turing instabilities
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
87
may be obtained. Outside this space, no such instabilities occur. In particular, we have no Turing instabilities, if κ− > κ+ . Substituting (24) we obtain: det(J) (25) δ < 2 ⇒ no Turing instabilities j11
Condition (26) can be formulated stricter: by derivating Re(λ(κ)) with respect to κ, we find a maximum at 1 DX j22 + DY j11 κe = (27) 2π 2DX DY
So for the MS-model we can state:
Because Re(λ(κ)) ≤ Re(λ(κe )) it is obvious that there is a Turing space only if Re(λ(κe )) > 0 demanding as a necessary condition:
Φ(X1∗ + K)[(Yˆ − 1)X1∗ − 2K)] δ< mK2 ⇒ no Turing instabilities
(26)
Eq. (26) shows, that the increase of Yˆ , Φ and DX leads to stabilization. On the other hand, increasing K, m and DY destabilizes the system. Altogether we can state that the occurrence of Turing instabilities is enhanced by a low inflow of nutrient into the system and a high loss of bacteria due to mortality. 12 κ
+
10
κ
8
κe
6 4 2κ 0 −5
−4
(a) 10
14 12
10
−3
10
Dx
κ
10 κ
8 6 4κ + 2 0 −2
(b) 10
κe
10
−1
0
10
DY
Fig. 5. κ− , κ+ and κe for different values of DX and DY , respectively. Outside the area limited by the thick lines, no Turing instabilities can be met.
κ− < κe < κ+
(28)
Since we can show that the inequality κe < κ+ holds on the entire parameter space, the essential condition is κ− < κe (see Fig. 5) leading to: δ<
det(J) − j12 j21 ⇒ no Turing instabilities 2 j11
(29)
3.3. Simulations on a discretized domain The analytical examination of the sediment model led to the important result, that an attracting homogeneous state of nonzero bacteria population may loose its stability, if diffusion is taken into account. Beyond the diffusion instability inhomogeneous distributions can occur which cannot be computed analytically. To illustrate the emergence of inhomogeneous distributions of bacteria and nutrient, we consider 41 boxes. We use the parameter set 1 and set Yˆ = 8, DX = 2 × 10−4 and δ = 5000. The initial values are the homogeneous profiles that are perturbed by a cosine function with five periods over the domain length (pattern type 0). Fig. 6 shows, that for parameter set 1 the initial spatial frequency of the profile’s perturbation is changing quickly in time and finally a pattern will form, consisting of two cosine cycles (pattern type 1). The transient dynamics is rather short, so that already after time t = 60 the final distribution exists. To illustrate how a change of bacteria’s mobility will affect the pattern formation process, we will change now the value of δ. For a merely small interval of δ ∈ [0, 5000] a homogeneous profile will appear in the longtime-run (see Fig. 7). This spatially constant profile will loose its stability when δ increases. An abrupt transition to two maxima occurs after crossing the Turing instability threshold. Choosing higher values 5000 < δ < 6000,
88
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
Fig. 6. Long-term behavior of a spatial cosine oscillation in the MS-model (parameter set 1, Yˆ = 8, DX = 2 × 10−4 , δ = 5000).
the amplitudes of the profile are growing rapidly until they reach a range in which the increase is moderate 6000 < δ < 13,000 (pattern type 1). At a value of δ = 13,000, the qualitative behavior of the MS-model’s long-term-solution changes abruptly again: instead of the spatial oscillation with two cosines, a two-modal profile occurs (pattern type 2). Fig. 8 illustrates the abrupt change of the stationary states. Additionally, it is interesting to note, that the small increase of δ causes the profile to show a short spatial interval, where X is extinct. This interval widens, if δ is increased further. The consequence of the appearance of this “dead-zone” in the middle of the domain is that at this point, the domain is separated into two subdomains with no diffusive mass transfer between them. This would mean that life of microorganisms in deeper layers of the sediment may be maintained only by nutrient inflow due to bioirrigation. Further analysis shows that the two different types of inhomogeneous profiles, shown in Fig. 8 coexist
for a broad range of δ-values. Thus, we can conclude, that the transition from pattern type 1 to pattern type 2 at δ = 13,000 in Fig. 8 is not due to a bifurcation but rather to a transition of the initial condition of our simulation from one basin of attraction to another one. To study, which stationary states exist for the MS-model and how they change, when the parameter values are modified, we use the program Content (Kuznetsov and Levitin, 1996) to follow equilibrium points under variation of DY . Content allows the analysis of one-dimensional PDEs at a user defined discretization of the spatial domain. We use parameter set 1 and chose a discretization (1000 boxes) fine enough that our results are stable with respect to the discretization length and the detection of equilibrium points that are numerical artifacts is avoided. Fig. 9 shows a diagram in which the biomass of bacteria X in the uppermost discretization-box is plotted versus DY . It illustrates the emergence of stationary states from the homogeneous profile indicated by the
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
89
Fig. 7. The increase of δ effects (abrupt) changes of the equilibrium states of bacteria population (subplots only differ by the location of viewpoints).
90
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
20
X
15
10
5
0
0
0,2
0,4
0,6
0,8
1
depth Fig. 8. Stable profiles at δ = 13,000 (thin line) and δ = 13,001 (thick line).
horizontal line at X(z = 0) = 6.88, Y(z = 0) = 1.14. At each bifurcation point on that line (marked by squares) a non-homogeneous profile emerges. Since these profiles are characterized by a certain wavenumber (at least in the neighborhood of the bifurcation) we can assign this characteristic wavenumber to the bifurcation points. Knowing that Re(λ(κ)) = 0 at the bifurcation points we can evaluate the location of the bifurcation points using (20) and deriving the critical values DYc (κ). We obtain: κ 1/2 1 3/2 2 5/2 ≥3 DYc () 6.1614 1.7712 1.0370 1.0076 5.0632 <0
Fig. 9. Projection of the MS-model bifurcation diagram on the DY /X[0]-plane. Squares symbolize bifurcations (not all plotted). Most bifurcation branches are labeled due to their main oscillation mode (number) and their orientation at the origin (sign). The branch of homogeneous profiles we marked by 0, and branches with extincted population of bacteria in the middle part by exc. There exist also a pair of branches representing unimodal profiles with bacteria extinction at the boundaries (exb). They superpose the 2- and 3/2- branches displayed very close to the abscissa.
At DY ≈ 1.0, we find two bifurcation points in close succession. They can be assigned to the wavenumbers 2 and 3/2, respectively. At DY ≈ 1.77, profiles of wavenumber 1 emerge, followed by a bifurcation with κ = 5/2 at DY ≈ 5.06. Profiles showing an half cosine oscillation κ = 1/2 emerge at DY ≈ 6.16. Note that solutions on the κ = 1 branch do not exist for DY > 3.6 . . . , since the branch “disappears” in a bifurcation with the κ = 2 branch. In the same way, the branches characterized by κ = 1/2 meet in a bifurcation on the 2+ set at DY ≈ 60 (not shown).
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
All other branches can be followed to arbitrary high values of DY . At [DY ≈ 2.4, X ≈ 6.88], the projection shows a superposition of homogeneous profiles and another curve exc. This second curve represents profiles with zones of extinction (as shown in Fig. 7). There are no equilibrium points on this curve showing homogeneous distributions of microorganisms. At [DY ≈ 2.4, X ≈ 6.88] we find two distinct equilibrium points (one homogeneous, the other not) that are superimposed by projection. A similar branch of inhomogeneous solutions representing profiles with extinction of bacteria at the boundaries and a unimodal distribution of biomass in the middle can also be followed, but can hardly be distinguished from the abscissa in the projection plot. Additionally to the location of equilibrium profiles their stability with respect to perturbations is of interest. Unfortunately, Content does not provide the necessary eigenvalues. Nevertheless, we were able to get some information on stability by employing the program Candys/QA (applying a much coarser spatial discretization) (Feudel and Jansen, 1992). We could derive a similar bifurcation diagram but got additional branches, which we consider as numerical artifacts (due to the coarse discretization). Furthermore, we check the stability for some values of DY and some wavenumbers κ by numerical experiments: we start with the homogeneous profile and perturb it with a small oscillation (wavenumber κpert ). By forward integration we test to which stationary state the perturbed profile converges. Table 2 summarizes the results of our numerical experiments. We see, that the homogeneous profile is loosing stability in favor of profiles with two cosine oscillations (2+ and the phase shifted 2−). Increasing DY further this attractor coexists with stable 3/2+ and 3/2− profiles. At higher values of DY , profiles with a zone of extinction (exc+ and exc-) will become attractive, instead of the 2+ and 2− distributions (cf. Fig. 7). On the 1/2±, 1± or 5/2±
91
Fig. 10. Location of stationary states in a MS-model with nonuniform bioirrigation. The 3/2+ branch (compare Fig. 9) could not be found.
branches, we found no attracting intervals, they seems to be completely unstable. 3.4. Effects of non-uniform bioirrigation So far we discussed a model with uniform bioirrigation acting as a spatially homogeneous forcing. With respect to natural systems this is not realistic. As mentioned above, the input of nutrient—due to bioirrigation—becomes smaller with increasing depth. This induces an inhomogeneous external forcing and thus we expect an imposed external depth-dependence. Let us check whether this externally imposed forcing destroys the pattern formation process. To take the depth dependence into account we have to analyze the Eq. (15) with a > 0. Fig. 10 shows the bifurcation diagram for a model with a = 0.1. So, the effective hydraulic conductivity at the top of the domain is Φ(z = 0) = Φ whereas it is Φ(z = 1) = exp(−0.1)Φ = 0.9Φ at the lower boundary. Of course this is a very moderate decrease
Table 2 Numerical experiments: the homogeneous profiles were perturbed by a (cosine) oscillation of wavenumber κpert and converge to the noted branch in the bifurcation diagram (cf. Fig. 9) DY κpert Converges to branch
1.3 3/2 3/2+
1.5 2 2+
2.0 2 2+
4.0 2 2+
4.0 3/2 3/2+
6.0 2 exc+
6.0 3/2 3/2+
10.0 2 exc+
10.0 3/2 3/2+
10.0 5/2 exc+
92
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
of bioirrigation effects, but it allows us to compare the diagram with Fig. 9. Surveying both diagrams we are able to identify most of the branches (except of a part of the 3/2+ branch, that is missing in Fig. 10). In contrast to the model with uniform forcing the diagram displayed in Fig. 10 shows no bifurcations connecting branches. Pitchfork bifurcations split up and lead to the disconnection of branches. The reason for this splitting is a symmetry breaking, which is introduced by the depth-dependent bioirrigation. Additionally, there is a shift of the branches to the left, so that we can say that coexistence is more easy to achieve (i.e. with a smaller DY ). This observation is consistent with our findings that the decrease of inflow of nutrient leads to a weakened stability of the homogeneous profile. The depth-dependence of bioirrigation applies already a certain spatial structure in form of an external forcing. Surprisingly, the emerging spatial distribution of species beyond the instability do not show only the imposed spatial profile but a much broader variety of different spatial distributions. Thus, depthdependent external forcing do not prevent the emergence of spatial patterns different from the imposed ones. 3.5. Effects of boundary fluxes Very similar to Fig. 10 is the bifurcation diagram we obtain by running the model with uniform bioirrigation but including an upper boundary flux of nutrient (Fig. 11). Following the same strategy as in the last subsection, we limit the effects of the changed boundary conditions very strictly to allow the comparison with the bifurcation plot in the uniform case (Fig. 9). Fig. 11 shows some common features with Fig. 10. In particular, we again obtain symmetry breaking and hence a splitting up of bifurcation points. The shift of branches towards lower values (in comparison to the uniform case) is neglectible, although we would expect a shift to the right due to the stabilizing effects of higher inflow of nutrients. An important change occurs looking at the spatial profiles with a “dead-zone”: the formerly isolated exc+/ exc- branch has now contact to other branches (of type 3/2). So, when changing parameters a “dead zone” may appear without a sudden change of the system’s behavior but by the expansion of a zero point.
Fig. 11. Bifurcation diagram for the spatial MS-model with uniform bioirrigation and a small upper boundary-flux Φsurf = 0.01.
4. Discussion We have studied a minimal sediment model, describing the interaction of a bacteria population and its nutrient. In this model we take into account that bacteria can live in an active state, in which they reproduce consuming nutrient, as well as in a passive state of dormancy. Furthermore, we assume that the activation of bacteria is enhanced by the size of population. This assumption is supported by experimental studies which show a stimulation of the growth of certain marine bacteria by signal molecules. The local dynamics of this model, which can be regarded as a two-species predator–prey system (with input of prey, predation and the mortality of the predator) is characterized by the coexistence of two equilibria over a large range of parameters. One equilibrium is a state, in which the bacteria population is extinct. This (trivial) equilibrium is stable in the whole parameter range considered. Assuming that the activation of dormant bacteria is fast and the input of nutrient is high, there is another stable equilibrium, characterized by a positive population density of bacteria. Varying the parameters of the system this steady state can loose its stability by crossing a Hopf-bifurcation or a turning point. In the first case, the basin of attraction of the equilibrium is contracting to the equilibrium point, while in the second case it collapses abruptly. Thus, one can conclude that the system is more sensitive to perturbations in parameter regions in which
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
the Hopf-bifurcation is close. The reason is that the Hopf-bifurcation is subcritical leading to an unstable limit cycle whose location determines the boundary of the basin of attraction for the nontrivial equilibrium state. The main focus of this paper is the study of pattern formation processes due to the interaction of nonlinear growth processes of the bacteria and diffusive transport. We have shown that in our simple minimal sediment model Turing patterns can emerge. Neglecting diffusion and nonconstant forcing the only spatial profiles, that are stationary and stable, are those, that correspond to the equilibria, found for the local model. When diffusion comes into play, the non-trivial equilibrium can loose its stability in favor of a heterogeneous distribution of species due to a diffusion instability. We have derived the conditions for the instability to occur: the nutrient has to be much more mobile than the bacteria which is usually fulfilled in the sediment due to the stickiness of bacteria to the sediment matrix. Additionally, a minimal supply of nutrient must be guaranteed and the mortality of bacteria must be relatively low. A further analysis of the model proved, that several stable heterogeneous profiles can coexist. In that case it depends on the perturbation that is added to the constant profiles to which the model converges. We depicted the location of equilibria in the state-space in a bifurcation diagram (see Fig. 9). Further analyses of the spatial model yield the result that the minimal sediment model exhibits multistability, i.e. it possesses several stable coexisting heterogeneous profiles. In that case it depends crucially on the initial condition to which heterogeneous distribution of nutrient and bacteria the model converges. Thus, there are different spatial profiles possible under the same environmental conditions. Finally, we have studied the influence of a nonconstant spatial forcing on the pattern formation process. This forcing describes decreasing bioirrigation with increasing depth and/or boundary fluxes. Both situations are more realistic in natural environments than constant forcings. Surprisingly, the externally imposed spatial structures did not destroy the pattern formation process. We again find multiple coexisting heterogeneous profiles different from the one induced by the non-constant forcing. Qualitative changes in
93
the bifurcation diagram are only observed when the non-uniformity of the forcing is pronounced. It is important to note that the non-constant forcing leads to a symmetry-breaking which results in a split up of bifurcation points. Furthermore, the original homogeneous profile is replaced by an inhomogeneous one. According to our studies the appearance of Turing patterns can be regarded as one possible mechanism of the formation of heterogeneous spatial distributions of nutrients and microorganisms in the sediment. But there are many other processes which may influence the emergence of spatial patterns in the sediment and which are not taken into account here. Such processes include advection, chemotaxis and the competition between different degradation pathways. Future work needs to take these processes into account.
Acknowledgements We thank H. Cypionka, W. Ebenhöh and T. Gross for discussions. This work was supported by the DFG research group BioGeoChemistry of Tidal Flats.
References Alonso, D., Bartumeus, F., Catalan, J., 2002. Mutual interference between predators can give rise to Turing spatial patterns. Ecology 83 (1), 28–34. Arrowsmith, D., Place, C., 1992. Dynamical Systems—Differential Equations, Maps and Chaotic Behavior. Chapman & Hall/CRC Press, London. Bruns, U., Cypionka, H., Overmann, J., 2002. Cyclic AMP and acryl homoserine lactones increase the cultivation efficiency of heterotrophic bacteria from the central baltic sea. Appl. Environ. Microbiol. 68 (8), 3978–3987. Bruns, A., Nübel, U., Cypionka, H., Overmann, J., 2003. Effect of signal compounds and incubation conditions on the culturability of freshwater bacterioplankton. Appl. Environ. Microbiol. 69 (4), 1980–1989. Camazine, S., Deneubourg, J.-L., Franks, N., Sneyd, J., Theraulaz, G., Bonabeau, E., 2001. Self-Organization in Biological Systems. Princeton University Press, Princeton. Castets, V., Dulos, E., Boissonade, J., De Kepper, P., 1990. Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern. Phys. Rev. Lett. 64 (24), 2953–2956. Feudel, U., Jansen, W., 1992. CANDYS/QA—a software system for the qualitative analysis of nonlinear dynamical systems. Int. J. Bifurc. Chaos 2, 773–794.
94
M. Baurmann, U. Feudel / Ecological Complexity 1 (2004) 77–94
Henry, B., Wearne, S., 2002. Existence of Turing instabilities in a two-species fractional reaction–diffusion system. SIAM J. Appl. Math. 62 (3), 870–887. Hunter, K., Wang, Y., VanCappellen, P., 1997. Kinetic modeling of microbially-driven redox chemistry of subsurface environments: coupling transport, microbial metabolism and geochemistry. J. Hydrol. 209, 53–80. Kitsunezaki, S., 1997. Interface dynamics for bacterial colony formation. J. Phys. Soc. Jpn. 66, 1544–1550. Kuznetsov, Y., Levitin, V., 1996. Content: a multiplatform environment for analyzing dynamical systems. Dynamical Systems Laboratory, Centrum voor Wiskunde en Informatica, Amsterdam. Kuznetsov, S., Mosekilde, E., Dewel, G., Borckmans, P., 1997. Absolute and convective instabilities in a one dimensional Brusselator flow modell. J. Chem. Phys. 106 (18), 7609–7616. Madani, S., Meysman, F., Middelburg, J., 2003. Biogeochemical modeling of sediments from the Santa Barbara Basin (California). In: Rullkötter, J. (Ed.), BioGeoChemistry of Tidal Flats, vol. 12. Forschungszentrum Terramare, Wilhelmshaven, pp. 91–93. Mudryk, Z., Podgórska, B., Ameryk, A., Bolałek, J., 2000. The occurrence and activity of sulfate-reducing bacteria in the bottom sediments of the gulf of Gda´nsk. Oceanologia 42 (1), 105–117. Murray, J., 1993. Mathematical Biology, second ed. Springer, Berlin. Nicolis, G., 1995. Introduction to Nonlinear Science. Cambridge University Press, Cambridge.
Nicolis, G., Prigogine, I., 1977. Self-Organization in Nonequilibrium Systems—From Dissipative Structures to Order through Fluctuations. Wiley, New York. Petrovskii, S., Li, B.-L., Malchow, H., 2003. Quantification of the spatial aspect of chaotic dynamics in biological and chemical systems. Bull. Math. Biol. 65 (3), 425– 446. Rovinsky, A., Menzinger, M., 1994. Differential flow instability in dynamical systems without an unstable (activator) subsystem. Phys. Rev. Lett. 72 (13), 2017–2020. Satnoianu, R., Merkin, J., Scott, S., 1998. Spatio-temporal structures in a differential flow reactor with cubic auto catalator kinetics. Physica D 124, 345–367. Satnoianu, R., Menzinger, M., Maini, P., 2000. Turing instabilities in general systems. Math. Biol. 41 (6), 493–512. Satnoianu, R., Maini, P., Menzinger, M., 2001. Parameter space analysis, pattern sensitivity and model comparison for Turing and stationary flow-distributed waves (FDS). Physica D 160, 79–102. Tsimring, L., Levine, H., Aranson, I., Ben-Jacob, E., Cohen, I., Shochet, O., Reynolds, W., 1995. Aggregation patterns in stressed bacteria. Phys. Rev. Lett. 75 (9), 1859– 1862. Turing, A., 1952. The chemical basis of morphogenesis. Philos. Trans. R. Soc. Lond. B 237 (B 641), 37–72. VanCappellen, P., Wang, Y., 1996. Cycling of iron and manganese in surface sediments: a general theory for the coupled transport and reaction of carbon oxygen, nitrogen, sulfur, iron and manganese. Am. J. Sci. 296, 197–243.