Model parameters

## 1. Introduction

Our world is changing at an unprecedented rate and we need to understand the nature of the change and to make predictions about the way in which it might affect systems of interest. In order to address questions about the impact of environmental change, and to understand what might be done to mitigate the predicted negative effects, ecologists need to develop the ability to project models into novel, future conditions. This will require the development of models based on understanding the (dynamical) processes that result in a system behaving the way it does. The majority of current ecological models are excellent at describing the way in which a system has behaved, but they are poor at predicting its future state, especially in novel conditions [1].

One fruitful model strategy in ecology has been the “biology-as-physics” way to approach ecosystems, i.e. setting up simple equations from which they could obtain precise answers [2]. This tradition of modelling has a long history that can be traced back to the work of the mathematical physicist Vito Volterra, who used simple differential equations to examine trends in prey and predatory fish populations in the Adriatic [3]. Such models have the advantage of mathematical tractability; often they can be solved analytically to give precise answers and can be easily interrogated to determine the sensitivity of the model to its parameters [1]. On the other hand, it seems obvious that simple models will never accurately reflect any particular system since they lack realism.

A common simplification of these mathematical ecological models is that they rely on the well-mixed assumption or, in the physics parlance, the mean-field (MF) approximation. It is well known that the MF assumption can simplify a complex *n*-species system by replacing all interactions for any one species with the average or effective interaction strength. While the MF assumption seems reasonable in the case of, for example, plankton dynamics, it seems hardly appropriate for other situations. One of such case is when there is spatial heterogeneity, produced for instance by a gradient in some parameter that controls the dynamics. This is the situation in the first problem that we will analyze here. Another case if in which MF breaks down is for example for sessile species whose individuals by definition may only interact with others in a limited neighborhood provided that their niches overlap. The latter is the case for forest trees (the second system we will analyze here). Relaxing the MF assumption requires us to take into account the extent to which the strength of interactions among many species changes with the relative distances between individuals in space [4]. Many features of ecological dynamics such as the patterns of diversity and spatial distributions of species can be fundamentally changed when abandoning the MF assumption [4, 5]. A common way of relaxing the MF assumption is by formulating a spatially explicit individual-based model (IBM), or multi-agent system, whose straightforward implementation is by means of a cellular automaton.

Therefore, in order to go beyond MF and taking into account spatial heterogeneity, we will consider the application of cellular automata (CA) to two different important problems in Ecology where the space introduces important information or it simply cannot be neglected.

The first problem is about getting spatio-temporal early warnings of *catastrophic regime shifts* in ecosystems [6]. There is increasing evidence that ecosystems can pass thresholds and go through regime shifts where sudden and large changes in their functions take place. An example of such a regime shift is lakes that suddenly switch from clear to turbid water due to algae blooms. These blooms are connected to *eutrophication* i.e. the overenrichment of aquatic ecosystems with nutrients, principally phosphorous [7]. This is a widespread environmental problem because when it occurs, many of the ecosystem services which humans derive from these systems, such as fisheries and places for recreation, can be lost. Furthermore, it is often difficult, costly and impossible to reverse these changes once a certain threshold has been crossed. This is why early warnings of these shifts are so important to ecosystem management. In the case of MF models the rising of the temporal variance for the nutrient concentration was shown that it works as an early warning signal [8]. Later on it was shown that in many cases if one takes into accounts explicitly the space the spatial variance provides an even earlier early warning [9,10]. Thus, in section 2, I will present a cellular automaton that models a lake as a square lattice with the phosphorous concentration as the dynamical variable defined on lattice cells.

The second problem represents a major challenge in ecology: to understand and predict the organization and spatial distribution of biodiversity using mechanistic models. Ecologists have long strived to understand the distribution of relative species abundance (RSA) as well as the species–area relationships (SAR) in different communities [11, 12]. These metrics provide critical information that together can help uncover the forces that structure and maintain ecological diversity [13, 14]. Competition between species is one of the main mechanisms proposed to explain the observed RSA and SAR in different communities. In its basic form, the dynamics of competition-driven communities result from the degree to which species have overlapping niches because of their sharing of similar resource needs [15]. Hence in section 3 I will introduce a simple microscopic spatially explicit model to address the biodiversity distribution and spatial patterns observed in natural communities. This model, in terms of local competitive interactions between sessile species (for example trees), requires both niche overlap and spatial proximity.

