- Email: [email protected]
Spatial enhancement of population uncertainty near the extinction threshold
Ecological Modelling 174 (2004) 191–201
Spatial enhancement of population uncertainty near the extinction threshold Yu Itoh, Kei-ichi Tainaka∗ , Tomoyuki Sakata, Tomomi Tao, Nariyuki Nakagiri Department of Systems Engineering, Shizuoka University, Hamamatsu 432-8561, Japan
Abstract Ecosystem dynamics can be studied using model populations. Perturbation experiments have often been applied to simulated populations in order to study the uncertainty in ecosystem dynamics. Most of these studies have predicted a stationary state. We report on the uncertainty in the dynamics near extinction. In particular, we explore fluctuation enhancements, i.e., enhanced variability in dynamics of paths to extinction. We examine two dynamic models on a two-dimensional lattice of finite size: (1) the contact process (CP) in which interactions are restricted to occur between adjacent lattice points, and (2) mean-field simulation (MFS), where interactions occur globally, between any pair of lattice points. Computer simulation reveals that, for both CP and MFS, the random drift of density about a stationary state increases with the decrease of steady-state density. Drift is much more pronounced in the vicinity of the critical mortality rate, at the transition to extinction. Simulation demonstrates that MFS shows greater variability when relative mortality rate is low whereas CP shows much more pronounced variability when mortality is near the critical threshold. The CP process shows wider fluctuations while the MFS process shows minimal increases following perturbations that lead to extinction. Because interactions are local for CP, there are a variety of different paths to extinction. © 2004 Elsevier B.V. All rights reserved. Keywords: Lattice model; Perturbation experiment; Uncertainty; Small population
1. Introduction Under human management, ecosystems suffer perturbations or disturbances. The investigation of perturbation experiments is essential to conserve species and habitat (Paine, 1966; May, 1973; Tilman and Downing, 1994). The most familiar approach in a perturbation experiment is the press perturbation, where one or more quantities are altered and held at higher or lower levels. We focus on the effects on dynamics of increased mortality rate of a target species. There ∗ Corresponding author. Tel.: +81-53-478-1228; fax: +81-53-478-1228. E-mail address: [email protected] (K.-i. Tainaka).
are a variety of real cases of such perturbations. For example, (i) human development: destruction or fragmentation of habitat often drives species to extinction (World Commission on Dams, 2000); (ii) biological control: people reduce a target species by the use of non-indigenous enemies (Simberloff and Stiling, 1996); (iii) removal experiment: e.g., Paine (1966) periodically removed starfish in experimental plots. In general, it is very difficult to predict the response of ecosystems to applied perturbation. Many authors have reported the indeterminacy of ecological response (Yodzis, 1988; Pimm, 1993; Schoener, 1993; Tainaka, 1994; Schmitz, 1997). Most studies have concentrated on the change of stationary state caused by perturbations. The indeterminacy is thought
0304-3800/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2004.01.004
192
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
to originate from non-linearity or indirect effects (Nakagiri et al., 2001; Abarca-Arenas and Ulanowicz, 2002; Wootton, 2002). The latter are particularly important for our understanding of overall ecosystem functioning (Krivtsov et al., 2000; Krivtsov, 2001, this issue). In the present paper, we do not focus on a stationary state, but examine the transient process to extinction: we describe the uncertainty in the dynamic process. Such uncertainty comes from the property of extinction (phase transition). It is therefore very difficult to predict the extinction of species. We study random drift fluctuation in a small lattice population. In particular, we determine whether or not the so-called fluctuation enhancement (Kubo et al., 1973) occurs. Here fluctuation enhancement, which is observed in the field of physics (Tsuchiya and Horie, 1985), is defined as an increased variance in the process, i.e., wider variability in paths during a phase transition. Since extinction in an ecosystem can be regarded as a phase transition, we may be able to gain a better understanding of extinction by examining this phase transition process in detail. In the present paper, we demonstrate fluctuation enhancement in a lattice population. Local interaction is shown to be essential for the appearance of this phenomenon. To date, the dynamics of small populations of a single species has been studied by various methods. Typical examples are stochastic differential equations and discrete birth–death processes (Renshaw, 1991; Nisbet and Gurney, 1982; Soulé, 1987). When the population size becomes small, the risk of extinction increases. Several authors have estimated the risk of extinction, using the following stochastic differential equation (Lande, 1995): x x˙ = Rx 1 − (1) + Fnoise (x), K where the dot represents the derivative with respect to the time, and the variable x is the population size. The first term in the right hand side of (1) denotes the logistic equation, which includes two parameters R and K, and the last term denotes a noise (random drift) that usually depends on x. If the population size x becomes small, the noise term becomes large. In this approach, however, the intensity of noise is not uniquely determined. On the other hand, the discreteness of the birth–death process can generate noise
(Renshaw, 1991; Soulé, 1987). However, this approach ignores the spatial distribution of individuals. In the present paper, we introduce a spatially explicit model (lattice model). Such a model is more general because it includes local interactions. In recent years, lattice models have been widely applied in the field of ecology (Tainaka, 1988; Matsuda et al., 1992; Nowak et al., 1994; Harada and Iwasa, 1994; Sato et al., 1994; Kubo et al., 1996; Durrett and Levin, 1994; Nakagiri et al., 2001). In the present paper, we do not apply the approach of stochastic equations such as (1), instead we present a simulation model that directly implements random drift. Our model uses a finite lattice system. Two models are implemented. The first, the contact process (CP), a lattice version of the logistic equation, has been extensively studied from mathematical (Harris, 1974; Durrett, 1988; Liggett, 1985) and physical (Katori and Konno, 1991; Marro and Dickman, 1999) aspects. Considerable information on CP has been accumulated, but results are mainly related to a stationary state in an infinitely large domain. However, we focus on dynamic processes in a finite lattice. In the next section, we describe our models and methods in detail. There are two simulation models: the CP, in which interactions occur between neighboring lattice sites, and the mean-field simulation (MFS), where interactions occur globally, between any pair of lattice points. Moreover, we utilize the technique of perturbation experiments. Section 3 reports the results of random drift about a stationary state. In Section 4, results of dynamic processes due to perturbation experiments are described for both CP and MFS. In the final section, we discuss (i) the close relationship between fluctuation enhancement in dynamic processes and random drift about a stationary state, and (ii) the reason why enhancement is stronger in the case of local interactions.
2. Model We consider a single species distributed over a square lattice. Each lattice site is either empty (O) or occupied (X). The site being occupied (X) means that an individual or a sub-population (treated as a single unit) exists at that node. Birth and death processes
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
are given by X + O → 2X, m
→ O. X−
(2a) (2b)
The processes (2a) and (2b) simulate reproduction and death, respectively. The parameter m represents the mortality rate of an individual or the extinction rate of the sub-population. The reactions (2a) and (2b) are carried out in two ways: MFS or CP. We first describe the simulation method for CP:
193
(i) Two lattice sites are randomly and independently chosen; if these sites are X and O or O and X, then the O is changed into X. Now we explain the extinction perturbation experiment. We carry out computer experiments of phase transition. It is known that the CP exhibits a phase transition at the critical point m = mc ; when m increases and exceeds the value mc , the steady-state density changes (continuously) from non-zero value to zero (Marro and Dickman, 1999). Our experiment is performed as follows: the lattice is populated with an initial distribution and an initial mortality rate is chosen that is below the critical value, mc , and the simulation is initiated. Once this initial population has reached a steady state (t = 100), the mortality rate is changed from m = m1 to a new value, m = m2 , that is greater than mc , and is held at this value thereafter. Following this perturbation, the species eventually goes extinct. We repeat this experiment many times, and record the time dependence of species density to determine the level of uncertainty associated with this process.
(1) Initially, we distribute individuals on the square lattice; the initial distribution is not important, since the system evolves into a stationary state. (2) The reactions (2a) and (2b) are performed in the following two steps: (i) we perform two-body reaction (2a): choose one lattice site randomly, and then randomly specify one of four neighboring sites. If the pair of sites are (X, O) or (O, X), then O is changed into X. (ii) we perform one-body reaction (2b). Choose one lattice point randomly; if the site is occupied by X, the site will become O by the rate m. In a real simulation, the maximum mortality max{m} = 2.1. Therefore, we set mortality = m/3 and then perform mortality reaction three times. For example, when m = 2.1, we run mortality reaction three times with m = 2.1/3 = 0.7. (3) Repeat step (2) L × L = 10,000 times, where L × L is the total number of lattice points. This is the Monte Carlo step (Tainaka, 1988). (4) Repeat step (3) until the system reaches a stationary state. In the following work, we use 100 repeats to obtain a stationary state.
We investigate the noise (random drift) about a stationary state. The system approaches a stationary state, where the population size of X fluctuates around the steady-state density (Fig. 1). Fig. 2 shows the typical spatial patterns of stationary states for two values of m; black and white denote the lattice sites of X and O, respectively. We obtain the average As and variance Vs of density x(t) in a stationary state as T2 1 As = x(t) dt, (3) T1 − T2 T1 T2 1 Vs = [x(t) − As (t)]2 dt, (4) T1 − T2 T1
Here we have employed periodic boundary conditions (replacing any index that exceeds the boundary, L, by using the modulo operator on the index; L + 1 is replaced by 1). Next, we describe the method for the MFS in which long-ranged (global) interactions are allowed. MFS is very similar to CP, but step (i) in (2) for CP is replaced with the following:
where T1 and T2 are arbitrary but well separated times chosen after the system has reached a stationary state (T1 << T2 ). In Fig. 3, the average density As and variance Vs are plotted against m; grey and black circles represent MFS and CP, respectively. The steady-state density (K) can be approximated by the average As . Phase transition (between stable state equilibrium and extinction)
3. Random drift in stationary state
194
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
Fig. 1. Typical population dynamics for the contact process (CP). The system, which has 1002 lattice sites, approaches a stationary state where the population size of X fluctuates around the steady-state density (m = 0.7). Similar dynamics are obtained for the mean-field simulation (MFS), although the steady-state densities are different between CP and MFS.
occurs at m = mc . In the case of MFS, we will show that the steady-state density is given by (2 − m)/2 (see the next section); the phase transition point is thus given by mc = 2. On the other hand, for CP, the critical value is approximately mc ∼ 1.21 (Marro and Dickman, 1999). Fig. 3 reveals that the variance Vs increases as the steady-state density, As , decreases. In particular, the variance diverges near the critical point mc (Fig. 4). This phenomenon is associated with the so-called “critical slowing-down” (Marro and Dickman, 1999) which is the divergence of the relaxation time (Binder, 1979) near the critical point, mc . Here the relaxation time is defined in Appendix A. For MFS, we will prove the critical slowing-down (Appendix A). Note that the variance or the intensity of noise is uniquely determined for both MFS and CP, so long as both values of death rate m and the total number of lattice points are given. From Fig. 4, the relation between Vs and As in low densities can be approximated by Vs ∝ A−α s ,
(5)
where α is a positive constant. We have α ∼ 0.5 for MFS and α ∼ 0.8 for CP. Thus, the variance (also referred to as random drift) for CP grows more rapidly than that for MFS, near extinction.
Fig. 2. Snapshots of typical stationary patterns for CP: (a) m = 0.3, and (b) m = 1.0. If m < mc , the species X survives. Black and white denote the lattice sites of X and O, respectively. These mortality rates are below the critical point mc , so that the species X will survive (mc ≈ 1.21). In the case of (a), the spatial distribution is nearly random. In the case of (b), individuals are spatially clumped.
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
195
Fig. 3. The average density As and variance Vs from simulation results at stationary states (L × L = 1002 ), for various values of the relative death rate, m/mc . Grey and black denote MFS and CP, respectively. Both As and Vs are averaged over 200 < t < 500. mc = 2.0 for MFS and while for CP, the exact value of mc is unknown, we assume mc = 1.21.
4. Results of perturbation experiments 4.1. Result of MFS Perturbation experiments are repeatedly performed from m1 to m2 , where m2 > mc and m1 < mc . We prepare N initial spatial patterns (ensembles); each of
them has a distinct initial density, xi (0) (i = 1, 2, . . . , N), but these initial conditions do not effect the results. Each simulation proceeds with m = m1 for 100 iterations; the system evolves into a stationary state. The value of xi (100) is very close to the steady-state density at m = m1 . After the perturbation (t > 100), m = m2 . We calculate the density xi (t) in order to
196
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
sites [(X and O) or (O and X)]. The second term of (8) denotes the death process. Eq. (8) leads to the following logistic equation: x x˙ = Rx 1 − , (9) K where the quantities R and K are defined by R = 2 − m,
Fig. 4. Log–log relationship between the variance, Vs , and the steady-state density, As . These data are the same as those presented in Fig. 3; Grey represent MFS and black represent CP values.
obtain the time dependencies of the average A(t) and the variance V(t) which are defined by 1 A(t) = xi (t), (6) N i
V(t) =
1 [xi (t) − A(t)]2 . N
(7)
i
We want to emphasize that the quantities A(t) and V(t) are averaged over different ensembles, while both As and Vs denote the time averages defined by (3) and (4). As demonstrated in the previous section, the random drift about the stationary state increases with decrease of the steady-state density; especially near the extinction point mc , the random drift becomes extremely high. It is, therefore, expected that indeterminacy occurs in the transient state of the dynamic process. To confirm the indeterminacy, we explore the fluctuation enhancement, which constitutes an extreme increase in the value of V(t). If the size of system is infinitely large (L → ∞), then the population dynamics for MFS is given by x˙ = 2x(1 − x) − mx,
(8)
where the dot represents the time derivative which is measured by the Monte Carlo step (Tainaka, 1988) and x is the density of X. The first term of right hand side of (8) represents the birth process; the factor (1−x) in this term is the density of empty site, and the coefficient 2 means that there are two ways to select a pair of
K = 21 (2 − m).
The steady-state solution of the logistic Eq. (9) can be obtained by setting the time derivative in (9) to zero: the non-trivial density in the stationary state is given by x = K = (2 − m)/2. This is consistent with the result of MFS (Fig. 3): the average density As obtained by simulation is almost equal to the steady-state density K derived by mean-field theory. We carry out perturbation experiments of extinction in MFS. It should be emphasized that the system size L is finite (L = 100). In this case, the dynamics becomes stochastic. In Fig. 5, a typical result of a perturbation experiment is displayed; the time dependencies of both average density A(t) and the variance V(t), which are calculated by simulation, are plotted. It is found from computer simulation that the variance V(t) never increases in the dynamic process of perturbation experiment. Clear fluctuation enhancement never takes place for extreme values of m1 and m2 ; especially when m2 is much larger than mc ; both A(t) and V(t) rapidly decrease with time. 4.2. Result of CP The basic equation of CP for L → ∞ becomes; x˙ = (PXO − POX ) − mx,
(10)
where Pij is the probability density of finding a state i at a site and a state j at a nearest neighbor of the former site (i, j = X, O). Note that the following relations hold: Pij = Pji , Pij = Pi , (11) j
where PX = x, PO = (1 − x). If we assume the mean-field theory [Pij = Pi Pj or PXO = x(1 − x)], Eq. (10) becomes (8). Unfortunately, the basic Eq. (10) for CP cannot be solved. Nevertheless, it is known that the similar extinction occurs as predicted by the mean-field theory (8).
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
197
Fig. 5. The results of 100 repeated perturbation experiments in the mean-field simulation (MFS). Different initial patterns were applied to a fixed lattice size, L × L = 1002 . The time dependencies of both average A(t) and variance V(t) (defined by Eqs. (6) and (7)) are displayed. The solid line represents A(t) and the points represent V(t). V(t) values are multiplied by 106 . At time t = 100, the death rate m changes from m1 = 1.0 to m2 = 2.1 (mc = 2.0). It can be seen that the variance V(t) never increases.
We carry out perturbation experiments for CP: when t = 100, the mortality rate m is raised from m1 to m2 (m1 < mc < m2 ). In Fig. 6, the time dependencies of both average density A(t) and the variance V(t) are plotted, where m1 = 0.5 and (a) m2 = 1.3, (b) m2 = 2.0 and (c) m2 = 5.0. Fig. 6a clearly demonstrates that the fluctuation enhancement takes place. Even though the average density decreases, the variance V(t) strongly increases during the intermediate stage of phase transition. Note that such a phenomenon is not always observed. In Fig. 6b and c, no fluctuation enhancement is observed. A fluctuation enhancement is not always caused by a critical slowing-down; it is not observed in the mean-field limit (see Fig. 5). We associate enhancement with the formation of aggregates or clusters of individuals (Fig. 2b): this clumping behavior causes variation to increase. To explore the degree of clumping, we obtain the ratio RXX =
PXX , x2
(12)
near the stationary state (Tainaka, 1994). When the distribution of individuals is random, RXX = 1. When RXX > 1 the distribution is clustered, whilst when RXX < 1, the distribution is uniform. The basic
Eq. (10) in stationary state becomes PXO = 21 mK,
(13)
because we put x˙ = 0. On the other hand, the second equation of (11) implies that PXX = K − PXO .
(14)
From (12)–(14), we have RXX ∝ K−1 . Near the critical point (m → mc ), the steady-state density K satisfies K ∝ (mc − m)β , where β is a positive constant (Marro and Dickman, 1999). It follows that RXX ∝ (mc − m)−β .
(15)
This relation represents the clumping behavior. When m approaches mc , the degree of clumping of X increases rapidly (RXX → ∞). Such clumping behavior (15) may cause the abrupt increase of random drift near extinction (Fig. 4). Similarly, the enhancement of V(t) can be caused by clumping behavior. These results are suggesting that when m2 is near the critical state mc , then the transition to extinction is slow and the large variance is exposed during the slow transition.
198
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
5. Discussion Simulated perturbation experiments are useful for the studying of ecosystem response. Study of the process to extinction is very important for conservation biology. In the present article, we study perturbation experiments of phase transition (extinction) to report the uncertainty (fluctuation enhancement). In these experiments, the mortality rate m of individuals is suddenly increased from m1 to m2 , where m1 < mc < m2 . The density x of individuals is thus changed from a positive value to zero. We estimate the fluctuation in population size (density). Simulations are carried out by two different methods: the CP and the MFS. In the former case, interaction is restricted to occur between adjacent lattice points, whereas in the latter case, interaction globally occurs between any pair of lattice points. For a CP, the fluctuation enhancement is clearly observed as illustrated in Fig. 6a. The fluctuation enhancement means that there are a variety of routes to extinction. It is not easy to predict when extinction will occur because of this uncertainty associated with this process. It is emphasized for the CP that the clearest enhancement appears when m2 takes a value near the critical point mc . If the mortality rate m is increased to a point far above mc , then no enhancement occurs. In these cases, it is relatively easy to predict the time to extinction (see Fig. 6b and c). However, such an extreme increase of mortality rate is very rare in real ecosystems. Some intermediate level of removal of a species is more common, so that the prediction of extinction becomes very hard. Fluctuation enhancement is associated with not only the random drift but also the clumping behavior (Eqs. (12) and (15)) in stationary state. It is found that random drift increases with decrease in steady-state density (Figs. 3 and 4). Such increase may be related to critical slowing-down (see Appendix A). Moreover, we find in the vicinity of extinction that the increase of random drift for a CP is more pronounced than that for a MFS. This is caused by the spatial distribution of individuals: when m approaches the
199
extinction point mc , the degree of clumping of species X in stationary state is extreme (see Eq. (15)). When m2 ∼ mc , the CP system in dynamical process more or less has the property of critical behavior (extreme increase of variance Vs near the extinction threshold) in stationary state. For this reason, the enhancement may be clearly observed for CP. In the Appendix we attempt to explain why fluctuation enhancement is not clear in the MFS system. It should be emphasized that MFS ignores the spatial distribution of individuals. To date, fluctuation enhancement has been observed in some physical systems (Kubo et al., 1973; Tsuchiya and Horie, 1985). In such systems, values of order parameters (in this case “density”) are known to increase from zero to a positive value. In contrast, we have shown just the opposite case: even though density is decreased, the fluctuation enhancement occurs. In general, if the density would increase, the variance V(t) often increased. An example is the exponential growth; the density difference at initial condition always increases. Our results suggest that both critical slowing-down and the clumping behavior of individuals may be essential to fluctuation enhancement. Almost all species form a spatially clumped distribution. Especially, when a species becomes endangered, the degree of clumping usually increases. This is due to the inherent nature of biological species: an offspring is produced in the neighborhood of parents. Here we discuss fluctuation enhancement. The definition should be made clearer, because increases of variance, V(t), have been observed for other reasons. For example, it has been shown that the variance is large for the stochastic logistic equation when R = 0 and K takes a large value (Renshaw, 1991). This large variance originates as a result of Brownian drift. In contrast, the fluctuation enhancement presented here is a new source of variability. We have shown that fluctuation enhancement is due to the spatial characteristics of the process; only local interactions give rise to increased fluctuations. The clumping behavior of individuals (Eq. (15)) near the extinction point appears to play an important role.
Fig. 6. Results of 100 perturbation experiments for the contact process (CP). The solid line represents A(t) and the points represent V(t). The critical value of the mortality parameter (point of phase transition) is mc ∼ 1.21. The death rate m is changed from m1 = 0.5 to m2 . (a) m2 = 1.3, (b) m2 = 2.0, (c) m2 = 5.0, (a) clearly exhibits this fluctuation enhancement (increase of V(t)). When m2 takes large values, no fluctuation enhancement occurs as can be seen in (b) and (c).
200
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
Finally, we propose the following experiment. Consider the cultivation of bacteria. Culture bacteria in two nutrient media: agar and a fluid medium. The growth of bacteria in the agar medium parallels the CP in which the growth process (2a) occurs locally: offspring are reproduced at a sites adjacent to their parents. In the fluid medium, the solution is continuously stirred, so that offspring are randomly dispersed in the medium. Thus, interactions in this medium correspond to a MFS. The reaction (2b) corresponds to the death processes of bacteria in the actual experiment; (2b) is realized using UV light or radioactive rays. The parameter m represents the intensity of light (radiation). Our results predict that there are a variety of different processes to extinction in the case of agar medium (CP) and that the stirred liquid medium will show a smooth transition to extinction.
Acknowledgements
x(t) − K ∝ exp(−(mc − m)t).
Appendix A. Mean-field theory We prove the critical slowing-down, and explain why fluctuation enhancement almost disappears in mean-field theory. Assume that the total number of sites is infinitely large (L → ∞); the dynamic equation for MFS becomes the logistic Eq. (9). The steady-state solution can be obtained by setting the time derivative in (9) to zero. The non-trivial solution (x = K) for 0 < m < 2 is stable, and phase transition occurs at mc = 2. Consider the density x(t) at time t, and assume that the system initially stays near the stationary state (x(0) ∼ K). Expand x around the steady-state density (K), x(t) = K + X(t), where X takes small values (X 1). Inserting this equation into (9) and collecting first order terms to obtain (A.1)
(A.2)
From this equation, the relaxation time is defined by 1/(mc − m). Hence, we prove that the relaxation time diverges as m → mc . This is called the critical slowing-down: the variance (V) in stationary state may diverge near m = mc . Next, we explain why fluctuation enhancement almost disappears in the mean-field limit. It is known in the field of physics that there is a close relation between critical slowing-down and fluctuation enhancement (Kubo et al., 1973; Tsuchiya and Horie, 1985). We carry out the same perturbation experiment N times, with initial densities given by xi (0) (i = 1, 2, . . . , N). After expanding and rearranging terms in the definition of V(t), we have 2 1 1 V(t) = xi (t)2 − xi (t) N N i
The authors are grateful to Professors Jin Yoshimura for valuable discussions. They also thank to referees for many useful comments and criticisms.
˙ = −(mc − m)X. X
From this equation it follows that
1 = [xi (t) − xj (t)]2 . N
i
i>j
Taking the derivative with respect to time, t, it follows that 2 V˙ = 2 (xi − xj )(˙xi − x˙ j ). (A.3) N i>j
The dynamics for t > 0 is assumed to be determined by the logistic Eq. (9), and following the perturbation, both parameters R and K are negative (m2 > mc ). Let xi (t) and xj (t) be two solutions of (9) whose initial values are slightly different from each other. Without loss of generality, we assume that xi (0) > xj (0). From (9), we have xi + x j x˙ i − xj = R(xi − xj ) 1 − . (A.4) K Hence, the right-hand side of (A.4) is always negative, and the difference xi (t) − xj (t) decreases as the time proceeds. It is therefore proved from (A.4) that the variance V(t) always decreases. Strictly speaking, this proof is insufficient, since the dynamics for a finite lattice (finite value of L) is described by a stochastic differential Eq. (1). Nevertheless, the logistic Eq. (9)
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
explains well how (i) critical slowing-down occurs, and why (ii) enhancement of fluctuation almost disappears in the mean-field limit. References Abarca-Arenas, L.G., Ulanowicz, R.E., 2002. The effects of taxonomic aggregation on network analysis. Ecol. Model. 149, 285–296. Binder, K. (Ed.), 1979. Applications of Monte Carlo Methods in Statistical Physics. Springer-Verlag, Berlin. Durrett, R., 1988. Lecture Notes on Particle Systems and Percolation. Wadsworth and Brooka/Cole Advanced Books & Software, California. Durrett, R., Levin, S.A., 1994. The importance of being discrete (and spatial). Theor. Popul. Biol. 46, 363–394. Harada, Y., Iwasa, Y., 1994. Lattice population dynamics for plants with dispersing seeds and vegetative propagation. Res. Popul. Ecol. 36, 237–249. Harris, T.E., 1974. Contact interaction on a lattice. Ann. Prob. 2, 969–988. Katori, M., Konno, N., 1991. Upper bounds for survival probability of the contact process. J. Stat. Phys. 63, 115–130. Krivtsov, V., 2001. Study of cause-and-effect relationships in the formation of biocenoses: their use for the control of eutrophication. Russian J. Ecol. 32, 230–234. Krivtsov, V., Corliss, J., Bellinger, E., Sigee, D., 2000. Indirect regulation rule for consecutive stages of ecological succession. Ecol. Model. 133, 73–82. Kubo, R., Matsuo, K., Kitahara, K., 1973. Fluctuation and relaxation of macrovariables. J. Stat. Phys. 9, 51–96. Kubo, T., Iwasa, Y., Furumoto, N., 1996. Forest spatial dynamics with gap expansions: total gap area and gap size distribution. J. Theor. Biol. 180, 229–246. Lande, R., 1995. Mutation and conservation. Conserv. Biol. 9, 782–791. Liggett, T.M., 1985. Interacting Particle Systems. Springer-Verlag, New York. Marro, J., Dickman, R., 1999. Nonequilibrium Phase Transition in Lattice Models. Cambridge University Press, Cambridge. Matsuda, H., Ogita, N., Sasaki, A., Sato, K., 1992. Statistical mechanics of population: the lattice Lotka–Volterra model. Prog. Theor. Phys. 88, 1035–1049.
201
May, R.M., 1973. Stability and Complexity in Model Ecosystems. Princeton University Press, Princeton. Nakagiri, N., Tainaka, K., Tao, T., 2001. Indirect relation between species extinction and habitat destruction. Ecol. Model. 137, 109–118. Nisbet, R.M., Gurney, W.S.C., 1982. Modelling Fluctuating Populations. Wiley, New York. Nowak, M.A., Bonhoeffer, S., May, R.M., 1994. More spatial games. Int. J. Bifurcation Chaos 4, 33–56. Paine, R.T., 1966. Food web complexity and species diversity. Am. Nat. 100, 65–75. Pimm, S.L., 1993. Discussion: understanding indirect effects: is it possible? In: Kawanabe, H., Cohen, J.E., Iwasaki, K. (Eds.), Mutualism and Community Organization. Oxford University Press, Oxford, pp. 199–209. Renshaw, E., 1991. Modelling Biological Populations in Space and Time. Cambridge University Press, Cambridge. Sato, K., Matsuda, H., Sasaki, A., 1994. Pathogen invasion and host extinction in lattice structured populations. J. Math. Biol. 32, 251–268. Schmitz, O.J., 1997. Press perturbations and the predictability of ecological interactions in a food web. Ecology 78, 55–69. Schoener, T.W., 1993. On the relative importance of direct vs. indirect effects in ecological communities. In: Kawanabe, H., Cohen, J.E., Iwasaki, K. (Eds.), Mutualism and Community Organization. Oxford University Press, Oxford, pp. 365–411. Simberloff, D., Stiling, P., 1996. How risky is biological control? Ecology 77, 1965–1974. Soulé, M.E., 1987. Viable Populations for Conservation. Cambridge University Press, Cambridge. Tainaka, K., 1988. Lattice model for the Lotka–Volterra system. J. Phys. Soc. Jpn. 57, 2588–2590. Tainaka, K., 1994. Intrinsic uncertainty in ecological catastrophe. J. Theor. Biol. 166, 91–99. Tilman, D., Downing, J.A., 1994. Biodiversity and stability in grasslands. Nature 367, 363–365. Tsuchiya, Y., Horie, S., 1985. Evolution process of the Williams domain in a nematic liquid crystal. J. Phys. Soc. Jpn. 54, 1–4. Wootton, J.T., 2002. Indirect effects in complex ecosystems: recent progress and future challenges. J. Sea Res. 48, 157–172. World Commission on Dams, 2000. Dams and Development: A New Framework for Decision-Making. Earthscan, London. Yodzis, P., 1988. The indeterminacy of ecological interactions as perceived through perturbation experiments. Ecology 69, 508– 515.
Spatial enhancement of population uncertainty near the extinction threshold Yu Itoh, Kei-ichi Tainaka∗ , Tomoyuki Sakata, Tomomi Tao, Nariyuki Nakagiri Department of Systems Engineering, Shizuoka University, Hamamatsu 432-8561, Japan
Abstract Ecosystem dynamics can be studied using model populations. Perturbation experiments have often been applied to simulated populations in order to study the uncertainty in ecosystem dynamics. Most of these studies have predicted a stationary state. We report on the uncertainty in the dynamics near extinction. In particular, we explore fluctuation enhancements, i.e., enhanced variability in dynamics of paths to extinction. We examine two dynamic models on a two-dimensional lattice of finite size: (1) the contact process (CP) in which interactions are restricted to occur between adjacent lattice points, and (2) mean-field simulation (MFS), where interactions occur globally, between any pair of lattice points. Computer simulation reveals that, for both CP and MFS, the random drift of density about a stationary state increases with the decrease of steady-state density. Drift is much more pronounced in the vicinity of the critical mortality rate, at the transition to extinction. Simulation demonstrates that MFS shows greater variability when relative mortality rate is low whereas CP shows much more pronounced variability when mortality is near the critical threshold. The CP process shows wider fluctuations while the MFS process shows minimal increases following perturbations that lead to extinction. Because interactions are local for CP, there are a variety of different paths to extinction. © 2004 Elsevier B.V. All rights reserved. Keywords: Lattice model; Perturbation experiment; Uncertainty; Small population
1. Introduction Under human management, ecosystems suffer perturbations or disturbances. The investigation of perturbation experiments is essential to conserve species and habitat (Paine, 1966; May, 1973; Tilman and Downing, 1994). The most familiar approach in a perturbation experiment is the press perturbation, where one or more quantities are altered and held at higher or lower levels. We focus on the effects on dynamics of increased mortality rate of a target species. There ∗ Corresponding author. Tel.: +81-53-478-1228; fax: +81-53-478-1228. E-mail address: [email protected] (K.-i. Tainaka).
are a variety of real cases of such perturbations. For example, (i) human development: destruction or fragmentation of habitat often drives species to extinction (World Commission on Dams, 2000); (ii) biological control: people reduce a target species by the use of non-indigenous enemies (Simberloff and Stiling, 1996); (iii) removal experiment: e.g., Paine (1966) periodically removed starfish in experimental plots. In general, it is very difficult to predict the response of ecosystems to applied perturbation. Many authors have reported the indeterminacy of ecological response (Yodzis, 1988; Pimm, 1993; Schoener, 1993; Tainaka, 1994; Schmitz, 1997). Most studies have concentrated on the change of stationary state caused by perturbations. The indeterminacy is thought
0304-3800/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2004.01.004
192
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
to originate from non-linearity or indirect effects (Nakagiri et al., 2001; Abarca-Arenas and Ulanowicz, 2002; Wootton, 2002). The latter are particularly important for our understanding of overall ecosystem functioning (Krivtsov et al., 2000; Krivtsov, 2001, this issue). In the present paper, we do not focus on a stationary state, but examine the transient process to extinction: we describe the uncertainty in the dynamic process. Such uncertainty comes from the property of extinction (phase transition). It is therefore very difficult to predict the extinction of species. We study random drift fluctuation in a small lattice population. In particular, we determine whether or not the so-called fluctuation enhancement (Kubo et al., 1973) occurs. Here fluctuation enhancement, which is observed in the field of physics (Tsuchiya and Horie, 1985), is defined as an increased variance in the process, i.e., wider variability in paths during a phase transition. Since extinction in an ecosystem can be regarded as a phase transition, we may be able to gain a better understanding of extinction by examining this phase transition process in detail. In the present paper, we demonstrate fluctuation enhancement in a lattice population. Local interaction is shown to be essential for the appearance of this phenomenon. To date, the dynamics of small populations of a single species has been studied by various methods. Typical examples are stochastic differential equations and discrete birth–death processes (Renshaw, 1991; Nisbet and Gurney, 1982; Soulé, 1987). When the population size becomes small, the risk of extinction increases. Several authors have estimated the risk of extinction, using the following stochastic differential equation (Lande, 1995): x x˙ = Rx 1 − (1) + Fnoise (x), K where the dot represents the derivative with respect to the time, and the variable x is the population size. The first term in the right hand side of (1) denotes the logistic equation, which includes two parameters R and K, and the last term denotes a noise (random drift) that usually depends on x. If the population size x becomes small, the noise term becomes large. In this approach, however, the intensity of noise is not uniquely determined. On the other hand, the discreteness of the birth–death process can generate noise
(Renshaw, 1991; Soulé, 1987). However, this approach ignores the spatial distribution of individuals. In the present paper, we introduce a spatially explicit model (lattice model). Such a model is more general because it includes local interactions. In recent years, lattice models have been widely applied in the field of ecology (Tainaka, 1988; Matsuda et al., 1992; Nowak et al., 1994; Harada and Iwasa, 1994; Sato et al., 1994; Kubo et al., 1996; Durrett and Levin, 1994; Nakagiri et al., 2001). In the present paper, we do not apply the approach of stochastic equations such as (1), instead we present a simulation model that directly implements random drift. Our model uses a finite lattice system. Two models are implemented. The first, the contact process (CP), a lattice version of the logistic equation, has been extensively studied from mathematical (Harris, 1974; Durrett, 1988; Liggett, 1985) and physical (Katori and Konno, 1991; Marro and Dickman, 1999) aspects. Considerable information on CP has been accumulated, but results are mainly related to a stationary state in an infinitely large domain. However, we focus on dynamic processes in a finite lattice. In the next section, we describe our models and methods in detail. There are two simulation models: the CP, in which interactions occur between neighboring lattice sites, and the mean-field simulation (MFS), where interactions occur globally, between any pair of lattice points. Moreover, we utilize the technique of perturbation experiments. Section 3 reports the results of random drift about a stationary state. In Section 4, results of dynamic processes due to perturbation experiments are described for both CP and MFS. In the final section, we discuss (i) the close relationship between fluctuation enhancement in dynamic processes and random drift about a stationary state, and (ii) the reason why enhancement is stronger in the case of local interactions.
2. Model We consider a single species distributed over a square lattice. Each lattice site is either empty (O) or occupied (X). The site being occupied (X) means that an individual or a sub-population (treated as a single unit) exists at that node. Birth and death processes
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
are given by X + O → 2X, m
→ O. X−
(2a) (2b)
The processes (2a) and (2b) simulate reproduction and death, respectively. The parameter m represents the mortality rate of an individual or the extinction rate of the sub-population. The reactions (2a) and (2b) are carried out in two ways: MFS or CP. We first describe the simulation method for CP:
193
(i) Two lattice sites are randomly and independently chosen; if these sites are X and O or O and X, then the O is changed into X. Now we explain the extinction perturbation experiment. We carry out computer experiments of phase transition. It is known that the CP exhibits a phase transition at the critical point m = mc ; when m increases and exceeds the value mc , the steady-state density changes (continuously) from non-zero value to zero (Marro and Dickman, 1999). Our experiment is performed as follows: the lattice is populated with an initial distribution and an initial mortality rate is chosen that is below the critical value, mc , and the simulation is initiated. Once this initial population has reached a steady state (t = 100), the mortality rate is changed from m = m1 to a new value, m = m2 , that is greater than mc , and is held at this value thereafter. Following this perturbation, the species eventually goes extinct. We repeat this experiment many times, and record the time dependence of species density to determine the level of uncertainty associated with this process.
(1) Initially, we distribute individuals on the square lattice; the initial distribution is not important, since the system evolves into a stationary state. (2) The reactions (2a) and (2b) are performed in the following two steps: (i) we perform two-body reaction (2a): choose one lattice site randomly, and then randomly specify one of four neighboring sites. If the pair of sites are (X, O) or (O, X), then O is changed into X. (ii) we perform one-body reaction (2b). Choose one lattice point randomly; if the site is occupied by X, the site will become O by the rate m. In a real simulation, the maximum mortality max{m} = 2.1. Therefore, we set mortality = m/3 and then perform mortality reaction three times. For example, when m = 2.1, we run mortality reaction three times with m = 2.1/3 = 0.7. (3) Repeat step (2) L × L = 10,000 times, where L × L is the total number of lattice points. This is the Monte Carlo step (Tainaka, 1988). (4) Repeat step (3) until the system reaches a stationary state. In the following work, we use 100 repeats to obtain a stationary state.
We investigate the noise (random drift) about a stationary state. The system approaches a stationary state, where the population size of X fluctuates around the steady-state density (Fig. 1). Fig. 2 shows the typical spatial patterns of stationary states for two values of m; black and white denote the lattice sites of X and O, respectively. We obtain the average As and variance Vs of density x(t) in a stationary state as T2 1 As = x(t) dt, (3) T1 − T2 T1 T2 1 Vs = [x(t) − As (t)]2 dt, (4) T1 − T2 T1
Here we have employed periodic boundary conditions (replacing any index that exceeds the boundary, L, by using the modulo operator on the index; L + 1 is replaced by 1). Next, we describe the method for the MFS in which long-ranged (global) interactions are allowed. MFS is very similar to CP, but step (i) in (2) for CP is replaced with the following:
where T1 and T2 are arbitrary but well separated times chosen after the system has reached a stationary state (T1 << T2 ). In Fig. 3, the average density As and variance Vs are plotted against m; grey and black circles represent MFS and CP, respectively. The steady-state density (K) can be approximated by the average As . Phase transition (between stable state equilibrium and extinction)
3. Random drift in stationary state
194
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
Fig. 1. Typical population dynamics for the contact process (CP). The system, which has 1002 lattice sites, approaches a stationary state where the population size of X fluctuates around the steady-state density (m = 0.7). Similar dynamics are obtained for the mean-field simulation (MFS), although the steady-state densities are different between CP and MFS.
occurs at m = mc . In the case of MFS, we will show that the steady-state density is given by (2 − m)/2 (see the next section); the phase transition point is thus given by mc = 2. On the other hand, for CP, the critical value is approximately mc ∼ 1.21 (Marro and Dickman, 1999). Fig. 3 reveals that the variance Vs increases as the steady-state density, As , decreases. In particular, the variance diverges near the critical point mc (Fig. 4). This phenomenon is associated with the so-called “critical slowing-down” (Marro and Dickman, 1999) which is the divergence of the relaxation time (Binder, 1979) near the critical point, mc . Here the relaxation time is defined in Appendix A. For MFS, we will prove the critical slowing-down (Appendix A). Note that the variance or the intensity of noise is uniquely determined for both MFS and CP, so long as both values of death rate m and the total number of lattice points are given. From Fig. 4, the relation between Vs and As in low densities can be approximated by Vs ∝ A−α s ,
(5)
where α is a positive constant. We have α ∼ 0.5 for MFS and α ∼ 0.8 for CP. Thus, the variance (also referred to as random drift) for CP grows more rapidly than that for MFS, near extinction.
Fig. 2. Snapshots of typical stationary patterns for CP: (a) m = 0.3, and (b) m = 1.0. If m < mc , the species X survives. Black and white denote the lattice sites of X and O, respectively. These mortality rates are below the critical point mc , so that the species X will survive (mc ≈ 1.21). In the case of (a), the spatial distribution is nearly random. In the case of (b), individuals are spatially clumped.
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
195
Fig. 3. The average density As and variance Vs from simulation results at stationary states (L × L = 1002 ), for various values of the relative death rate, m/mc . Grey and black denote MFS and CP, respectively. Both As and Vs are averaged over 200 < t < 500. mc = 2.0 for MFS and while for CP, the exact value of mc is unknown, we assume mc = 1.21.
4. Results of perturbation experiments 4.1. Result of MFS Perturbation experiments are repeatedly performed from m1 to m2 , where m2 > mc and m1 < mc . We prepare N initial spatial patterns (ensembles); each of
them has a distinct initial density, xi (0) (i = 1, 2, . . . , N), but these initial conditions do not effect the results. Each simulation proceeds with m = m1 for 100 iterations; the system evolves into a stationary state. The value of xi (100) is very close to the steady-state density at m = m1 . After the perturbation (t > 100), m = m2 . We calculate the density xi (t) in order to
196
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
sites [(X and O) or (O and X)]. The second term of (8) denotes the death process. Eq. (8) leads to the following logistic equation: x x˙ = Rx 1 − , (9) K where the quantities R and K are defined by R = 2 − m,
Fig. 4. Log–log relationship between the variance, Vs , and the steady-state density, As . These data are the same as those presented in Fig. 3; Grey represent MFS and black represent CP values.
obtain the time dependencies of the average A(t) and the variance V(t) which are defined by 1 A(t) = xi (t), (6) N i
V(t) =
1 [xi (t) − A(t)]2 . N
(7)
i
We want to emphasize that the quantities A(t) and V(t) are averaged over different ensembles, while both As and Vs denote the time averages defined by (3) and (4). As demonstrated in the previous section, the random drift about the stationary state increases with decrease of the steady-state density; especially near the extinction point mc , the random drift becomes extremely high. It is, therefore, expected that indeterminacy occurs in the transient state of the dynamic process. To confirm the indeterminacy, we explore the fluctuation enhancement, which constitutes an extreme increase in the value of V(t). If the size of system is infinitely large (L → ∞), then the population dynamics for MFS is given by x˙ = 2x(1 − x) − mx,
(8)
where the dot represents the time derivative which is measured by the Monte Carlo step (Tainaka, 1988) and x is the density of X. The first term of right hand side of (8) represents the birth process; the factor (1−x) in this term is the density of empty site, and the coefficient 2 means that there are two ways to select a pair of
K = 21 (2 − m).
The steady-state solution of the logistic Eq. (9) can be obtained by setting the time derivative in (9) to zero: the non-trivial density in the stationary state is given by x = K = (2 − m)/2. This is consistent with the result of MFS (Fig. 3): the average density As obtained by simulation is almost equal to the steady-state density K derived by mean-field theory. We carry out perturbation experiments of extinction in MFS. It should be emphasized that the system size L is finite (L = 100). In this case, the dynamics becomes stochastic. In Fig. 5, a typical result of a perturbation experiment is displayed; the time dependencies of both average density A(t) and the variance V(t), which are calculated by simulation, are plotted. It is found from computer simulation that the variance V(t) never increases in the dynamic process of perturbation experiment. Clear fluctuation enhancement never takes place for extreme values of m1 and m2 ; especially when m2 is much larger than mc ; both A(t) and V(t) rapidly decrease with time. 4.2. Result of CP The basic equation of CP for L → ∞ becomes; x˙ = (PXO − POX ) − mx,
(10)
where Pij is the probability density of finding a state i at a site and a state j at a nearest neighbor of the former site (i, j = X, O). Note that the following relations hold: Pij = Pji , Pij = Pi , (11) j
where PX = x, PO = (1 − x). If we assume the mean-field theory [Pij = Pi Pj or PXO = x(1 − x)], Eq. (10) becomes (8). Unfortunately, the basic Eq. (10) for CP cannot be solved. Nevertheless, it is known that the similar extinction occurs as predicted by the mean-field theory (8).
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
197
Fig. 5. The results of 100 repeated perturbation experiments in the mean-field simulation (MFS). Different initial patterns were applied to a fixed lattice size, L × L = 1002 . The time dependencies of both average A(t) and variance V(t) (defined by Eqs. (6) and (7)) are displayed. The solid line represents A(t) and the points represent V(t). V(t) values are multiplied by 106 . At time t = 100, the death rate m changes from m1 = 1.0 to m2 = 2.1 (mc = 2.0). It can be seen that the variance V(t) never increases.
We carry out perturbation experiments for CP: when t = 100, the mortality rate m is raised from m1 to m2 (m1 < mc < m2 ). In Fig. 6, the time dependencies of both average density A(t) and the variance V(t) are plotted, where m1 = 0.5 and (a) m2 = 1.3, (b) m2 = 2.0 and (c) m2 = 5.0. Fig. 6a clearly demonstrates that the fluctuation enhancement takes place. Even though the average density decreases, the variance V(t) strongly increases during the intermediate stage of phase transition. Note that such a phenomenon is not always observed. In Fig. 6b and c, no fluctuation enhancement is observed. A fluctuation enhancement is not always caused by a critical slowing-down; it is not observed in the mean-field limit (see Fig. 5). We associate enhancement with the formation of aggregates or clusters of individuals (Fig. 2b): this clumping behavior causes variation to increase. To explore the degree of clumping, we obtain the ratio RXX =
PXX , x2
(12)
near the stationary state (Tainaka, 1994). When the distribution of individuals is random, RXX = 1. When RXX > 1 the distribution is clustered, whilst when RXX < 1, the distribution is uniform. The basic
Eq. (10) in stationary state becomes PXO = 21 mK,
(13)
because we put x˙ = 0. On the other hand, the second equation of (11) implies that PXX = K − PXO .
(14)
From (12)–(14), we have RXX ∝ K−1 . Near the critical point (m → mc ), the steady-state density K satisfies K ∝ (mc − m)β , where β is a positive constant (Marro and Dickman, 1999). It follows that RXX ∝ (mc − m)−β .
(15)
This relation represents the clumping behavior. When m approaches mc , the degree of clumping of X increases rapidly (RXX → ∞). Such clumping behavior (15) may cause the abrupt increase of random drift near extinction (Fig. 4). Similarly, the enhancement of V(t) can be caused by clumping behavior. These results are suggesting that when m2 is near the critical state mc , then the transition to extinction is slow and the large variance is exposed during the slow transition.
198
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
5. Discussion Simulated perturbation experiments are useful for the studying of ecosystem response. Study of the process to extinction is very important for conservation biology. In the present article, we study perturbation experiments of phase transition (extinction) to report the uncertainty (fluctuation enhancement). In these experiments, the mortality rate m of individuals is suddenly increased from m1 to m2 , where m1 < mc < m2 . The density x of individuals is thus changed from a positive value to zero. We estimate the fluctuation in population size (density). Simulations are carried out by two different methods: the CP and the MFS. In the former case, interaction is restricted to occur between adjacent lattice points, whereas in the latter case, interaction globally occurs between any pair of lattice points. For a CP, the fluctuation enhancement is clearly observed as illustrated in Fig. 6a. The fluctuation enhancement means that there are a variety of routes to extinction. It is not easy to predict when extinction will occur because of this uncertainty associated with this process. It is emphasized for the CP that the clearest enhancement appears when m2 takes a value near the critical point mc . If the mortality rate m is increased to a point far above mc , then no enhancement occurs. In these cases, it is relatively easy to predict the time to extinction (see Fig. 6b and c). However, such an extreme increase of mortality rate is very rare in real ecosystems. Some intermediate level of removal of a species is more common, so that the prediction of extinction becomes very hard. Fluctuation enhancement is associated with not only the random drift but also the clumping behavior (Eqs. (12) and (15)) in stationary state. It is found that random drift increases with decrease in steady-state density (Figs. 3 and 4). Such increase may be related to critical slowing-down (see Appendix A). Moreover, we find in the vicinity of extinction that the increase of random drift for a CP is more pronounced than that for a MFS. This is caused by the spatial distribution of individuals: when m approaches the
199
extinction point mc , the degree of clumping of species X in stationary state is extreme (see Eq. (15)). When m2 ∼ mc , the CP system in dynamical process more or less has the property of critical behavior (extreme increase of variance Vs near the extinction threshold) in stationary state. For this reason, the enhancement may be clearly observed for CP. In the Appendix we attempt to explain why fluctuation enhancement is not clear in the MFS system. It should be emphasized that MFS ignores the spatial distribution of individuals. To date, fluctuation enhancement has been observed in some physical systems (Kubo et al., 1973; Tsuchiya and Horie, 1985). In such systems, values of order parameters (in this case “density”) are known to increase from zero to a positive value. In contrast, we have shown just the opposite case: even though density is decreased, the fluctuation enhancement occurs. In general, if the density would increase, the variance V(t) often increased. An example is the exponential growth; the density difference at initial condition always increases. Our results suggest that both critical slowing-down and the clumping behavior of individuals may be essential to fluctuation enhancement. Almost all species form a spatially clumped distribution. Especially, when a species becomes endangered, the degree of clumping usually increases. This is due to the inherent nature of biological species: an offspring is produced in the neighborhood of parents. Here we discuss fluctuation enhancement. The definition should be made clearer, because increases of variance, V(t), have been observed for other reasons. For example, it has been shown that the variance is large for the stochastic logistic equation when R = 0 and K takes a large value (Renshaw, 1991). This large variance originates as a result of Brownian drift. In contrast, the fluctuation enhancement presented here is a new source of variability. We have shown that fluctuation enhancement is due to the spatial characteristics of the process; only local interactions give rise to increased fluctuations. The clumping behavior of individuals (Eq. (15)) near the extinction point appears to play an important role.
Fig. 6. Results of 100 perturbation experiments for the contact process (CP). The solid line represents A(t) and the points represent V(t). The critical value of the mortality parameter (point of phase transition) is mc ∼ 1.21. The death rate m is changed from m1 = 0.5 to m2 . (a) m2 = 1.3, (b) m2 = 2.0, (c) m2 = 5.0, (a) clearly exhibits this fluctuation enhancement (increase of V(t)). When m2 takes large values, no fluctuation enhancement occurs as can be seen in (b) and (c).
200
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
Finally, we propose the following experiment. Consider the cultivation of bacteria. Culture bacteria in two nutrient media: agar and a fluid medium. The growth of bacteria in the agar medium parallels the CP in which the growth process (2a) occurs locally: offspring are reproduced at a sites adjacent to their parents. In the fluid medium, the solution is continuously stirred, so that offspring are randomly dispersed in the medium. Thus, interactions in this medium correspond to a MFS. The reaction (2b) corresponds to the death processes of bacteria in the actual experiment; (2b) is realized using UV light or radioactive rays. The parameter m represents the intensity of light (radiation). Our results predict that there are a variety of different processes to extinction in the case of agar medium (CP) and that the stirred liquid medium will show a smooth transition to extinction.
Acknowledgements
x(t) − K ∝ exp(−(mc − m)t).
Appendix A. Mean-field theory We prove the critical slowing-down, and explain why fluctuation enhancement almost disappears in mean-field theory. Assume that the total number of sites is infinitely large (L → ∞); the dynamic equation for MFS becomes the logistic Eq. (9). The steady-state solution can be obtained by setting the time derivative in (9) to zero. The non-trivial solution (x = K) for 0 < m < 2 is stable, and phase transition occurs at mc = 2. Consider the density x(t) at time t, and assume that the system initially stays near the stationary state (x(0) ∼ K). Expand x around the steady-state density (K), x(t) = K + X(t), where X takes small values (X 1). Inserting this equation into (9) and collecting first order terms to obtain (A.1)
(A.2)
From this equation, the relaxation time is defined by 1/(mc − m). Hence, we prove that the relaxation time diverges as m → mc . This is called the critical slowing-down: the variance (V) in stationary state may diverge near m = mc . Next, we explain why fluctuation enhancement almost disappears in the mean-field limit. It is known in the field of physics that there is a close relation between critical slowing-down and fluctuation enhancement (Kubo et al., 1973; Tsuchiya and Horie, 1985). We carry out the same perturbation experiment N times, with initial densities given by xi (0) (i = 1, 2, . . . , N). After expanding and rearranging terms in the definition of V(t), we have 2 1 1 V(t) = xi (t)2 − xi (t) N N i
The authors are grateful to Professors Jin Yoshimura for valuable discussions. They also thank to referees for many useful comments and criticisms.
˙ = −(mc − m)X. X
From this equation it follows that
1 = [xi (t) − xj (t)]2 . N
i
i>j
Taking the derivative with respect to time, t, it follows that 2 V˙ = 2 (xi − xj )(˙xi − x˙ j ). (A.3) N i>j
The dynamics for t > 0 is assumed to be determined by the logistic Eq. (9), and following the perturbation, both parameters R and K are negative (m2 > mc ). Let xi (t) and xj (t) be two solutions of (9) whose initial values are slightly different from each other. Without loss of generality, we assume that xi (0) > xj (0). From (9), we have xi + x j x˙ i − xj = R(xi − xj ) 1 − . (A.4) K Hence, the right-hand side of (A.4) is always negative, and the difference xi (t) − xj (t) decreases as the time proceeds. It is therefore proved from (A.4) that the variance V(t) always decreases. Strictly speaking, this proof is insufficient, since the dynamics for a finite lattice (finite value of L) is described by a stochastic differential Eq. (1). Nevertheless, the logistic Eq. (9)
Y. Itoh et al. / Ecological Modelling 174 (2004) 191–201
explains well how (i) critical slowing-down occurs, and why (ii) enhancement of fluctuation almost disappears in the mean-field limit. References Abarca-Arenas, L.G., Ulanowicz, R.E., 2002. The effects of taxonomic aggregation on network analysis. Ecol. Model. 149, 285–296. Binder, K. (Ed.), 1979. Applications of Monte Carlo Methods in Statistical Physics. Springer-Verlag, Berlin. Durrett, R., 1988. Lecture Notes on Particle Systems and Percolation. Wadsworth and Brooka/Cole Advanced Books & Software, California. Durrett, R., Levin, S.A., 1994. The importance of being discrete (and spatial). Theor. Popul. Biol. 46, 363–394. Harada, Y., Iwasa, Y., 1994. Lattice population dynamics for plants with dispersing seeds and vegetative propagation. Res. Popul. Ecol. 36, 237–249. Harris, T.E., 1974. Contact interaction on a lattice. Ann. Prob. 2, 969–988. Katori, M., Konno, N., 1991. Upper bounds for survival probability of the contact process. J. Stat. Phys. 63, 115–130. Krivtsov, V., 2001. Study of cause-and-effect relationships in the formation of biocenoses: their use for the control of eutrophication. Russian J. Ecol. 32, 230–234. Krivtsov, V., Corliss, J., Bellinger, E., Sigee, D., 2000. Indirect regulation rule for consecutive stages of ecological succession. Ecol. Model. 133, 73–82. Kubo, R., Matsuo, K., Kitahara, K., 1973. Fluctuation and relaxation of macrovariables. J. Stat. Phys. 9, 51–96. Kubo, T., Iwasa, Y., Furumoto, N., 1996. Forest spatial dynamics with gap expansions: total gap area and gap size distribution. J. Theor. Biol. 180, 229–246. Lande, R., 1995. Mutation and conservation. Conserv. Biol. 9, 782–791. Liggett, T.M., 1985. Interacting Particle Systems. Springer-Verlag, New York. Marro, J., Dickman, R., 1999. Nonequilibrium Phase Transition in Lattice Models. Cambridge University Press, Cambridge. Matsuda, H., Ogita, N., Sasaki, A., Sato, K., 1992. Statistical mechanics of population: the lattice Lotka–Volterra model. Prog. Theor. Phys. 88, 1035–1049.
201
May, R.M., 1973. Stability and Complexity in Model Ecosystems. Princeton University Press, Princeton. Nakagiri, N., Tainaka, K., Tao, T., 2001. Indirect relation between species extinction and habitat destruction. Ecol. Model. 137, 109–118. Nisbet, R.M., Gurney, W.S.C., 1982. Modelling Fluctuating Populations. Wiley, New York. Nowak, M.A., Bonhoeffer, S., May, R.M., 1994. More spatial games. Int. J. Bifurcation Chaos 4, 33–56. Paine, R.T., 1966. Food web complexity and species diversity. Am. Nat. 100, 65–75. Pimm, S.L., 1993. Discussion: understanding indirect effects: is it possible? In: Kawanabe, H., Cohen, J.E., Iwasaki, K. (Eds.), Mutualism and Community Organization. Oxford University Press, Oxford, pp. 199–209. Renshaw, E., 1991. Modelling Biological Populations in Space and Time. Cambridge University Press, Cambridge. Sato, K., Matsuda, H., Sasaki, A., 1994. Pathogen invasion and host extinction in lattice structured populations. J. Math. Biol. 32, 251–268. Schmitz, O.J., 1997. Press perturbations and the predictability of ecological interactions in a food web. Ecology 78, 55–69. Schoener, T.W., 1993. On the relative importance of direct vs. indirect effects in ecological communities. In: Kawanabe, H., Cohen, J.E., Iwasaki, K. (Eds.), Mutualism and Community Organization. Oxford University Press, Oxford, pp. 365–411. Simberloff, D., Stiling, P., 1996. How risky is biological control? Ecology 77, 1965–1974. Soulé, M.E., 1987. Viable Populations for Conservation. Cambridge University Press, Cambridge. Tainaka, K., 1988. Lattice model for the Lotka–Volterra system. J. Phys. Soc. Jpn. 57, 2588–2590. Tainaka, K., 1994. Intrinsic uncertainty in ecological catastrophe. J. Theor. Biol. 166, 91–99. Tilman, D., Downing, J.A., 1994. Biodiversity and stability in grasslands. Nature 367, 363–365. Tsuchiya, Y., Horie, S., 1985. Evolution process of the Williams domain in a nematic liquid crystal. J. Phys. Soc. Jpn. 54, 1–4. Wootton, J.T., 2002. Indirect effects in complex ecosystems: recent progress and future challenges. J. Sea Res. 48, 157–172. World Commission on Dams, 2000. Dams and Development: A New Framework for Decision-Making. Earthscan, London. Yodzis, P., 1988. The indeterminacy of ecological interactions as perceived through perturbation experiments. Ecology 69, 508– 515.