## 2. *Mendota* cellular automaton: catastrophic shift in lakes and its spatial early warnings

### 2.1. The Cellular Automaton

The *Mendota* cellular automaton (MCA) described in this section is proposed to analyze the catastrophic transition from clear to turbid water and discuss possible early warning signals to this shift. The model will be a spatial version of a mean field model introduced by Carpenter [16] consisting of the three differential equations for phosphorus densities in soil (*U*), in surface sediment (*M*) and in water (*P*). In fact *P* and *M*, who are attached to the lake, are local spatial variables, while *U*, which describes the surroundings of the lake is taken as a global (non-spatial) variable. The evolution equations for *U*(*t*), *P*(*x,y;t*) and *M*(*x,y;t*) are:

where

Parameters of the model are defined in Table 1. We have also included diffusion with a diffusion coefficient *D* = 0.1. Another modification, in order to incorporate the effect of mechanical stirring of the lake waters (wind, currents, animals) is that we consider that at each time *t*, *a*(*t*) ≡ *cU*(*t*) *a*(*x,y;t*) fluctuates locally, from point to point, around its average global value *a*(*t*) in the interval [*a*(*t*) - Δ, *a*(*t*) + Δ ] where we have taken Δ = 0.125 and have verified that the results do not depend much on this value.

The lake is represented by a square lattice of *L*×*L* cells each one identified by its integer coordinates (*i,j*). Of course lakes of arbitrary shape could be studied by embedding them into a square lattice like the one above, with appropriate boundary conditions. Another approximation is that the system is two-dimensional, there is no depth. That is, on each cell there are two local variables assigned: *P*(*i,j*) and *M*(*i,j*). Therefore, equations (1) to (3) lead to the following CA synchronous update rules in discrete time, where now *t* represent the time measured in years:

Symbol | Definition | Units | value ^{} |

b c | Permanent burial rate of sediment P P runoff coefficient | y_1 y_1 | 0.001 0.00115 |

F-H | Annual agricultural import minus export of P per unit lake area to the watershed | g_m_2y_1 | 0 or 13 |

h | Outflow rate of P | y_1 | 0.15 |

m | P density in the lake when recycling is 0.5 r | g_m_2 | 2.4 |

r | Maximum recycling rate of P | g_m_2y_1 | 0.019 |

q | Parameter for steepness of f(P) near m | Unitless | 8 |

s | Sedimentation rate of P | y_1 | 0.7 |

W0 | Nonagricultural inputs of P to the watershed before disturbance, per unit lake area | g_m_2_y_1 | 0.147 |

WD | Nonagricultural inputs of P to the watershed after disturbance, per unit lake area | g_m_2_y_1 | 1.55 |

### 2.2. Observables

Catastrophic shifts have characteristic fingerprints or ’wave flags’. Some of the standard

catastrophe flags are: modality (at least two well defined attractors), sudden jumps and a large or anomalous variance [17]. Therefore we will calculate the following corresponding observable quantities for the phosphorous density in water *P,* which is the relevant variable in our case.

1. The *spatial variance* of *P*(*x,y;t*), *σ* _{s }^{2 }, defined as

2. The *patchiness or cluster structure*.

The reason leading to the rise in *σ* _{s }is the onset of fluctuations in the spatial dependence of *P*(*x,y;t*) at a given time.

3. Similarly to the mean free model, the *temporal variance* *σ* _{t }^{2 }, which has been suggested by Brock and Carpenter [8] as an early warning signal. At an arbitrary point, say (*x,y*) = (0,0), is defined as for temporal bins of size τ.

### 2.3. Simulations and Results

Different simulations were run for 1500 years to illustrate changes of the system over time. The first 250 years of each simulation are a highly simplified representation of the history of Lake Mendota. The remaining years of each simulation illustrate contrasting management policies for soil phosphorus.

Following [16], simulations were initiated at stable equilibrium values calculated with *F* = *H* = 0 and *W* for undisturbed conditions. These represent presettlement conditions, and were maintained for years 0–100. For years 100–200, *W* was changed to the value for disturbed conditions (*W*D, Table 1), representing the advent of agriculture in the region. For years 200–250, *F* and *H* were increased to the values estimated for a period of intensive industrialized agriculture in the Lake Mendota watershed, and *W* was maintained at *W* _{D} (Table 1). That is, four different stages are considered:

years 0-100

*W=W*_{0 }=0.147,*F*-*H =*0*.*years 101-200

*W*=1.55,*F*-*H*= 0*.*years 201-250

*W*=1.55,*F*-*H*= 31.6-18.6=13*.*after year 250

*W*=1.55,*F*-*H*= 0*.*

After year 250, the simulations were different. Simulation 1 represents management to balance the phosphorus budget of agriculture. In this simulation, after year 250 *F* = *H* and *W* is maintained at *W* _{D}.

As one can see from Figure 1, the 1500 years history of lake Mendota produced by the Mendota CA is similar to the Carpenter's 2005 non spatial model [16]: a gradual first shift for *P* that start at year 250 and a second more drastic shift between years 440 and 485. Notice that the *σ* _{s }produced by the BCI CA reaches at year 425 three times the constant value it had along the first 400 years while it reaches its maximum value at year 467. The anticipation of the early warning depends on the convention chosen to establish when the shift actually happens. If we adopt the *Maxwell convention* [17] i.e. the shift coincides with the time when the variance reaches its maximum value [9]. If, on the other hand, we chose the Delay Convention the shift occurs when the system stabilizes in its final attractor i.e. at year 485 (see Figure 1). Therefore *σ* _{s }provides an early warning of the coming shift that varies between 40 and 60 years.

Concerning the practical use of *σ* _{s }as an early warning, Figure 2 shows that the time at which the signal becomes appreciable depends on the number of points of the grid (or lattice size *L*). As expected, there is a trade-off between anticipation and cost: the larger the number of points involved in the statistics, the sooner will be this early signal, until, for grids above 100×100 = 10,000 points, the curve doesn’t change very much. For instance, for a 20×20 grid the early warning becomes noticeable 20 years later than when computed for the entire lattice.

Figure 3 shows color maps representing different configurations when the lake moves towards a catastrophic shift from clear to turbid water. Figure 3(A) to (C)correspond, respectively, to pictures of the lake at t = 250 yrs (when still *F=H*), at t = 437 yrs (30 years before the peak in σ_{s} ) and at t = 467 yrs (just at this peak of σ_{s}). For *t* = 467 typically there occur patches covering a wide scale of sizes (Figure 3(C)).

Finally let us have a look to *σ* _{t }and compare it versus *σ* _{s }. Figure 4 shows both variances. The temporal variance exhibits a temporal delay of around 40 years leaving much less margin for eventual remedial management actions. This delay can be easily understood since *σ* _{t }employs data in times where the fluctuations are still small.

### 2.4. Discussion and some Remarks

We have proposed a CA that describes the evolution of the P concentration, which dominates the eutrophication process in lakes, with its parameters calibrated for Lake Mendota. We have considered different possible early warnings for eutrophication shifts in lakes. In the case of *σ* _{s }, by measuring samples of *P* on a grid of points on the lake surface, it was found that a grid containing few points might be sufficient for the purpose of extracting an appropriate signal, and that a significant growth in *σ* _{s }could serve as an early warning of an imminent transition. The spatial variance appears to have an advantage over the temporal one, as *σ* _{t }is delayed with respect to *σ* _{s }.

When studying the origin of the rise in *σ* _{s,} we found that it is connected to the appearance of spatial patterns, in the form of clusters of clear and turbid water. The spatial patchiness under routine conditions is quite different from spatial patchiness near a regime shift. This is because in the last case in some places the clear water state is realized while in others the turbid water state occurs. In mathematical terms, close to the regime shift there are two competing attractors or alternative states (clear and turbid water) one occurs in a given set of cells and the other in the set of remaining cells. This explains the spatial fluctuations in the concentration of P. Under routine conditions, either before or after the catastrophic shift, only one attractor remains.

We then conclude that the visualization of the onset of those patches, for example by aerial or satellite imaging of the lake surface, may be an effective way of anticipating an eutrophication transition.

However a note of caution is required, in thermally stratified real lakes two main layers can be distinguished. The hypolimnion is the dense, bottom layer of water in thermally-stratified lakes, which is isolated from surface wind-mixing. The epilimnion is top-most layer, occurring above the deeper hypolimnion and being exposed at the surface, it typically becomes turbulently mixed as a result of surface wind-mixing. On the one hand, the flow of phosphorus in recycling occurs first from sediment to and also from decomposition in the hypolimnion. Dissolved P in the hypolimnion tends to be rather well-mixed and homogeneous. Then vertical mixing of hypolimnetic and epilimnetic water occurs from time to time (roughly 10-15 day return time [18]. Horizontal dispersion of dissolved P is fast after a vertical mixing event, within a day or so, smoothing over the patchiness structure in the epilimnion. As a consequence patches could never survive during a year, which is the unit of time in our model. On the other hand, some lakes show important spatial differences in terms of total phosphorus related to the spatial distribution of submerged plants or cyanobacteria blooms. The biomass tends to accumulate in several zones, by the hydrodynamic and wind actions, and determine strong differences of algal biomass as well as total phosphorous [19]. In any case, since the spatial model we are analyzing is two-dimensional and there is no depth (the hypolimnion and the epilimnion are collapsed to a single layer), the above processes are not included. Indeed it seems that the patches produced by the spatial model, instead of long lasting structures, should be interpreted as frequent popping in and out clusters. This dynamic patchy pattern changes many times during the course of one year.

It is worth to remark that the quantitative details of our conclusions depend on the choice of parameter values employed in our model. We have verified that the qualitative behavior of our results do not depend strongly on those values. Furthermore our model as presented in this work is schematic in the sense that the quality of the lake’s water is dependent of a single parameter, the amount of Phosphorous in solution. In real cases other environmental factors might be playing an important role in the evolution, and the model and its predictions could be much more complicated. Nevertheless it appears that our main conclusions should hold in the more realistic case: spatial variance of critical quantities is an earlier signal than the temporal variance, and its associated cluster structure of the patterns formed in the eutrophication process could be the fastest detectable warning that a catastrophic change is about to occur.

## 3. *Barro colorado island* cellular automaton: species competition in physical and niche spaces

### 3.1. Niche Lotka-Volterra Competition Cellular Automaton

The model combines features from both niche-based models of interspecific competition [20] and from the neutral theory of biodiversity [21].

The CA considered in this section is based on the Lotka–Volterra competition equations for *N* _{sp} species [20]: *dns/*d*t* =*g* _{s }*n* _{s }(1 *−* _{r }*α* _{sr }*n* _{r }), where *s* and *r* run from 1 to *N*sp, *n* _{s }is the density of species *s*; *g* _{s }its maximum per capita growth rate, and the coefficients *α* _{sr }represent the competitive interaction of species *r* over species *s*. Estimating the interaction strengths among species is, however, far from trivial because of the large number of potential interactions in diverse natural communities, feedback effects and non-linearities in interaction-strength functions [22]. In fact the Lotka*-Volterra Competition Model* (LVCM) include *N* _{sp}×(*N* _{sp}-1) competition coefficients *α* _{sr }. In the case of a tropical forest, where there are typically some hundreds of different coexisting species, estimating the corresponding several thousand coefficients becomes definitely an impossible task. Thus community ecologists have resorted to different approaches to include these diverse interaction strengths into their dynamics. Following the commonly used MacArthur and Levins’ approach [23,20], we consider the simplest one-dimensional finite niche wherein the resource utilization function of each species *s* can be represented by a normal distribution *P* _{s }(*x*) = exp[*−*(*x − μ* _{s })^{2} */*2*σ* _{s }^{2}], with mean *μ* _{s }and a standard deviation *σ* _{s }that measures the niche width. Then for each pair of species *s* and *r*, the strength of its competition is determined by their niche overlap, i.e. the overlapping between *P* _{s }(*x*) and *P* _{r }(*x*) (see Figure 5). We call this model the *Niche Lotka-Volterra Competition Model* (NLVCM). Assuming a finite normalized niche (*x ∈* [0, 1]), the competition coefficients *α* _{sr }can thus be computed as:

*.*

On the other hand, inspired by neutral theory and to keep things as simple as possible, we assume two important simplifications. Firstly, that all species have the same maximum growth rate i.e. *g* _{s }= 1 and the same niche width i.e. *σ* _{s }= *σ*. The robustness of this simplification was recently analyzed for the NLVCM [24]. Therefore, each species’ growth rate depends both on its position in the niche axis and on the niche positions of the individuals of other species that make its competitive neighborhood at each time step. The (intrinsic) functional equivalence between species is consistent with the chosen normalization for the *α* _{rs }in (Fig. 5) that ensures that the matrix *α* is symmetric and allows it to be expressed as:

where ‘erf’ denotes the standard error rate function. Secondly, we consider that the community has a constant, finite size *N* [21] such that local reproductive and mortality events are tightly synchronized (see the used update rule below); this is equivalent to the ‘zero-sum assumption’ of Hubbell’s neutral model of biodiversity [21].

The model is defined as follows:

* Landscape. We work on a *L*× *L* = *N* square lattice with periodic boundary conditions to avoid border effects. Each lattice cell (i,j) is always occupied by one individual belonging to a certain species specified by its niche position *μ* _{s }(i,j), s =1,2,...,*N* _{sp}, where *N* _{sp} is the number of initially coexisting species.

* Initial conditions. Because the initial spatial distribution of individuals at BCI is of course unknowable, we run the model from 100 different initial spatial configurations of completely random assignation of species (i.e. corresponding to the maximum entropy).

* Dynamics and update rule. The dynamics is asynchronous. At each time step *t* a focal cell is chosen at random. The LVC dynamics is implemented by using the "copy the fittest" rule i.e. by replacing each focal individual by the individual having the maximum fitness in the Moore neighborhood (the focal individual and its eight surrounding neighbors), which eventually can be itself and thus no change occurs at the focal node.

We will compute in addition of RSA and SAR, aggregation properties.

### 3.2. Comparison with empirical data from tropical forests

In this subsection we will contrast the model against field data from a well-studied ecological community, the tropical forest at Barro Colorado Island (BCI), Panama [25,26,27]. Population structure and spatial patterns have been particularly important themes in the ecology of tropical trees [21]. The long-term research program coordinated by the Smithsonian Center for Tropical Forest Science [28] includes inventory data collected for the first time in 1982 and then every 5 years from 1985 to 2005, with measurements of all trees with diameter at breast height (dbh) *>* 1 cm for a plot of 50 ha in BCI [29].

*3.3. Fitting parameters and biodiversity indices*

We set the initial number of species *N*sp = 320. This is the maximum number of species found for BCI in the six censuses for 1982 and 1985, although this number steadily declined afterward from census to census until in 2005 the number of species recorded was 283 [29]. The number of individuals was kept fixed as *N* = 250,000 (*L* = 500) that is close to the maximum number of individuals observed in the set of BCI censuses. Therefore, the only free model parameters to be fitted are *σ*, controlling the intensity of interspecies competition, and the number of steps of simulation *t*.

This second parameter is needed because we are interested in transient states rather than in the steady state. The reason for this is simple: the set of censuses for tropical forests collected by the Center for Tropical Forest Science reveal that they are non-stationary systems: species richness decreases from one census to the next in *all* tropical forests and in many cases at a considerable pace [28]. Just to cite some examples: at BCI (Panama), roughly 3% of the number of species have been disappearing every 5 years; for Bukit Timah (Singapore), this percentage varies between 3% and 8% between consecutive censuses; for Edoro (Congo) 12% of the species disappeared between the 1994 census and the 2000 census, etc. The dynamics of the RSA between censuses in each forest are also far from being steady. It is noteworthy that most analyses of the BCI (and of other tropical forests) largely involve fitting separately several metrics of community structure and organization (e.g. SAR, RSA, indices of diversity) for each census rather than considering them as a dynamic sequence whose attributes needs be fitted together [25], [30],[31].

Moreover, the BCI CA reaches a steady state with four species after a long transient of the order of 3 billion time steps (data not shown). This a known result for MF [32] that was analyzed in detail in [24][33]. It is of course impossible to decide whether this steady state could be a reasonable result in a more or less distant future, and it is indeed absurd to run any ecological model for such a long time horizon without including the effect of speciation at this very long temporal scale. Rather the BCI CA aims to describe the transient corresponding to an entire set of censuses (rather than viewing them as separate snapshots as is generally done), which requires specifying the simulation time *t* as a parameter.

The decrease (although moderate) of Shannon entropy provides another way to verify the non-equilibrium nature of the BCI forest and this negative entropy production corresponds to the loss of species and the changes in community structure over time (see below).

We searched for the combination of *σ* and the number of time steps that produce the best agreement between the theoretical and the empirically observed biodiversity for the entire set of six BCI censuses. In addition to RSA we used two well known and often used indices of diversity, the Shannon entropy *S* and the Simpson index 1- *D* [34], to describe community structure for each census. Normalized versions of these indices are given by

where *fs* is the relative abundance of species *s*. For both indices, the greater the value, the larger the sample diversity.

The procedure to find the optimal *σ* and *t* is as follows:

First, for a given value of *σ*, each simulation is run until the entropy *S* becomes equal to the entropy corresponding to the BCI population distribution for the first census (1982), i.e. *S*1 = 0*.*694. It turns out that for most values of *σ* except for a narrow set of values centered around *σ* = 0*.*09, the RSA are clearly different from those of BCI in 1982. A useful representation in ecology to perform such comparison between theoretical and empirical RSA is provided by the commonly used *dominance–diversity* *curves*. These curves are obtained by ranking the species according to their population abundance and plotting the RSA (%) versus species rank. Specifically, we found that the theoretical dominance–diversity curve converges to the empirical one as *σ* approaches 0.09. The Simpson diversity index 1 *− D*1 = 0*.*949 observed for the first census served to narrow even further the range of values of *σ* around 0.09.

Second, for *σ* = 0*.*09, we look for the best fit to the entire sequence of six BCI censuses. We found that the optimal fit occurs for *t* _{1} = 37 million of time steps (Mts) for the first 1982 census and then there is a correspondence of 1 year *←→* 250 000 time steps. Therefore, the census years 1982, 1985, 1990, 1995, 2000 and 2005 correspond, respectively, to *t* = 37, 37.75, 39, 40.25, 41.5 and 42.75 Mts. The Barro Colorado cellular automaton (BCI CA) is then defined by *L* = 500, *N*sp = 320, *σ* = 0*.*09 and *t* _{1} = 37 Mts.

Figure 6 shows the close agreement between the *S* and 1 *− D* calculated for BCI and the predicted average values obtained for 100 simulations each starting from different initial conditions, for the sequence of six censuses. This figure also shows the MF results, i.e. numerical integration of LVC equations, indicating that this MF cannot simultaneously and accurately predict the values of both biodiversity indices.

In addition, the RSA curves for the MF approximation showed a poorer agreement with the empirical ones for each census (not shown).

*3.4. Dominance–diversity curves*

As in most biological communities, most trees species on BCI have few individuals: only five species *Hybanthus prunifolius*, *Faramea occidentalis*, *Trichilia tuberculata*, *Desmopsis* *panamensis* and *Alseis blackiana* individually had more than 5% of the total community size.

Figure 7 compares the predicted and observed dominance–diversity curves for the BCI 1995 and 2005 censuses. In both cases the predicted curve slightly underestimates the abundance (notice the logarithmic vertical scale) of rare species, i.e. those representing a RSA smaller than 0.05% which is a percentage well below 0.35% and 0.31%, the average values of RSA at BCI for 1995 and 2005, respectively (horizontal lines in figure 7). It turns out that for values of *σ* departing from 0.09 this agreement found for common species becomes worse.

The most abundant species that yield the BCI CA are those at both ends of the niche axis, i.e. *μ* _{s }0 or 1. The explanation is quite simple: species at the ends of the niche axis are exposed to less competition than those which depart from the ends and experiment competition from the two sides instead of from only one.

*3.5. Spatial patterns*

A prediction of BCI CA model is that local interactions in physical space introduce an interesting fact: individuals of different species close in niche space become spatially segregated. This is shown in Figure 8-A which corresponds to a 100 × 100 patch of the entire lattice depicting the spatial distribution of the three most abundant species at the right end of the niche axis (i.e. all with *μ* _{s }close to one) after 40 Mts. In order to contrast this against empirical data, Figure 8-B reproduces the spatial distribution found in 1982 for three species of quite similar trees with comparable population sizes, *Faramea occidentalis*, *Trichilia tuberculata* and *Alseis blackiana,* in a 141 m × 141 m [1] - patch of the BCI 50 ha plot [25,26,27]. Notice that these three species also show spatial segregation, although the empirical segregation seems to be lower than the theoretical one. A simple explanation for this is that the BCI CA is little too simplistic, only local replacement of species is taking into account. Indeed we have checked that by introducing a global migration parameter *m,* modeling dispersal of seeds by winds or birds, the agreement between theoretical and empirical segregation patterns considerably improves.

The message is then that, while individuals of these three species could potentially have strong competitive interactions because of similarities of their niches, their spatial segregation attenuates the strength of their interspecific competitive interactions whenever these are local. Spatial segregation is known to be a mechanism that allows the coexistence of competing species with similar niches even when space is introduced only implicitly in a MF model [35]. Furthermore, space is known to be important in the formation and maintenance of stable vegetation patterns [36,37]. Thus, the local competitive interactions might actually reinforce the *emergent neutrality* found in non-spatial competition models [38].

There are many options to quantify whether there is spatial segregation among these dominant species. One of them is to partition a square patch *L* _{P }of side multiple of 3 (e.g. *L* _{P }= 99) into 3×3 Moore neighborhoods and measure the fraction *f* _{c }of these (*L* _{P }/3)^{2} in which there is coexistence of species (at least two of the three species are present). Then by repeating the same procedure for exactly the same number of individuals belonging to each of these three species but *randomly* distributed, and calculating the fraction *fr* _{c }, the segregation index γ can be constructed as:

When γ < 1 (>1) there is (there is not) spatial segregation among the species. We found γ = 0.41, confirming then that these three dominant species, whose niches are tightly packed at the right end of the niche axis, are actually spatially segregated.

The SAR (average number of species in a sampled area) constitute one of the main metrics in spatial community ecology [14], [21]. To generate SAR, we proceed as in [34] by dividing the entire area into non-overlapping squares and rectangles and counting the number of species present in quadrats of 20 *×* 20, 25 *×* 20, 25 *×* 25, 50 *×* 25, 50 *×* 50, 100*×*50, 100*×*100, 250*×*100, 250*×*250 and 500*×*500. This procedure yields the mean number of species in each size quadrat that make up the SAR for each year. Figure 9 shows that after 39 Mts, the model produces a SAR that fits quite well to the one obtained from the 1990 BCI census [34] for large area plots (above 1 ha). However, as the area of plots decreases, the agreement with data worsens. This is expected, because we know that our model underestimates the abundance of the scarcest species as shown in the dominance-diversity curves (figure 7) whose presence becomes more and more significant as the plot sizes become smaller than in larger plots that typically contain more individuals of many different species. Another useful spatial metric is provided by species–individuals curves, which are obtained exactly as the SAR but using the mean number of individuals in quadrats of each area instead of sampled area in the horizontal axis. Again, there is a qualitatively good agreement between model predictions and field data (not shown here).

### 3.6. On the importance of the locality of the competition and possible extensions to other realms.

The BCI CA is a model of local competition that jointly considers organisms in both physical and niche space. It must be viewed as a minimalistic model of local competition, with basically only two parameters: *σ*, the species’ niche width controlling the intensity of the interspecific competition that drives the spatiotemporal dynamics, and the simulation time *t* (since it is focused on transients rather than in the equilibrium state). This ‘microscopic’ model of local competition can predict with reasonable accuracy the dynamic sequence of patterns of community structure, species-packing and the spatial distribution of forest species. Interestingly, the model shows that when the well-mixing assumption underlying the MF approximation breaks down, species that are clumped in niche space appear spatially segregated. Including space is well known to allow the long-term coexistence of strongly competing species and to permit the formation of stable patterns [5],[35],[37].

## 4. Conclusions

Each one of the two CA presented here address an open problem in ecology, namely how to anticipate catastrophic shifts in ecosystems and understanding the forces that shape community biodiversity, respectively. Both are fundamental scientific questions that go beyond ecology and require a major interdisciplinary effort and would have a significant impact on ecosystem management and conservation. The two CA represent minimal spatially explicit ecological models, that could be called "physicist's models" since they include the minimal number of adjustable parameters.

First, let us consider some guidelines for future work concerning the problem of eutrophication of lakes and the more general problem of the development of early warning signals of catastrophic shifts. The eutrophication is a result of a combination of natural processes and human impacts [39]. CA like Mendota could serve to disentangle the effect of these natural processes, like climate change [40], from anthropic effects. This could be done using paleolimnological analysis to estimate and connect historic input and run-off rates with temperature and humidity indices [41]. Another interesting extension of this CA could be from 2D to 3D, to take into account the depth, in order to describe the flow of phosphorus in recycling from sediment and also from decomposition in the hypolimnion. Going beyond lakes, the analysis of spatially explicit models is relevant, for example, to understand phenomena like clumping and spatial segregation in plant communities [42]. It was shown that vegetation patches, which have been extensively studied for arid lands, can be approached as a pattern formation phenomenon [43]-[44]. It has been hypothesized that this vegetation patchiness could be used as a signature of imminent catastrophic shifts between alternative states [36]. Evidences that the patch-size distribution of vegetation follows a power law were later found in arid Mediterranean ecosystems [37]. This implies that vegetation patches were present over a wide range of size scales, thus displaying scale invariance. It was also found that with increasing grazing pressure, the field data revealed deviations from power laws. Hence, the authors proposed that this power law behavior may be a warning signal for the onset of desertification. These spatial early warnings complement temporal ones like the variance of time series.

Second, in the case of tropical forests, these are extremely diverse ecosystems -many of which contain in 50 ha sample plots more tree species as occur in all of US and Canada combined- for which a huge volumes of data (covering 47 permanent plots in 21 Countries with 4.5 million trees of 8,500 species [28]) has been accumulated throughout the last thirty years. This mega diversity, together with the abundance of available data, makes tropical forests a paradigm for research on the interdisciplinary field of complex systems dynamics. The BCI CA parameters could be easily calibrated for modelling other tropical forests different from BCI and contrast with empiric data. Future natural extensions of this CA that are worth considering are the inclusion of a migration rate, due to animals or wind propagating seeds beyond local neighborhoods, and noise. These two factors would add a dose of stochasticity, making the model more realistic. Considering to other realms, beyond ecology, the type of competition of BCI CA could be applied to understand other natural as well as non-natural communities. For example in the case of social and economic systems [45], communities formed by different kinds of interacting firms [46] -for example stores, banks, restaurants, etc- wherein territorial competition between heterogeneous individuals (i.e. occupying different niche positions) occurring in distinct local neighborhoods is a key factor controlling their dynamics. Of particular interest would be study the role of space in the winner-takes-all markets [47], in which in principle slight differences in performance of the firms can lead to enormous differences in reward. Long familiar in sports and entertainment, this payoff pattern has increasingly permeated law, finance, fashion, publishing, and other fields [48].

## Acknowledgments

I thank financial support from ANII (SNI), Uruguay. I am indebted to Steven Carpenter and Néstor Mazzeo for discussion on the material presented in section 2 and to Pablo Inchausti since the material of section 3 is based on a recent paper we wrote [49].

## Notes

- Since we use a sqaure 500 ( 500 lattice CA to represent the rectangular 500 m ( 1000m BCI plot, it turns that the lattice spacing corresponds to (2 ( 1.41 m.