Statistical mechanics of ecological systems: Neutral theory and beyond

I. Introduction

“It is interesting to contemplate an entangled bank, clothed with many plants of many kinds, with birds singing on the bushes, with various insects flitting about, and with worms crawling through the damp earth, and to reflect that these elaborately constructed forms, so different from each other in so complex a manner, have been all produced by laws acting around us.” In this celebrated text from the Origin of Species (), Darwin eloquently conveys his amazement for the underlying laws of nature: despite the striking diversity of shapes and forms, it exhibits deep commonalities that have emerged over wide scales of space, time, and organizational complexity. For more than 50 years now, ecologists have collected census data for several ecosystems around the world from diverse communities such as tropical forests, coral reefs, plankton, etc. However, despite the contrasting biological and environmental conditions in these ecological communities, some macroecological patterns can be detected that reflect strikingly similar characteristics in very different communities (see Table ). This suggests that there are ecological mechanisms that are insensitive to the details of the systems and that can structure general patterns. Although the biological properties of individual species and their interactions retain their importance in many respects, it is likely that the processes that generate such macroecological patterns are common to a variety of ecosystems and they can therefore be considered to be universal. The question then is to understand how these patterns arise from just a few simple key features shared by all ecosystems. Contrary to inanimate matter, living organisms adapt and evolve through the key elements of inheritance, mutation, and selection.

TABLE I.

Descriptions of the Macroecological patterns.
Macroecological pattern	Definition
Alpha diversity:	The number of species found in a given area, regardless of their abundance and spatial distributions. Sometimes it is more appropriate to assess the diversity of a local community by taking into account their abundances as well. For this purpose two alternative indexes have been used, the Shannon (H) and Simpson (D) indexes ().
Beta diversity:	The probability that two individuals at distance $𝑟$ belong to the same species (i.e., they are conspecific). Several alternative definitions have been used in the literature, including the pair (or two-point) correlation function (PCF), the Sørensen similarity index, or the Jaccard similarity index. However, the general purpose of beta diversity is to describe the turnover of multiple species in space.
Pair correlation function (PCF):	The correlation in species’ abundance between pairs of samples as a function of their distance.
Relative species abundance (RSA):	The probability that a species has $𝑛$ individuals in a given region. When multiplied by the total number of species in the region, this gives the number of species with $𝑛$ individuals (see Sec. ). This is sometimes called the species-abundance distribution (SAD).
Species (time) turnover distribution (STD):	The probability density function that the ratio of the future to the current population sizes of any species has a value $𝜆$ in a given ecological community (see Sec. ).
Persistence or lifetime distribution (SPT):	The probability density function of the time interval between the emergence and local extinction of any species within a given area (see Sec. ).
Species-area relationship (SAR):	The function that relates the mean number of species $𝑆$ to the area $𝐴$ they live in. On a relatively large range of spatial scales it is well approximated by a power law $𝑆 (𝐴) = 𝑐 𝐴 𝑧$ , where $0 < 𝑧 < 1$ (typical values of $𝑧$ are 0.2 or 0.3). It has been shown that the SAR in the log-log scale has three qualitatively different behaviors from local to continental spatial scales: approximately linearlike at very small and very large scales, and power-law-like for intermediate scales. This is referred to as a triphasic SAR (see Sec. ).
Endemic-area relationship (EAR):	The mean number of species that are present in the area $𝐴$ but not outside it (see Sec. ).

This fascinating intellectual challenge fits perfectly into the way physicists approach scientific problems and their style of inquiry. Statistical physics and thermodynamics have taught us an important lesson that not all microscopic ingredients are equally important if a macroscopic description is all one desires. Consider, for example, a simple system like a gas. In the case of an ideal gas, the assumptions are that the molecules behave as pointlike particles that do not interact and that only exchange energy with the walls of the container in which they are kept at a given temperature. Despite its vast simplifications, the theory yields amazingly accurate predictions of a multitude of phenomena, at least in a low-density regime and/or at not too low temperatures. Just as statistical mechanics provides a framework to relate the microscopic properties of individual atoms and molecules to the macroscopic or bulk properties of materials, ecology needs a theory to relate key biological properties at the individual scale, with macroecological properties at the community scale. Nevertheless, this step is more than a mere generalization of the standard statistical mechanics approach. Indeed, in contrast to inanimate matter, for which particles have a given identity with known interactions that are always at play, in ecosystems we deal with entities that evolve, mutate, and change, and that can turn on or off as well as tune their interactions with partners. Thus the problem at the core of the statistical physics of ecological systems is to identify the key elements one needs to incorporate in models in order to reproduce the known emergent patterns and eventually discover new ones.

Historically, the first models defining the dynamics of interacting ecological species were those of Lotka and Volterra, which describe asymmetrical interactions between predator-prey or resource-consumers systems. The Lotka and Volterra equations have provided much theoretical guidance. For instance, developed a model for studying interactions among consumers which exploit common resources. By making use of different time scales, he showed how resources can be included in the Lotka and Volterra equations. He also derived the formulas for the competition matrix, which provided hints about how it can be measured empirically (). Seemingly, this suggested a viable way to measure niche overlap between species. Soon after, studied the problem of how similar competing species can be and yet coexist in a community. According to Gause’s competitive exclusion principle (), two species cannot occupy the same niche in the same environment for a long time (see Table ). They found that environmental fluctuations limit niche overlap and therefore species’ similarity. All these studies, and many other variations and generalizations, provided a robust theoretical basis for understanding ecosystems. Their limitations were identified as more research was carried out (). These approaches share the idea that species can be well described by deterministic models which are shaped by the fundamental concept of niche, which, however, can be appropriately defined only a posteriori. This theoretical focus is not eminently suited for studying a wide range of empirical patterns. In fact, all such models have several drawbacks: (1) They are mostly deterministic models and often do not take into account stochastic effects in the demographic dynamics (). (2) As the number of species in the system increases, they become analytically intractable and computationally expensive. (3) They have a lot of parameters that are difficult to estimate from ecological data or experiments. (4) It is very difficult to draw generalizations that include spatial degrees of freedom. And (5) while time series of abundance are easily analyzed, it remains challenging to analytically study the macroecological patterns they generate and thus, their universal properties.

TABLE II.

Ecological glossary.
Ecological terminology	Definition
Trophic level:	The set of all species belonging to the same level in the food chain. Individuals of species belonging to higher levels feed upon those in the lower ones, while individuals belonging to the same trophic level compete for the same pool of resources. Neutral theory is an ecological theory for species in one specific trophic level, whereas Niche theory (see below) can also deal with species at different trophic levels.
Gause’s competitive exclusion principle:	A pair of species cannot stably coexist if they feed upon exactly the same resources under the same environmental conditions ().
Niche theory:	Species can stably coexist in an ecological community if their characteristics (or traits) allow them to specialize on one particular set of resources or environment conditions (niches) in which they are superior to their competitors. In other words, they occupy different niches (). Such niche separation is deemed to enhance trade-offs and facilitate coexistence, even though there is no a priori method to identify the correct niches that favor coexistence. The underlying rationale of the theory is Gause’s exclusion principle, and the classical model is the Lotka-Volterra set of differential equations. In contrast to niche theory, neutral theory claims that niche differences are not essential for coexistence.
Species-level neutral models:	These models assume that all species are equivalent as they all have the same probability of immigration, extinction, and speciation. The only state variable of these models is the number of species in a community and, thus, it cannot predict the distribution of the population sizes across species (see RSA) ().
Individual-level neutral models:	These models assume that all species are equivalent at the individual level, having the same birth, death, immigration, and speciation rates regardless of their identity. Therefore, although species can have different abundances, they are competitively equivalent. The state variable of these models is the population of any species in a region and therefore, such models can be used to understand species richness and abundance ().
Symmetric models:	Any model whose outcomes are invariant when exchanging species identities. The family of symmetric models is larger than the neutral ones: for instance, effects of environmental stochasticity or density dependence of individuals on the per-capita birth and death rates can be accommodated in symmetric (but not neutral) models. Some do not make a distinction between symmetric and neutral models.
Dispersal limited process:	Any process that constrains offspring to disperse in the vicinity of its parents.
Multispecies voter model (MVM) with speciation:	A spatially explicit neutral model in which each individual is located in a regular lattice and belongs to a species. Given a spatial configuration at time $𝑡$ , the configuration at time $𝑡 + 1$ is obtained as follows: an individual is chosen at random and is replaced by a copy of one of its nearest neighbors (chosen at random) with a probability $1 - 𝜈$ , or by an individual belonging to a new species (not already present in the lattice) with probability $𝜈$ . $𝜈$ is the speciation parameter, although it may incorporate immigration effects.

A pioneering attempt to explain macroecological patterns as a dynamic equilibrium of basic and universal ecological processes—and that also implicitly introduced the concept of neutrality in ecology—was made by in the famous monograph titled “The theory of island biogeography.” In this work, they proposed that the number of species present on an island (and forming a local community) changes as the result of two opposing forces: on the one hand, species not yet present on the island can reach the island from the mainland (where there is a metacommunity); and on the other hand, the species already present on the island may become extinct. MacArthur and Wilson’s model implies a radical departure from the then main current of thought among contemporary ecologists for at least three reasons: (1) Their theory stresses that demographic and environmental stochasticity can play a role in structuring the community as part of the classical principle of competitive exclusion. (2) The number of coexisting species is the result of a dynamic balance between the rates of immigration and extinction. (3) No matter which species contribute to this dynamic balance between immigration and extinction on the island, all the species are treated as identical. Therefore, they introduced a model that is neutral at the level of species (see Table ), even though they did not think of it as such.

Just a few years later, the ecologist H. Caswell proposed a model in which the species in a community are essentially a collection of noninteracting entities and their abundance is driven solely by random migration and immigration. In contrast to the mainstream vision of niche community assembly, where species persist in the community because they adapt to the habitat, Caswell stressed the importance of random dispersal in shaping ecological communities. Although the model was unable to correctly describe the empirical trends observed in a real ecosystem, it is important because it pictured ecosystems as an open system, within which various species have come together by chance, past history, and random diffusion.

In 2001, greatly inspired by the theory of island biogeography and the dispersal limitation concept (see Table ), Hubbell published an influential monograph titled “The Unified Neutral Theory of Biodiversity and Biogeography” [however, the debate dates back to 1979 ()]. Unlike the niche theory and the approach adopted by Lotka and Volterra, the neutral theory (NT) aims to model only species on the same trophic level (monotrophic communities, see Table ), species that therefore compete with each other because they all feed on the same pool of limited resources. For instance, competition arises among plant species in a forest because all of them place demands on similar resources like carbon, light, or nitrate. Other examples include species of corals, bees, hoverflies, butterflies, birds, and so on. The NT is an ecological theory within which organisms of a community have identical per-capita probabilities of giving birth, dying, migrating, and speciating, regardless of the species they belong to. Thus, from an ecological point of view, the originality of Hubbell’s NT lies in the combination of several factors: (i) it assumes competitive equivalence among interacting species; (ii) it is an individual-based stochastic theory founded on mechanistic assumptions about the processes controlling the origin and interaction of biological populations at the individual level (i.e., speciation, birth, death, and migration); (iii) it can be formulated as a dispersal limited sampling theory; and (iv) it is able to describe several macroecological patterns through just a few fundamental ecological processes, such as birth, death, and migration (; ; ; ). Although the theory has been highly criticized by many ecologists as being unrealistic (; ; ; ; ; ), it does provide very good results when describing observed ecological patterns and it is simple enough to allow analytical treatment (, , ; ; ). However, such precision does not necessarily imply that communities are truly neutral and indeed non-neutral models can also produce similar patterns (; ). Yet the NT does call into question approaches that are either more complex or equally unrealistic (; ). Moreover, NT is not only a useful tool to reveal universal patterns but also it is a framework that provides valuable information when it fails. One of the strengths of NT is that one can, in fact, falsify one or more of its assumptions and thereby actually test the theory. Few models in community ecology meet this gold standard. These features have made NT an important approach in the study of biodiversity (; ; ; ; ; ; ; ).

From a physicist’s perspective, NT is appealing as it represents a “thermodynamic” theory of ecosystems. Similar to the kinetic theory of ideal gases in physics, NT is a basic theory that provides the essential ingredients to further explore theories that involve more complex assumptions. Indeed, NT captures the fundamental approach of physicists, which can be summarized by Einstein’s celebrated quote “Make everything as simple as possible, but not simpler.” Finally, it should be noted that the NT of biodiversity is basically the analog of the theory of neutral evolution in population genetics () and indeed several results obtained in population genetics can be mapped to the corresponding ecological case ().

Statistical physics is contributing decisively to our understanding of biological and ecological systems by providing powerful theoretical tools and innovative steps (, ) to understand empirical data about emerging patterns of biodiversity. The aim of this review is not to present a complete and exhaustive summary of all the contributions to this field in recent years—a goal that would be almost impossible in such an active and broad interdisciplinary field—but rather, we want to introduce this exciting new field to physicists that have no background in ecology and yet are interested in learning about NT. Thus, we focus on what has already been done and what issues must be addressed most urgently in this nascent field, that of statistical physics applied to ecological systems. A nice feature of this field is the availability of ecological data that can be used to falsify models and highlight their limitations. At the same time, we see how the development of a quantitative theoretical framework will enable one to better understand the multiplicity of empirical experiments and ecological data.

This review is organized into six main sections. Section is a review of several important results that have been obtained by solving neutral models at stationarity. In particular, we present the theoretical framework based on Markovian assumptions to model ecological communities, where different models may be seen as the results of different NT ensembles. We also show how NT, despite its simplicity, can describe patterns observed in real ecosystems. In Sec. , we present more recent results on dynamic quantities related to NT. In particular, we discuss the continuum limit approximation of the discrete Markovian framework, paying special attention to boundary conditions, a subtle aspect of the time-dependent solution of the NT. In Sec. , we provide examples of how space plays an essential role in shaping the organization of an ecosystem. We discuss both phenomenological and spatially implicit and explicit NT models. A final subsection is devoted to the modeling of environmental fragmentation and habitat loss. In Secs. and we propose some emerging topics in this fledgling field and present the problems currently being faced. Finally, we close the review with a section dedicated to conclusions.

II. Neutral Theory at Stationarity

Neutral theory deals with ecological communities within a single trophic level, i.e., communities whose species compete for the same pool of resources (see Table ). This means that neutral models will generally be tested on data describing species that occupy the same position in the food chain, like trees in a forest, breeding birds in a given region, butterflies in a landscape, plankton, etc. Therefore, ecological food webs with predator-prey-type interactions are not suitable to be studied with standard neutral models.

As explained in the Introduction, ecologists have been studying an array of biodiversity descriptors over the last 60 years (see Table ), including RSA, SAR, and spatial PCF (see Fig. ). For instance, the RSA represents one of the most commonly used static measures to summarize information on ecosystem diversity. The analysis of this pattern reveals that the RSA distributions in tropical forests share similar shapes, regardless of the type of ecosystem, geographical location, or the details of species interactions (see Fig. ). Therefore, the functional form of the RSA [see seminal papers by and for a theoretical explanation of its origins] has been one of the great problems studied by ecologists. Indeed, much attention has been devoted to the precise functional forms of these patterns. A meta-analysis revealed that RSA distributions have basically three shapes: more often unimodal (log-normal like) in fully sampled communities, and either power law or without any mode (log-series like) within incompletely censused regions (). Also, it is difficult to tease apart, from the RSA alone, the nature of the basic processes driving communities (), because neutral and non-neutral mechanisms coexist in nature.

FIG. 1.

Visual scheme of two important macroecological patterns (see Table ): RSA and SAR. The functional shape of the RSA depends on the spatial scale considered, while the SAR generally displays a triphasic behavior (see Sec. ). There is a growing appreciation that the various descriptors of biodiversity are intrinsically interrelated, and substantial efforts have been devoted to understand the links between them. From .

FIG. 2.

Tree relative species-abundance data from the Barro Colorado Island (BCI), Yasuni, Pasoh, Lambir, Korup, and Sinharaja plots for trees that are 10 cm in stem diameter at breast height. The frequency distributions are plotted using Preston’s binning method as described by and the bars are the observed number of species binned into $l o g (2)$ abundance categories. The first histogram bar represents $⟨ 𝜙 1 ⟩ / 2$ , the second bar $⟨ 𝜙 1 ⟩ / 2 + ⟨ 𝜙 2 ⟩ / 2$ , the third bar $⟨ 𝜙 2 ⟩ / 2 + ⟨ 𝜙 3 ⟩ + ⟨ 𝜙 4 ⟩ / 2$ , the fourth bar $⟨ 𝜙 4 ⟩ / 2 + ⟨ 𝜙 5 ⟩ + ⟨ 𝜙 6 ⟩ + ⟨ 𝜙 7 ⟩ + ⟨ 𝜙 8 ⟩ / 2$ , and so on. Here $⟨ 𝜙 𝑛 ⟩$ is the number of species with an abundance $𝑛$ . As examples, we show the fits of a density-dependent symmetric model (black line), which will be studied in greater detail in Sec. . From .

It is instructive to derive some of these functional forms by starting with a very simple but extreme neutral model, which assumes that species are independent and randomly distributed in space. This null model tells us what we should expect when the observed macroecological patterns are driven only by randomness, with no underlying ecological mechanisms. If the density of individuals in a very large region is $𝜌$ , then the probability that a species has $𝑛$ individuals within an area $𝑎$ is well approximated by a Poisson distribution with a mean $𝜌 𝑎$ . As defined in Table , this is the RSA for the area $𝑎$ . However, the empirical data are not well described by a Poisson distribution and as we shall see later on better fits are usually given by log-series, gamma, or log-normal distributions.

The SAR curve can also be calculated within a slightly more accurate model, which still assumes that species are independent and randomly situated in space. Let us now suppose that a region with area $𝐴 0$ contains $𝑆 t o t$ species in total ( $𝛼$ diversity; see Table ), and that the species $𝑖$ has $𝑛 𝑖$ individuals in $𝐴 0$ . If we consider a smaller area $𝐴$ within the region, then the probability that an individual will not be found in such area is $1 - 𝐴 / 𝐴 0$ , while the probability that the whole species $𝑖$ is not present therein is $(1 - 𝐴 / 𝐴 0) 𝑛 𝑖 = 1 - 𝑝 𝑖$ . If we now consider the random variable $𝐼 𝑖 (𝐴)$ , which is 1 when the species $𝑖$ is found within the area $𝐴$ and 0 if not, then $⟨ 𝐼 𝑖 (𝐴) ⟩ = 𝑝 𝑖 = 1 - (1 - 𝐴 / 𝐴 0) 𝑛 𝑖$ , because $𝐼 𝑖 (𝐴)$ is a Bernoulli random variable with an expectation value $𝑝 𝑖$ . Therefore, the mean number of species in the area $𝐴$ (i.e., the SAR) is simply $𝑆 (𝐴) = \sum 𝑖 ⟨ 𝐼 𝑖 (𝐴) ⟩$ , which is

𝑆 (𝐴) = 𝑆 t o t - 𝑆 t o t \sum 𝑖 = 1 (1 - 𝐴 / 𝐴 0) 𝑛 𝑖 .

(1)

Although this model was originally studied by , we now know that it significantly overestimates species diversity at almost all spatial scales ().

Beta diversity (see Table ) can be estimated under the assumptions we have mentioned. Now, regardless of the spatial distance between two individuals, the probability that two of them belong to the same species $𝑖$ is $𝑛 𝑖 (𝑛 𝑖 - 1) / 𝑁 (𝑁 - 1)$ , where $𝑁$ is the total number of individuals in the community, i.e., $𝑁 = \sum 𝑖 𝑛 𝑖$ . Therefore, the probability to find any pair of conspecific individuals is $(1 / 𝑆 t o t) \sum 𝑆 t o t 𝑖 = 1 𝑛 𝑖 (𝑛 𝑖 - 1) / 𝑁 (𝑁 - 1)$ . This means that the random placement model with independent species predicts that beta diversity should not depend on the distance between two individuals. Again, we now have clear evidence that the probability that two individuals at distance $𝑟$ belong to the same species is a decaying function of $𝑟$ ().

The failure of the random placement model to capture the RSA, SAR, and beta diversity is a clear indication that ecological patterns are driven by nontrivial mechanisms that need to be appropriately identified. Thus, we shall assess to what extent the NT at stationarity can provide predictions in agreement with empirical data.

There are two related but distinct analytical frameworks that have been used to mathematically formulate the NT of biodiversity at stationarity for both local and metacommunities: the first fixes the total population, whereas the second fixes the average total population of a community. From an ecological point of view, a local community is defined as a group of potentially interacting species sharing the same environment and resources. Mathematically, when modeling a local community the total community population abundance remains fixed. Alternatively, a metacommunity can be considered a set of interacting communities that are linked by dispersal and migration phenomena. In this case, it is the average total abundance of the whole metacommunity that is held constant. From the physical point of view, these roughly correspond to the microcanonical (fixed total abundance) and the grand-canonical (fixed average total abundance) ensembles, respectively. The microcanonical ensemble or so-called zero-sum dynamics when death and birth events always occur as a pair originates from the sampling frameworks in population genetics pioneered by Warren Ewens and Ronald Fisher (). Note that even though a fixed size sample is one way to analyze available data, for the majority of cases (apart from very small size samples) the grand-canonical ensemble approach is that used routinely in statistical physics and it provides a precise yet largely simplified description of the system. The ultimate reason for this lies in the surprising accuracy of the asymptotic expansion of the gamma function (the mathematical framework heavily uses combinatorials and factorials, etc.). The Stirling approximation can be used for very large values of the gamma function; nevertheless, it is quite accurate even for values of the arguments of the order of 20. The advantages of the master equation (ME) (see Sec. ) and the grand-canonical ensemble approach stem from their computational simplicity, which make the results more intuitively transparent.

We now introduce the mathematical tools of stochastic processes that will be used extensively in the rest of the article.

A. Markovian modeling of neutral ecological communities

1. The master equation of birth and death

Let $𝑐$ be a configuration of an ecosystem that could be as detailed as the characteristics of all individuals in the ecosystem, including their spatial locations or as minimal as the abundance of a specified subset of species. Let $𝑃 (𝑐, 𝑡 | 𝑐 0, 𝑡 0)$ be the probability that a configuration $𝑐$ is seen at time $𝑡$ , given that the configuration at time $𝑡 0$ was $𝑐 0$ [referred to as $𝑃 (𝑐, 𝑡)$ for simplicity]. For our applications, the configurations $𝑐$ are typically species abundance (denoted by $𝑛$ ).

Assuming that the stochastic dynamics are Markovian, the time evolution of $𝑃 (𝑐, 𝑡)$ is given by the ME (; )

\partial 𝑃 (𝑐, 𝑡) \partial 𝑡 = \sum 𝑐' {𝑇 [𝑐 | 𝑐'] 𝑃 (𝑐', 𝑡) - 𝑇 [𝑐' | 𝑐] 𝑃 (𝑐, 𝑡)},

(2)

where $𝑇 [𝑐 | 𝑐']$ is the transition rate from the configuration $𝑐'$ to configuration $𝑐$ . Under suitable and plausible conditions (), $𝑃 (𝑐, 𝑡)$ approaches a stationary solution at long times $𝑃 𝑠 (𝑐)$ , which satisfies

\sum 𝑐' {𝑇 [𝑐 | 𝑐'] 𝑃 𝑠 (𝑐') - 𝑇 [𝑐' | 𝑐] 𝑃 𝑠 (𝑐)} = 0 .

(3)

Equation is typically intractable with analytical tools, because it involves a sum of all the configurations. If each term in the summation is zero, i.e.,

𝑇 [𝑐 | 𝑐'] 𝑃 𝑠 (𝑐') - 𝑇 [𝑐' | 𝑐] 𝑃 𝑠 (𝑐) = 0,

(4)

detailed balance is said to hold. A necessary and sufficient condition for the validity of detailed balance is that for all possible cycles in the configuration space, the probability of walking through it in one direction is equal to the probability of walking through it in the opposite direction. Given a cycle ${𝑐 1, 𝑐 2, \dots, 𝑐 𝑛 - 1, 𝑐 1}$ , detailed balance holds if and only if, for every such cycle,

𝑇 [𝑐 1 | 𝑐 2] 𝑇 [𝑐 2 | 𝑐 3] \dots 𝑇 [𝑐 𝑛 - 2 | 𝑐 𝑛 - 1] 𝑇 [𝑐 𝑛 - 1 | 𝑐 1] = 𝑇 [𝑐 1 | 𝑐 𝑛 - 1] 𝑇 [𝑐 𝑛 - 1 | 𝑐 𝑛 - 2] \dots 𝑇 [𝑐 3 | 𝑐 2] 𝑇 [𝑐 2 | 𝑐 1] .

(5)

This condition evidently corresponds to a time-reversible condition.

Now let us apply these mathematical tools to the study of community dynamics that is driven by random demographical drift. Consider a well-mixed local community. This is equivalent to saying that the distribution of species in space is not relevant, which should hold for an ecosystem with a linear size smaller than or of the same order as the seed dispersal range. In this case one can use $𝑐 = (𝑛 1, \dots, 𝑛 𝑆) = 𝐧$ , where $𝑛 𝑖$ is the population of the $𝑖$ th species and the system contains $𝑆$ species. We can rewrite Eq. in the following way:

\partial 𝑃 (𝐧, 𝑡) \partial 𝑡 = \sum 𝐧' \neq 𝐧 𝑇 [𝐧 | 𝐧'] 𝑃 (𝐧', 𝑡) - \sum 𝐧' \neq 𝐧 𝑇 [𝐧' | 𝐧] 𝑃 (𝐧, 𝑡) .

(6)

In this case, $𝑇 [𝐧' | 𝐧]$ can take into account birth and death, as well as immigration from a metacommunity. The simplest hypothesis is that $𝑇 [𝐧' | 𝐧]$ is the result of $𝑆$ elementary birth and death processes that occur independently for each of the $𝑆$ species, i.e.,

𝑇 [𝐧' | 𝐧] = 𝑆 \sum 𝑘 = 1 \prod 𝑖 \neq 𝑘 𝛿 𝑛' 𝑖, 𝑛 𝑘 [𝛿 𝑛' 𝑘, 𝑛 𝑘 + 1 𝑇 𝑘 (𝑛 + 1 | 𝑛) + 𝛿 𝑛' 𝑘, 𝑛 𝑘 - 1 𝑇 𝑘 (𝑛 - 1 | 𝑛)],

(7)

where $𝛿 𝑘', 𝑘$ is a Kronecker delta,

𝑇 𝑘 (𝑛 + 1 | 𝑛) = 𝑏 (𝑛, 𝑘)

(8)

is the birth rate, and

𝑇 𝑘 (𝑛 - 1 | 𝑛) = 𝑑 (𝑛, 𝑘)

(9)

the death rate. This particular choice corresponds to a sort of mean-field (ME) approach (; , ; ; ; ; ).

Our many-body ecological problem can also be formulated in a language more familiar to statistical physicists, where we consider the distribution of balls into urns. The “urns” are the species and the “balls” are the individuals. Birth and death processes correspond to adding or removing a ball to or from one of the urns using a rule as dictated by Eqs. and . Equation with the choice can be simplified if one assumes that the initial condition is factorized as $𝑃 (𝐧, 𝑡 = 0) = \prod 𝑆 𝑘 = 1 𝑃 𝑘 (𝑛 𝑘, 𝑡 = 0)$ . In this case the solution is again factorized as $𝑃 (𝐧, 𝑡) = \prod 𝑆 𝑘 = 1 𝑃 𝑘 (𝑛 𝑘, 𝑡)$ where each $𝑃 𝑘$ satisfies the following ME in 1 degree of freedom $𝑛 𝑘$ :

\partial 𝑃 𝑘 (𝑛 𝑘, 𝑡) \partial 𝑡 = 𝑃 𝑘 (𝑛 𝑘 - 1, 𝑡) 𝑏 (𝑛 𝑘 - 1, 𝑘) + 𝑃 𝑘 (𝑛 𝑘 + 1, 𝑡) 𝑑 (𝑛 𝑘 + 1, 𝑘) - [𝑏 (𝑛 𝑘, 𝑘) + 𝑑 (𝑛 𝑘, 𝑘)] 𝑃 𝑘 (𝑛 𝑘, 𝑡) .

(10)

The stationary solution of Eq. is easily seen to satisfy the detailed balance (; ).

First we note that because of the neutrality hypothesis species are assumed to be demographically identical and therefore, we can drop the $𝑘$ dependence factor from Eqs. and . In other words, we can concentrate on the probability $𝑃 (𝑛)$ that a given species (urn) has an abundance $𝑛$ (balls). In this case, species do not interact and thus in our calculation we can follow a particular species (the urns are taken to be independent). Following Eq. , it is easy to see that the solution of the birth-death ME that satisfies the detailed balance condition and that thus corresponds to equilibrium is (; )

𝑃 𝑠 (𝑛) = 𝑃 0 𝑛 \prod 𝑧 = 1 𝑏 (𝑧 - 1) 𝑑 (𝑧),

(11)

where $𝑃 0$ can be calculated by the normalization condition $\sum 𝑛 𝑃 𝑠 (𝑛) = 1$ , and assuming that all the rates are positive. In particular, $𝑏 (0) > 0$ and $𝑑 (0) = 0$ .

When there are $𝑆$ urns, all satisfying the same birth-death rules, the general and unique equilibrium solution is

𝑃 𝑠 (𝑛 1, 𝑛 2, \dots, 𝑛 𝑆) = 𝑆 \prod 𝑘 = 1 𝑃 𝑠 (𝑛 𝑘) .

(12)

Depending on the functional form of $𝑏 (𝑛)$ and $𝑑 (𝑛)$ , one can readily work out the desired $𝑃 𝑠 (𝐧)$ . We start with some cases that are familiar to physicists () and then move onto more ecologically meaningful cases (, ).

2. Physics ensembles

The random walk and Bose-Einstein distribution.: If one chooses $𝑏 (𝑛) = 𝑏 0$ and $𝑑 (𝑛) = 𝑑 1$ for $𝑛 \geq 1$ , and $𝑑 (𝑛) = 0$ otherwise [or alternatively, $𝑏 (𝑛) = (𝑛 + 1) 𝑏 0$ and $𝑑 (𝑛) = 𝑑 1 𝑛$ ], one obtains a pure exponential distribution $𝑃 (𝑛) = 𝑟 𝑛 (1 - 𝑟)$ where $𝑟 = 𝑏 0 / 𝑑 1$ . Note that in both these cases, $𝑏 (𝑛 - 1) / 𝑑 (𝑛)$ is the same. Substituting this into Eq. gives us the well-known Bose-Einstein distribution for nondegenerate energy levels, i.e.,

𝑃 (𝑛 1, 𝑛 2, \dots, 𝑛 𝑆) = 𝑟 𝑁 (1 - 𝑟) 𝑆,

(13)

where $𝑁 = \sum 𝑘 𝑛 𝑘$ . Here $𝑟$ corresponds to $𝑒 - 𝛽 (𝜖 - 𝜇)$ in the grand-canonical ensemble.

The Fermi-Dirac distribution: If $𝑑 (𝑛) = 𝑑 1 𝑛$ and $𝑏 (𝑛) = 0$ for any $𝑛$ other than 0, and equal to $𝑏 0$ otherwise, then $𝑃 (𝑛) = 𝑟 𝑛 / (1 + 𝑟)$ for $𝑛 = 0$ or $𝑛 = 1$ and $𝑃 (𝑛) = 0$ for other values of $𝑛$ , and accordingly the Fermi-Dirac distribution is achieved

𝑃 (𝑛 1, 𝑛 2, \dots, 𝑛 𝑆) = 𝑟 \sum 𝑘 𝑛 𝑘 (1 + 𝑟) - 𝑆 for 𝑛 𝑘 = 0 or 𝑛 𝑘 = 1 \forall 𝑖,

(14)

𝑃 (𝑛 1, 𝑛 2, \dots, 𝑛 𝑆) = 0 otherwise .

(15)

Boltzmann counting: If $𝑑 (𝑛) = 𝑑 1 𝑛$ and $𝑏 𝑛 = 𝑏 0$ for all $𝑛$ , then one obtains a Poisson distribution $𝑃 (𝑛) = 𝑒 - 𝑟 𝑟 𝑛 / 𝑛!$ , and this leads to

𝑃 (𝑛 1, 𝑛 2, \dots, 𝑛 𝑆) \propto 𝑟 𝑁 \prod 𝑆 𝑘 = 1 𝑛 𝑘! .

(16)

This is the familiar grand-canonical ensemble Boltzmann counting in physics, where $𝑟$ plays the role of fugacity. It is noteworthy that, unlike the conventional classical treatment (), where an additional factor of $𝑁!$ is obtained, here one gets the correct Boltzmann counting and thereby avoids the well-known Gibbs paradox (). Thus, if one were to ascribe energy values to each of the boxes and enforce a fixed average total energy, the standard Boltzmann result would be obtained whereby the probability of occupancy of an energy level $𝜖$ is proportional to $𝑒 - 𝛽 𝜖$ , where $𝛽$ is proportional to the inverse of the temperature.

3. Ecological ensembles

Density independent dynamics: We now consider the dynamic rules of birth, death, and speciation that govern the population of an individual species. The most simple ecologically meaningful case is to consider $𝑑 (𝑛) = 𝑑 𝑛$ and $𝑏 (𝑛) = 𝑏 𝑛$ for $𝑛 > 0$ , and $𝑏 (0) > 0$ , $𝑑 (0) = 𝑏 (- 1) = 0$ . We define $𝑟 = 𝑏 / 𝑑$ . Moreover, in order to ensure that the community will not become extinct at longer times, speciation may be introduced by ascribing a nonzero probability of the appearance of an individual from a new species, i.e., $𝑏 (0) = 𝑏 0 = 𝜈$ . In this case the probability for a species of having $𝑛$ individuals at stationarity is

𝑃 (𝑛) = ˜ 𝜈 [1 - ˜ 𝜈 l n (1 - 𝑟)] - 1 𝑟 𝑛 / 𝑛,

(17)

where $˜ 𝜈 = 𝜈 / 𝑏$ and $𝑛 > 0$ and $𝑃 (0) = 1 / [1 - ˜ 𝜈 l n (1 - 𝑟)]$ .

The RSA $⟨ 𝜙 (𝑛) ⟩$ is the average number of species with a population $𝑛$ , and this is simply ()

⟨ 𝜙 (𝑛) ⟩ = 𝑆 𝑃 (𝑛) = 𝜃 𝑟 𝑛 / 𝑛 .

(18)

This is the celebrated Fisher log-series distribution, i.e., the distribution Fisher proposed as that describing the empirical RSA in real ecosystems (). The parameter $𝜃 = 𝑆 ˜ 𝜈 / [1 - ˜ 𝜈 l n (1 - 𝑟)]$ is known as the Fisher number or biodiversity parameter.

In 1948, F. W. Preston, published a paper () challenging Fisher’s point of view. He showed that the log series is not a good description for the data from a large sample of birds. In fact, he observed an internal mode in the RSA that was absent in a log-series distribution. In particular, Preston introduced a way to plot the experimental RSA data by octave abundance classes (i.e., [ $2 𝑘$ , $2 𝑘 + 1$ ], for $𝑘 = 1, 2, 3, \dots$ ), showing that a good fit of the data was represented by a log-normal distribution. Indeed, there are several examples of RSA data that display this internal mode feature (see Fig. ). The intuition of Preston was that the shape of the RSA must depend on the sampling intensity or size of the community. Conventionally, when studying ecological communities, ecologists separate them into two distinct classes: small local communities (e.g., on an island) and metacommunities of much larger communities or those composed of several smaller local communities (see Fig. ). The neutral modeling schemes for these two cases, that we denote by the subscripts $𝐿$ and $𝑀$ , respectively, are not the same as the ecological processes involved differ. In fact, the immigration rate ( $𝑚$ ) in a local community is a crucial parameter as the community is mainly structured by dispersal limited mechanisms, and the speciation rate ( $𝜈$ ) can be neglected. On the other hand, in a metacommunity immigration does not occur (species colonize within the community) and the community is shaped by birth-death-competition processes, although speciation $𝜈$ also plays an important role. A metacommunity can also be thought of as consisting of many small semi-isolated local communities, each of which receives immigrants from other local communities. When considering metacommunity dynamics, the natural choice is to put a soft constraint on $𝐽 𝑀$ , i.e., the total number of individuals is free to fluctuate around the average population $⟨ 𝐽 𝑀 ⟩$ . Indeed, one finds that at the largest $𝐽 𝑀$ limit, the results obtained with a hard constraint on $𝐽 𝑀$ are equivalent to those with a soft constraint. On the other hand, when considering a local community, it is safer to place a hard constraint on the total population $𝐽 𝐿$ . In both these cases, $𝑆$ represents the total number of species that may potentially be present in the community, while the average number of species observed in the community is denoted by $⟨ 𝑆 ⟩$ .

FIG. 3.

Different neutral models of community ecology. In all these models, $𝑛 𝑖$ represents the abundance of the species $𝑖$ in the community, and $𝐽 𝐿$ and $𝐽 𝑀$ the total abundance in the local and metacommunities, respectively. (a) Hubbell’s zero-sum neutral model (). In the local community, each death is immediately followed by a birth or an immigration event. Speciation (or, equivalently, immigration) enables diversity to be maintained in the metacommunity. (b) Local community with immigrants from a metacommunity (). The local community now interacts with the metacommunity through a migration process ( $𝑚$ ). (c) Coalescent-type approach (), where community members are traced back to the ancestors that once immigrated into the community. (d) Joint RSA of many local communities (). The whole metacommunity RSA distribution is built by considering the joint RSA distributions of multiple local communities.

There are several ecological meaningful mechanisms that can generate a bell-shaped Preston-like RSA. The first of these involves density-dependent effects on birth and death rates. The second involves considering a Fisher log series as the RSA of a metacommunity acting as a source of immigrants to a local community embedded within it. The dynamics of the local community are governed by births, deaths, and immigration, whereas the metacommunity is characterized by births, deaths, and speciation. This leads to a local community RSA with an internal mode (; ; , ). A third way is to incrementally aggregate several local communities (see Appendix ).

Local dynamics with density-dependent rates: One of the most fundamental and long debated questions in community ecology is to understand what maintains species diversity within communities (; ). Mechanisms which are able to promote diversity include competitive trade-offs among species because of species-specific traits, balance between speciation and extinction, frequency or density dependence, and environmental variability. For instance, processes that hold the abundance of a common species in check inevitably lead to rare-species advantages, given that the space or resources freed up by density-dependent death can be exploited by less-common species. Therefore, interspecies frequency dependence is the community-level consequence of intraspecies density dependence, and, thus, they may be thought of as two different manifestations of the same phenomenon (). Because these processes are capable of driving community-level patterns such as SAR or RSA, one might hope to identify mechanisms by delving into observed patterns. However, communities with different governing processes can unexpectedly show similar patterns, thus suggesting that their shape cannot be, in general, used to identify specific mechanisms (; ).

In this review we focus on two of the most prominent hypotheses that explain species coexistence through frequency and density dependence: the Janzen-Connell (; ) and the Chesson-Warner hypotheses (). These mechanisms generally predict the reproductive advantage of a rare species due to ecological factors and they can be readily captured in a common mathematical framework that is presented later.

The Janzen-Connell hypothesis postulates that host-specific pathogens or predators act in the vicinity of the maternal parent. Thus, seeds that disperse farther away from the mother are more likely to escape mortality. This spatially structured mortality effect suppresses the uncontrolled population growth of locally abundant species relative to uncommon species, thereby producing a reproductive advantage to a rare species. The Chesson-Warner storage hypothesis explores the consequences of a variable external environment and it relies on three empirically validated observations: species respond in a species-specific manner to the fluctuating environment, there is a covariance between the environment and intraspecies and interspecies competition, or life history stages buffer the growth of population against unfavorable conditions. Such conditions prevail when species have similar per-capita rates of mortality but they reproduce asynchronously and there are overlapping generations.

We begin by noting that the mean number of species with $𝑛$ individuals $⟨ 𝜙 𝑛 ⟩$ is not determined by the absolute rates of birth or death but rather by their ratio $𝑟 𝑖, 𝑘 = 𝑏 (𝑖, 𝑘) / 𝑑 (𝑖 + 1, 𝑘)$ . This follows from the observation that $⟨ 𝜙 𝑛 ⟩$ is proportional to $⟨ 𝑟 1, 𝑘 𝑟 2, 𝑘 \dots 𝑟 𝑛 - 1, 𝑘 ⟩ 𝑘$ , where the average $⟨ \dots ⟩ 𝑘$ is obtained from all the species. This simple formulation () is sufficiently general to represent the communities of either symmetric species (in which all the species have the same demographic birth and death rates) or the case of asymmetric or distinct species. The more general asymmetric situation captures niche differences and/or differing immigration fluxes that might arise from the different relative abundances of distinct species in the surrounding metacommunity (see Table ).

The density dependence arising from the Janzen-Connell effect is spatially explicit, in plant communities at least, because it is caused by interactions among neighboring individuals. Here instead we formulate a density dependence which is spatially implicit and only abundance driven. We therefore introduce a modified symmetric theory that incorporates a rare-species advantage or common species disadvantage by making $𝑟 𝑛$ a decreasing function of abundance. The equations of density dependence in the per-capita birth and death rates for an arbitrary species of abundance $𝑛$ are

𝑏 (𝑛) 𝑛 = 𝑏 [1 + 𝑏 1 𝑛 + 𝑜 (1 𝑛 2)]

and

𝑑 (𝑛) 𝑛 = 𝑑 [1 + 𝑑 1 𝑛 + 𝑜 (1 𝑛 2)],

for $𝑛 > 0$ as the leading term of a power series in ( $1 / 𝑛$ ),

𝑏 (𝑛) 𝑛 = 𝑏 \infty \sum 𝑙 = 0 𝑏 𝑙 𝑛 - 𝑙

and

𝑑 (𝑛) 𝑛 = 𝑑 \infty \sum 𝑙 = 0 𝑑 𝑙 𝑛 - 𝑙,

where $𝑏 𝑙$ and $𝑑 𝑙$ are constants. This expansion captures the essence of density dependence by ensuring that the per-capita rates decrease and approach a constant value for a large $𝑛$ , given that the higher order terms are negligible. As noted earlier, the quantity that controls the RSA distribution is the ratio $𝑏 𝑛 / 𝑑 𝑛 + 1$ . Thus, the birth and death rates $𝑏 𝑛$ and $𝑑 𝑛$ can be defined up to $𝑓 (𝑛 + 1)$ and $𝑓 (𝑛)$ , respectively, where $𝑓$ is any arbitrary well-behaved function.

Strikingly, any relative abundance data can be considered as arising from effective density-dependent processes in which the birth and death rates are given by these expressions. Thus, one would expect that the per-capita birth rate or fecundity drops as the abundance increases, whereas mortality ought to increase with abundance. Indeed, the per-capita death rate can be arranged to be an increasing function of $𝑛$ , as observed in nature, by choosing an appropriate function $𝑓$ and appropriately adjusting the birth rate so that the ratio $𝑏 𝑛 / 𝑑 𝑛 + 1$ remains the same. This then ensures that the RSA does not vary.

The mathematical formulation of the density dependence may seem unusual to ecologists familiar with the logistic or Lotka-Volterra systems of equations, wherein the density dependence is typically described as a polynomial expansion of powers of $𝑛$ truncated at the quadratic level. Such an expansion is valid when the characteristic scale of $𝑛$ is determined by a fixed carrying capacity. Conversely, here the range of $𝑛$ is from 1 to an arbitrarily large value and not to some carrying capacity. Therefore, an expansion in terms of the powers of $1 / 𝑛$ is more appropriate. For this symmetric model (), and bearing in mind that

⟨ 𝜙 𝑛 ⟩ = 𝑆 𝑃 0 𝑛 - 1 \prod 𝑖 = 0 𝑏 (𝑖) 𝑑 (𝑖 + 1),

one readily arrives at the following relative species-abundance relationship:

⟨ 𝜙 𝑛 ⟩ = 𝜃 𝑥 𝑛 𝑛 + 𝑐,

(19)

where $𝑥 = 𝑏 / 𝑑$ and, for parsimony, we make the simple assumption that $𝑏 1 = 𝑑 1 = 𝑐$ and the higher order coefficients, $𝑏 2$ , $𝑑 2$ , $𝑏 3$ , $𝑑 3$ , $\dots$ , are all 0. The biodiversity parameter $𝜃$ is the normalization constant that ensures the total number of species in the community is $𝑆$ , and it is given by $𝜃 = 𝑆 [(1 + 𝑐) / 𝑐 𝑥] / 𝐹 (1 + 𝑐, 2 + 𝑐, 𝑥)$ , where $𝐹 (1 + 𝑐, 2 + 𝑐, 𝑥)$ is the standard hypergeometric function. The parameter $𝑐$ measures the strength of the symmetric density and frequency dependence in the community, and it controls the shape of the RSA distribution. This simple model () does a good job of matching the patterns of abundance distribution observed in the tropical forest communities throughout the world (see Fig. ). Note that when $𝑐 \to 0$ (the case of no density dependence), one obtains the Fisher log series. In this case $𝜃$ captures the effects of speciation.

Various analytical solutions for the Markovian (ME) approach have been suggested in the literature. Based on pioneering work by , an analytical solution was proposed for the ME when the size of both the local community and the metacommunity was fixed (; ; ). Here we refer to the work of , in whose model the birth and death rates in the metacommunity are given by

𝑇 (𝑛 + 1 | 𝑛) = 𝑏 𝑀 (𝑛) = 𝐽 𝑀 - 𝑛 𝐽 𝑀 𝑛 𝐽 𝑀 - 1 (1 - 𝜈),

(20)

𝑇 (𝑛 - 1 | 𝑛) = 𝑑 𝑀 (𝑛) = 𝑛 [𝐽 𝑀 - 𝑛 + (𝑛 - 1) 𝜈] 𝐽 𝑀 (𝐽 𝑀 - 1) .

(21)

During a unit of time, the population ( $𝑛 \neq 0, 𝐽 𝑀$ ) of the species under consideration can increase by one individual if we uniformly pick at random an individual that does not belong to that species [this occurs with probability $(𝐽 𝑀 - 𝑛) / 𝐽 𝑀$ ], we substitute it with an individual randomly picked from the considered species [with probability $𝑛 / (𝐽 𝑀 - 1)$ ], and, finally, there is no speciation (probability $1 - 𝜈$ ). Multiplying all these terms, we get the birth rate in Eq. . On the other hand, in the same unit of time the population of the species under consideration can decrease by one individual if we uniformly pick at random an individual that belongs to that species (probability $𝑛 / 𝐽 𝑀$ ) and we substitute it either with a completely new individual (speciation with probability $𝜈$ ) or with an individual which does not belong to the considered species [with probability $(𝐽 𝑀 - 𝑛) / (𝐽 𝑀 - 1)$ ]. These terms provide the death rate in Eq. .

Let $⟨ 𝜙 𝑛 (𝑡) ⟩ 𝑀$ designate the average number of species of abundance $𝑛$ in the metacommunity. Then in the stationary limit ( $𝑡 \to \infty$ ), the species that contribute to $⟨ 𝜙 𝑛 (𝑡) ⟩ 𝑀$ are those which entered the system by mutation at time $𝑡 - 𝜏$ and that have reached size $𝑛$ at time $𝑡$ (; ). As speciation is a Poissonian event (of rate $𝜈$ ) and due to neutrality, all species have the same probability $𝑝 (𝜏) 𝑑 𝜏 = 𝜈 𝑑 𝜏$ of appearing between the time interval $[𝜏, 𝜏 + 𝑑 𝜏]$ . Thus, the time evolution of $⟨ 𝜙 𝑛 (𝑡) ⟩$ is then simply given by $⟨ 𝜙 𝑛 (𝑡) ⟩ = \int 𝑡 0 𝑃 (𝑛, 𝑡 - 𝜏) 𝑝 (𝜏) 𝑑 𝜏 = 𝜈 \int 𝑡 0 𝑃 (𝑛, 𝜏) 𝑑 𝜏$ , where $𝑃 (𝑛, 𝑡)$ is the solution of the ME with transition rates given by Eqs. and and with the initial condition $𝑃 (𝑛, 0) = 𝛿 𝑛, 1$ ( $𝛿$ represents the Kronecker delta). Using Laplace transforms, the stationary RSA for the metacommunity of size $𝐽 𝑀$ is ()

⟨ 𝜙 𝑘 ⟩ 𝑀 = 𝜃 Γ (𝐽 𝑀 + 1) Γ (𝐽 𝑀 + 𝜃 - 𝑘) 𝑘 Γ (𝐽 𝑀 + 1 - 𝑘) Γ (𝐽 𝑀 + 𝜃),

(22)

where $𝜃 = (𝐽 𝑀 - 1) 𝜈 / (1 - 𝜈)$ .

Therefore the density of species of relative abundance $𝜔 = 𝑘 / 𝐽 𝑀$ in the metacommunity can be written as

𝑔 𝑀 (𝜔) = 𝜃 (1 - 𝜔) 𝜃 - 1 𝜔,

(23)

where $𝑔 𝑀 (𝜔) = l i m 𝐽 𝑀 \to \infty 𝐽 𝑀 ⟨ 𝜙 𝑘 ⟩ 𝑀$ and with $𝜔$ and $𝜃$ kept fixed (implying that $𝜈$ tends to zero and $𝑘$ tends to infinity).

Let us now consider a local community (of size $𝐽 𝐿$ ) in contact with the aforementioned metacommunity. As $𝐽 𝑀 ≫ 𝐽 𝐿$ , then mutations within the local community can be neglected ( $𝜈 \to 0$ ) and the equilibrium distribution in the metacommunity is negligibly modified by migration from or toward the community. Nevertheless, there is migration from the metacommunity to the local community and, thus, a migration rate $𝑚$ must be included in the transition rates for the local community dynamics:

𝑏 𝐿 (𝑛) = 𝐽 𝐿 - 𝑛 𝐽 𝐿 [𝑛 𝐽 𝐿 - 1 (1 - 𝑚) + 𝑚 𝜔],

(24)

𝑑 𝐿 (𝑛) = 𝑛 𝐽 𝐿 [𝐽 𝐿 - 𝑛 𝐽 𝐿 - 1 (1 - 𝑚) + 𝑚 (1 - 𝜔)],

(25)

where $𝜔 = 𝑘 / 𝐽 𝑀$ is the relative abundance of the species in the metacommunity and $𝑘$ is the abundance of the tracked species within the metacommunity.

The stationary solution $𝑃 𝑛 (𝜔)$ of the corresponding ME can be calculated exactly using Eq. , considering $𝑚$ as the immigration rate ():

𝑃 𝑛 (𝜔) = (𝐽 𝐿 𝑛) (𝜇 𝜔) 𝑛 [𝜇 (1 - 𝜔)] 𝐽 𝐿 - 𝑛 (𝜇) 𝐽 𝐿,

(26)

where $(𝐽 𝐿 𝑛)$ is the binomial coefficient, $(𝑎) 𝑛 = Γ (𝑎 + 𝑛) / Γ (𝑎)$ is the Pochammer symbol, and $𝜇 = (𝐽 𝐿 - 1) 𝑚 / (1 - 𝑚)$ .

Finally, the average number of species with abundance $𝑛$ in the local community $⟨ 𝜙 𝑛 ⟩ 𝐿$ is given by

⟨ 𝜙 𝑛 ⟩ 𝐿 = 𝐽 𝑀 \sum 𝑘 = 1 𝑃 𝑛 (𝑘 / 𝐽 𝑀) ⟨ 𝜙 𝑘 ⟩ 𝑀,

(27)

which displays an internal mode for appropriate values of the model parameters. A generalization of this result to a system of two communities of arbitrary yet fixed sizes that are subject to both speciation and migration () has also been carried out.

suggested that the log-series solution for the species-abundance distribution in the metacommunity is applicable only for species-rich communities, and that it does not adequately describe species-poor metacommunities. Application of the ME approach to the asymmetric species case was considered () and it was further demonstrated that the species-abundance distribution has exactly the same sampling formula for both zero-sum and nonzero-sum models within the neutral approximation (; ). The simplest mode of speciation, a point mutation, has been the one most commonly used to derive the species-abundance distribution of the metacommunity. A Markovian approach incorporating various modes of speciation such as random fission was also presented by , . An innovative way to model speciation was used by . They introduced protracted speciation, i.e., a gradual process whereby new species are created with a few individuals, instead of being an instantaneous process starting with exactly one individual. This gradual speciation improves predictions of species lifetimes, speciation rates, and the number of rare species.

B. Coalescent approach to neutral theory

An alternative approach to compute the distribution of individuals in a community is based on coalescent theory. The idea of this approach consists of tracing each member of the community back to their ancestors that first immigrated into the community (; ). A succinct comparison between the coalescent approach and other methods is presented by and .

Within this framework, the local community consists of a fixed number of individuals $𝐽$ that undergo zero-sum dynamics. At each time step a randomly chosen individual dies and is replaced by a new individual that is an offspring of a randomly chosen individual in the community (of size $𝐽 - 1$ ), or by one of $𝐼$ potential immigrants from the metacommunity. The immigration rate is thus $𝑚 = 𝐼 / (𝐼 + 𝐽 - 1)$ , where the parameter $𝐼$ is merely a proxy for the immigration rate. The metacommunity is governed by the same dynamics, with speciation replacing immigration.

Two observations are crucial to derive the sampling formula for the species-abundance distribution. First, because there is no speciation in the local community, each individual is either an immigrant or a descendent of an immigrant. Thus, the information pertaining to the species-abundance distribution is specified by considering “the ancestral tree” (tracing back each individual to its immigrant ancestor, which is somewhat similar to the phylogenetic tree construction in genetics) and the species composition of the set of all the ancestors. Second, this set of all ancestors can be considered as a random sample from the metacomunity and its species-abundance distribution is provided by the Ewens’ sampling formula (), corresponding to the Fisher’s log series in the limit of large sample size that describes neutral population genetics models with speciation and no immigration. As we do not know the ancestor of each individual, the final formula for the species-abundance composition of the data comes from summing the probabilities of all possible combinations of ancestors that give rise to the observed data.

In a tour de force calculation, the resulting sampling formula is found to be ()

𝑃 [𝐷 | 𝜃, 𝑚, 𝐽] = 𝐽! \prod 𝑆 𝑖 = 1 𝑛 𝑖 \prod 𝐽 𝑗 = 1 Φ 𝑗! 𝜃 𝑆 (𝐼) 𝐽 𝐽 \sum 𝐴 = 𝑆 𝐾 (𝐷, 𝐴) 𝐼 𝐴 (𝜃) 𝐴,

(28)

where $𝑃 [𝐷 | 𝜃, 𝑚, 𝐽]$ is the probability of observing the empirical data $𝐷$ consisting of $𝑆$ species with abundances $𝑛 𝑖$ , $𝑖 \in [1, 𝑆]$ given the value of the fundamental biodiversity number $𝜃$ of the metacommunity, the immigration rate $𝑚$ , and the local community size $𝐽$ . $Φ 𝑗$ is the number of species with abundance $𝑗$ ,

𝐾 (𝐷, 𝐴) = \sum {𝑎 1, \dots, 𝑎 𝑆 | \sum 𝑆 𝑖 = 1 𝑎 𝑖 = 𝐴} 𝑆 \prod 𝑖 = 1 ¯ 𝑠 (𝑛 𝑖, 𝑎 𝑖) ¯ 𝑠 (𝑎 𝑖, 1) ¯ 𝑠 (𝑛 𝑖, 1),

(29)

where $¯ 𝑠$ is the unsigned Stirling number of the first kind. Finally, $(𝑥) 𝑦$ is the Pochhammer symbol $\prod 𝑦 𝑖 = 1 (𝑥 + 𝑖 - 1)$ .

These expressions are simplified versions of the original results presented and as expected they provide a comparable fit to the ME approach for the data from the tropical forests. It is heartening that equivalent results can be derived using very different approaches.

III. Dynamics in Neutral Theory

So far we have focused on the implications of neutral theory when models describing neutral patterns reach a steady state. The stationary condition allows one to take advantage of a variety of different mathematical techniques to obtain analytical expressions of ecological patterns. However, stationarity is not always a good approximation, either because the ecosystem is still in a state of flux or because the assumption may hide different and important processes that lead to the same final steady state. Furthermore, one can calculate time-dependent correlation functions (); see, for example, Eq. . It is therefore essential to understand the temporal behavior of ecosystems in order to discriminate between ecological processes that would otherwise be indistinguishable at stationarity. Statistical comparisons between time-dependent patterns are usually more difficult than those between stationary patterns because they require more data and long empirical time series that are rarely available. In addition, although time-dependent solutions facilitate stronger tests when confronted with data, they are more difficult to obtain and this is the reason why only a few studies have investigated the temporal behavior of neutral models.

An important method to study the time dependence is van Kampen’s system size expansion. Assuming that the total population of the system $𝑁$ is very large and the initial population is of the order of $𝑁$ , one expects that, at short times at least, the probability distribution of the system $𝑃 (𝑛, 𝑡)$ would peak sharply around the macroscopic value $𝑛 = 𝑁 𝜙 (𝑡)$ with a width of order $\sim \sqrt 𝑁$ . The function $𝜙 (𝑡)$ is chosen to describe the evolution of the peak. These considerations lead to the ansatz $𝑛 = 𝑁 𝜙 (𝑡) + \sqrt 𝑁 𝜉$ , where $𝜉$ is a new random variable. When expressing the ME in terms of $𝜉$ and expanding it in $\sqrt 𝑁$ , at order $\sqrt 𝑁$ one recovers the deterministic evolution of the macroscopic state of the system, and at order $𝑁 0$ one obtains a linear Fokker-Planck (FP) equation whose coefficients depend on time through the function $𝜙 (𝑡)$ . Therefore, fluctuations around the peak are Gaussian in a first approximation [details of the method can be found in ]. This expansion has some attractive features but it is usually only a good approximation of the temporal evolution of the original ME at short or intermediate temporal scales, since the Gaussian behavior is usually lost over longer periods. The time dependence of mean-field neutral models was investigated by means of van Kampen’s system size expansion in , and a good agreement with simulations was found at early times.

Some neutral models can also be formulated in terms of discrete-time Markov processes that have the form $𝑝 𝑖 (𝑡 + 1) = \sum 𝑗 𝑄 𝑖 𝑗 𝑝 𝑗 (𝑡)$ , where $𝑝 𝑖 (𝑡)$ is the probability that a species has $𝑖$ individuals at time $𝑡$ , and $𝑄 𝑖 𝑗$ is the probability that a species changes its abundance from $𝑗$ to $𝑖$ in one time step (i.e., the transition matrix). If the community has $𝑁$ individuals in total and either one birth or one death is allowed in one time step, then $𝑄$ is a tridiagonal matrix with dimension $𝑁 + 1$ . Solving these models in time is basically equivalent to finding the eigenvalues along with the left and right eigenvectors of the matrix $𝑄$ , which is usually a nontrivial task. Using a theorem derived in population genetics, was able to obtain the eigendecomposition of the transition matrix of a neutral model that describes the behavior of the local community in the original Hubbell model (). Also, he was able to show that the first two eigenvectors are sufficient to provide a good approximation of the full time-dependent solution.

Others have focused on specific temporal patterns that can sometimes be calculated more easily. For instance, a measure that has been used to quantify the diversity of a community is the Simpson index () $𝐷$ , which is defined as the probability that two individuals drawn randomly from a well-mixed community belong to the same species. If the individuals are drawn with replacement, then one obtains $𝐷 = \sum 𝑖 𝑛 2 𝑖 / 𝑁 2$ , where $𝑛 𝑖$ is the number of individuals of species $𝑖$ (we have already encountered this index in Sec. when we calculated the beta diversity in the random placement model). This index is also a good indicator of the relative importance of spatial effects in a community. The dynamic evolution of $𝐷$ is usually simple enough to allow analytical calculations (). Other interesting dynamic properties of neutral models have been studied, where species’ extinction and monodominance were investigated by looking into the first-passage properties and fixation times of the process (; ). Empirical power-law relationships between the temporal mean and variance of population fluctuations (the so-called Taylor’s power law) are also in good agreement with neutral predictions ().

Although it is not mandatory from the neutral assumption, models often assume that demographic stochasticity is the main source of fluctuations in stochastic neutral dynamics (; ). Demographic randomness originates from the intrinsic stochastic nature of birth and death events within a discrete population of individuals. However, other sources may be important, such as environmental stochasticity that, by contrast, encompasses effects of abiotic and biotic environmental variables on all individuals [strictly speaking, we should more correctly refer to symmetric, instead of neutral, models (see Table ) when these incorporate environmental fluctuations, because the per-capita vital rates vary across individuals; however, many authors do not appreciate this distinction]. There is theoretical and empirical evidence that forest dynamics exhibit signatures of environmental variability (; ; ). Usually, neutral models based only on demographic stochasticity tend to overestimate the expected times to extinction for abundant species (), whose temporal fluctuations will also be underestimated (). Incorporating an environmental source of randomness as well as more realistic forms of speciation have made some dynamic aspects of NT more realistic ().

Here we focus on the dynamic aspects of the SAD () which under appropriate assumptions can be studied in detail (; ) and whose predictions can be benchmarked against empirical data. We also review recent progress in modeling dynamic patterns, including the species turnover distribution () as well as species’ persistence-times (or lifetimes) distributions (; ; ; ); see Table .

A. The continuum limit of the master equation

The microscopic description of a system starts by correctly identifying the variables that define all the possible configurations of the system. Having decided how to describe the states of the system, the next step is to consider the transition rates among different states. The configurations of an ecological community can be described by different variables according to different levels of coarsening of the spatiotemporal scales. However, for simplicity, we focus on ecosystems comprising $𝑆$ species and $𝑁$ individuals in total, in which the configurations are specified fully by the variable $𝐧 = (𝑛 1, \dots, 𝑛 𝑆)$ , where $𝑛 𝑖$ is the population of the $𝑖$ th species. Therefore, for the time being, we ignore any spatial effects, which will be considered in Sec. .

If the ecological community at time $𝑡 0$ was in the configuration $𝐧 0$ , then the probability that the system will be in configuration $𝐧$ at time $𝑡 > 𝑡 0$ is $𝑃 (𝐧, 𝑡 | 𝐧 0, 𝑡 0)$ [ $𝑃 (𝐧, 𝑡)$ for brevity]. Assuming a Markovian stochastic dynamics, the evolution of $𝑃 (𝐧, 𝑡)$ is governed by the ME with transition rates $𝑇 [𝐧 | 𝐧']$ from the states $𝐧'$ to $𝐧$ . As seen in Sec. , the specification of a model is tantamount to defining the transition rates. This equation can be recast in another form that allows a simple expansion. Interactions among species change the variable $𝐧$ in steps of size $𝐫 ℓ$ , where $ℓ = 1, 2, \dots, 𝐿$ and $𝐿$ is the total number of distinct possible changes or reactions. For instance, if a species can either decrease or increase its population by $𝑘$ in a given time step, then $𝐿 = 2$ . Species that change by 1 or that remain constant in a given time interval are instead described by $𝐿 = 3$ , and when the total population is conserved, we have $\sum 𝐿 ℓ = 1 𝑟 𝑖, ℓ = 0$ . Of all the possible reactions, however, only a few are usually significant, ecologically relevant, or meaningful. Moreover, unlike chemical reactions these meaningful reactions can be considered irreversible. Therefore, instead of summing over states $𝐧'$ , as in Eq. , we can sum over $𝐿$ different reactions. If $𝑡 ℓ [𝐧]$ is the transition rate of the $ℓ$ th reaction which involves the jumps $𝐫 ℓ$ when species’ populations are $𝐧$ , then we can recast Eq. in the form

\partial 𝑃 (𝐧, 𝑡) \partial 𝑡 = 𝐿 \sum ℓ = 1 {𝑡 ℓ [𝐧 - 𝐫 ℓ] 𝑃 (𝐧 - 𝐫 ℓ, 𝑡) - 𝑡 ℓ [𝐧] 𝑃 (𝐧, 𝑡)} .

(30)

Assuming that $𝑡 ℓ [𝐧]$ is only a function of $𝐱 = 𝐧 / 𝑁$ [this corresponds to Eqs. and with $𝑁 = 𝐽 𝑀$ and Eqs. and with $𝑁 = 𝐽 𝐿$ and $𝑁 ≫ 1$ in both cases], we can write

\partial 𝑃 (𝐱, 𝑡) \partial 𝑡 = 𝐿 \sum ℓ = 1 {𝑡 ℓ [𝐱 - 𝐫 ℓ 𝑁] 𝑃 (𝐱 - 𝐫 ℓ 𝑁, 𝑡) - 𝑡 ℓ [𝐱] 𝑃 (𝐱, 𝑡)},

(31)

where $𝐧$ has been replaced by $𝑁 𝐱$ . This form of the ME suggests that $| 𝐫 ℓ | / 𝑁 ≪ 1$ when $𝑁 ≫ 1$ and hence, a Taylor expansion around $𝐱$ of $𝑡 ℓ$ and $𝑃$ is legitimate, at least far from the boundaries. This is known as the Kramers-Moyal expansion (). Truncating to second order in $𝑁$ , we obtain a Fokker-Planck equation of the form

\partial 𝑃 (𝐱, 𝜏) \partial 𝜏 = - 𝑆 \sum 𝑖 = 1 \partial \partial 𝑥 𝑖 [𝐴 𝑖 (𝐱) 𝑃 (𝐱, 𝜏)] + 1 2 𝑁 𝑆 \sum 𝑖, 𝑗 = 1 \partial 2 \partial 𝑥 𝑖 \partial 𝑥 𝑗 [𝐵 𝑖 𝑗 (𝐱) 𝑃 (𝐱, 𝜏)] + 𝑂 (𝑁 - 2),

(32)

where the time has been rescaled ( $𝜏 = 𝑡 / 𝑁$ ), the functions $𝐴 𝑖 (𝐱)$ have been defined as

𝐴 𝑖 (𝐱) \equiv 𝐿 \sum ℓ = 1 𝑟 𝑖, ℓ 𝑡 ℓ [𝐱],

(33)

and the entries of the matrix $𝐵 (𝐱)$ are

𝐵 𝑖 𝑗 (𝐱) \equiv 𝐿 \sum ℓ = 1 𝑟 𝑖, ℓ 𝑟 𝑗, ℓ 𝑡 ℓ [𝐱] .

(34)

The matrix $𝐵 (𝐱)$ is a $𝑆 \times 𝑆$ symmetric matrix that is positive semidefinite (). For one-step jumps, the entries of the matrices $𝐫 ℓ \otimes 𝐫 ℓ$ can only be 1, $- 1$ , or 0. Unlike the ME , the Fokker-Planck equation governs the evolution of continuous stochastic variables that represent the relative species’ populations under the so-called diffusion approximation. As the total population extends to infinity, we recover the deterministic evolution given by the Liouville equation

\partial 𝑃 (𝐱, 𝜏) \partial 𝜏 = - 𝑆 \sum 𝑖 = 1 \partial \partial 𝑥 𝑖 [𝐴 𝑖 (𝐱) 𝑃 (𝐱, 𝜏)],

(35)

whose solution, for initial condition $𝑃 (𝐱, 0) = 𝛿 (𝐱 - 𝐱 0)$ , has the simple form, where $𝛿$ is the Dirac delta function. This form of $𝑃$ implies that $𝐱 (𝜏)$ , i.e., the relative species’ populations, are the solutions of the deterministic system $˙ 𝑥 𝑖 (𝜏) = 𝐴 𝑖 (𝐱)$ , where $𝑖 = 1, 2, \dots, 𝑆$ ().

So far we have assumed that the system has a fixed and finite total population $𝑁$ and therefore each population cannot exceed $𝑁$ , despite its fluctuations. However, sometimes it is useful to relax this constraint and allow each population to fluctuate within the interval $[0, \infty)$ . In this case, a parameter may or may not exist that allows one to recover the deterministic limit. However, if we assume that $𝑡 ℓ [𝐧]$ are smooth functions of $𝐧$ , and when $𝐧 \to 𝐧 + 𝐫 ℓ$ , they vary little with respect to a characteristic scale of the model (e.g., the total average population), and then, we can treat $𝐧$ as continuous variables and expand $𝑡 ℓ [𝐧 - 𝐫 ℓ] 𝑃 (𝐧 - 𝐫 ℓ, 𝑡)$ around $𝐧$ . Starting from the ME and using $𝐱$ to indicate the continuous populations, we obtain

\partial 𝑃 (𝐱, 𝑡) \partial 𝑡 = - 𝑆 \sum 𝑖 = 1 \partial \partial 𝑥 𝑖 [𝐴 𝑖 (𝐱) 𝑃 (𝐱, 𝑡)] + 12 𝑆 \sum 𝑖, 𝑗 = 1 \partial 2 \partial 𝑥 𝑖 \partial 𝑥 𝑗 [𝐵 𝑖 𝑗 (𝐱) 𝑃 (𝐱, 𝑡)],

(36)

where time has not been rescaled and $𝐴 𝑖 (𝐱)$ and $𝐵 𝑖 𝑗 (𝐱)$ are as in Eqs. and . If we are within a neutral framework with one-step jumps, then $𝑆 = 1$ , $𝐿 = 2$ , and $𝑟 1 = 1, 𝑟 2 = - 1$ . Setting $𝑡 1 [𝑥] \equiv 𝑏 (𝑥)$ (birth rate) and $𝑡 2 [𝑥] \equiv 𝑑 (𝑥)$ (death rate), from Eq. one gets

\partial 𝑃 (𝑥, 𝑡) \partial 𝑡 = \partial \partial 𝑥 [𝑑 (𝑥) - 𝑏 (𝑥)] 𝑃 (𝑥, 𝑡) + 12 \partial 2 \partial 𝑥 2 [𝑑 (𝑥) + 𝑏 (𝑥)] 𝑃 (𝑥, 𝑡) .

(37)

Starting from the discrete formulation of the birth and death rates that we have seen previously, one can derive the expressions of the corresponding continuous rates. Following our earlier discussion, we can write the following general expansion:

𝑏 (𝑥) 𝑥 = 𝑏 1 + 𝑏 0 𝑥 + \dots,

(38a)

𝑑 (𝑥) 𝑥 = 𝑑 1 + 𝑑 0 𝑥 + \dots,

(38b)

where we dropped higher order terms in

1 / 𝑥

. Clearly, Eq. can be applied when species have sufficiently large populations. The constants

𝑏 0

and

𝑑 0

produce a density-dependence effect that causes a rare-species advantage (disadvantage) when

𝑏 0 > 𝑑 0

(

𝑏 0 < 𝑑 0

). This effect has its roots in ecological mechanisms, such as the Janzen-Connell (; ) or Chesson-Warner effects (). However, the presence of net immigration or speciation in a local community () can also produce such density dependence, which is captured here by a mean-field approach. Usually, the skewness of the RSA of various tropical forests and coral reefs indicates a rare-species advantage (), and thus we will use

𝑏 0 > 𝑑 0

. In addition, to simplify the analytical treatment and for parsimony, we choose

𝑏 0 = - 𝑑 0

. Substituting the rates in Eq. into Eq. and setting

𝑏 0 = - 𝑑 0 > 0

, we obtain

\partial 𝑃 (𝑥, 𝑡) \partial 𝑡 = \partial \partial 𝑥 [(𝜇 𝑥 - 𝑏) 𝑃 (𝑥, 𝑡)] + 𝐷 \partial 2 \partial 𝑥 2 [𝑥 𝑃 (𝑥, 𝑡)],

(39)

where $𝐷 = (𝑑 1 + 𝑏 1) / 2$ , $𝜇 = 1 / 𝜏 = 𝑑 1 - 𝑏 1 > 0$ , and $𝑏 = 2 𝑏 0 > 0$ . $𝑃 (𝑥, 𝑡) \equiv 𝑃 (𝑥, 𝑡 | 𝑥 0, 0)$ is the probability density function of finding $𝑥$ individuals in the community at time $𝑡$ , given that at time $𝑡 = 0$ there were $𝑥 0$ . Therefore $\int 𝑛 + Δ 𝑛 𝑛 𝑃 (𝑥, 𝑡) 𝑑 𝑥$ is the fraction of species with a population between $𝑛$ and $𝑛 + Δ 𝑛$ at time $𝑡$ . Equation defines our model for species-rich communities and it governs the time evolution of a community population under the neutral approximation. The deterministic term drives the population to the stationary mean population per species $𝑏 / 𝜇$ and therefore, within our model, the population and the number of species can also fluctuate at stationarity with the ratio of the population to the number of species being fixed on average. Thus, this model is more flexible than the original model by in which “zero-sum dynamics” that fixes the total number of individuals in a community was assumed.

The link established between the FP equation and the ME provides a useful interpretation of the coefficients: $𝜇 = 𝑑 1 - 𝑏 1$ is the imbalance between the per-capita death and birth rates that inexorably drives the ecosystem to extinction in the absence of immigration or speciation. In this model $𝜏$ is the characteristic time scale associated with perturbations to the steady state. When $𝐷 𝜏 ≫ 1$ (as for the tropical forests we have analyzed) we have $𝑑 1 = (2 𝐷 + 𝜏 - 1) / 2 ≃ 𝐷$ , and so $𝐷$ can be thought of as an individual death rate. Finally, $𝑏$ plays the role of an effective immigration (or speciation) rate that prevents the community from becoming extinct.

As expected, different choices of $𝑏 0$ and $𝑑 0$ lead to very good fits of the RSA of various tropical forests. In particular, it is possible to see that the RSA fits are readily improved for large $𝑥$ when $𝑏 0$ and $𝑑 0$ are arbitrary parameters. However, if one introduces a fourth free parameter, the analytical treatment of the dynamics is much more involved. The Fokker-Planck equation with just three parameters is an ideal starting point to understand the dynamics governing species-rich ecosystems in a simplified fashion. Indeed, the agreement between empirical data and the macroecological properties predicted by the model is consistent with the simplifying assumption that tropical forests are close to their steady state. In addition, because the model is not only neutral, but also noninteracting, one may speculate that surviving species are those that were able to reach a steady state of coexistence by minimizing interspecies interactions.

Further insights into the nature of stochasticity can be achieved by writing down Eq. in the equivalent Langevin form, which is an equation for the state variables themselves. Within the Itô prescription (), Eq. is equivalent to the following Langevin equation:

˙ 𝑥 𝑖 (𝑡) = 𝐴 𝑖 (𝐱) + 𝑆 \sum 𝑗 = 1 𝑏 𝑖 𝑗 (𝐱) 𝜉 𝑗 (𝑡),

(40)

where $𝜉 𝑗 (𝑡)$ is Gaussian white noise with a zero mean and correlation $⟨ 𝜉 𝑖 (𝑡) 𝜉 𝑗 (𝑡') ⟩ = 𝛿 𝑖 𝑗 𝛿 (𝑡 - 𝑡')$ , and the matrix $𝑏 (𝐱)$ is defined by $𝐵 (𝐱) = 𝑏 (𝐱) 𝑏 𝑇 (𝐱)$ , where the entries of $𝐵 (𝐱)$ are given in Eq. . This shows that $𝑏 (𝐱)$ may be thought of as the “square root” of $𝐵 (𝐱)$ . However, because any orthogonal transformation of the noise $𝜉 𝑗 (𝑡)$ leaves the mean and the correlation unchanged, the matrix $𝑏 (𝐱)$ is not uniquely determined.

For the present model, the Fokker-Planck equation is equivalent to

˙ 𝑥 (𝑡) = 𝑏 - 𝑥 (𝑡) 𝜏 + \sqrt 𝐷 𝑥 (𝑡) 𝜉 (𝑡),

(41)

where we have instead assumed that the time correlation is $⟨ 𝜉 (𝑡) 𝜉 (𝑡') ⟩ = 2 𝛿 (𝑡 - 𝑡')$ . Equation shows how demographic stochasticity accounts for random ecological drift. In fact, given the relatively large number of individuals of any species, one expects that the detailed nature of the stochastic noise is not important, such that fluctuations are proportional to $\sqrt 𝑥 (𝑡)$ due to the central limit theorem (). Thus, the multiplicative noise in the Langevin equation is not introduced ad hoc and can be justified on the basis of more general considerations.

B. Boundary conditions

The stochastic dynamics described by Eq. governs the evolution of the population when it is strictly positive. Assuming that $𝑏 > 0$ and starting at some positive value, the process cannot reach negative values because the drift brings it to a positive value $𝑏 𝑑 𝑡$ in the time interval $𝑑 𝑡$ . Therefore, the process is always non-negative with the origin acting as a singular boundary, which must be specified on the basis of ecological considerations, i.e., one should define what happens when the random ecological drift occasionally leads to a vanishing population. The most frequent ecological situations are given by the following boundaries.

• We use reflecting boundary conditions at $𝑥 = 0$ to describe ecological communities in which a vanishing population is immediately replaced by a new individual (belonging to either a new or old species). For instance, this may happen when an ecosystem is coupled to a metacommunity that extends across large spatial scales. With these boundary conditions, one can describe ecosystems in which biodiversity is continuously sustained by the net immigration of new individuals of new or old species. This prevents an ecosystem from becoming extinct and can finally achieve a nontrivial steady state. In the case of the model defined by Eq. , reflecting boundaries are obtained by setting the flux of the probability distribution as zero at $𝑥 = 0$ , although there are some equations for which such boundaries cannot be arbitrarily set, such as Eq. ().
• A second possibility arises when a community can lose individuals without any replacement or a net emigration flux of individuals from the ecosystem exists. One can then describe the system by introducing absorbing boundary conditions at $𝑥 = 0$ . These constraints force ecosystems to march inexorably to extinction and the final steady state corresponds to the complete extinction of species, i.e., $l i m 𝑡 \to \infty 𝑃 (𝑥, 𝑡) = 𝛿 (𝑥)$ . Usually, absorbing boundaries are obtained by setting the probability distribution at $𝑥 = 0$ equal to zero, yet some equations do have norm decreasing solutions, such as Eq. (), which do not vanish at $𝑥 = 0$ .

It is possible to define other kinds of boundaries, according to different ecological behaviors when the population consists of just a few individuals. However, the reflecting and absorbing boundaries capture the most interesting situations in ecology.

As we have alluded, the reflecting and absorbing solutions of the Fokker-Planck equation have rather unusual properties. This is essentially due to the fact that the diffusion term $𝐷 𝑥$ vanishes at $𝑥 = 0$ . Obviously, it is always possible to rewrite a vanishing diffusion term at zero into a nonvanishing one with a suitable change of variables, yet the same change of variables leads to a singular drift term at the boundary. Thus, equations with a vanishing diffusion term or with a singularity of the drift term at the boundary in general share similar unusual properties. In fact, there exist regions in the parameter space of Eq. in which one cannot arbitrarily fix a boundary at zero, and the solution is completely defined once the initial condition is prescribed. This is because the probability of accessing the origin depends on the values of the parameters. The problem can be studied by introducing the probability to reach the origin for the first time. It is possible to show () that this first-passage probability to the origin is zero when $𝑏 / 𝐷 > 1$ and, therefore, in this case one cannot obtain any solution with absorbing boundaries and the solution can only be reflecting at $𝑥 = 0$ . However, a (unique) norm decreasing solution exists when $0 < 𝑏 / 𝐷 < 1$ , which corresponds to the regular absorbing solution but that does not vanish at $𝑥 = 0$ . Also, in the same region of parameters, a reflecting solution exists with an integrable singularity at $𝑥 = 0$ . All ecosystems that we have studied so far are well described by $𝑏$ and $𝐷$ within this parameter region.

C. Stationary and time-dependent solutions

An ecosystem described by Eq. can reach a nontrivial steady state when setting reflecting boundary conditions at $𝑥 = 0$ . In this case, the normalized stationary solution is the gamma distribution (), i.e.,

𝑃 0 (𝑥) = (𝐷 𝜏) - 𝑏 / 𝐷 Γ (𝑏 / 𝐷) 𝑥 𝑏 / 𝐷 - 1 𝑒 - 𝑥 / 𝐷 𝜏,

(42)

where $Γ (𝑥)$ is the gamma function (). Note the distribution exists only when $𝑏 > 0$ , i.e., when the ecosystem has a net immigration rate of species. $𝑃 0 (𝑥)$ is the probability density function (pdf) of the relative species abundance at stationarity and it gives the probability that a given arbitrary species present in the ecosystem has $𝑥$ individuals. Fitting the empirical RSAs allows one to estimate two combinations of the parameters, i.e., $𝑏 / 𝐷$ and $𝐷 𝜏$ , with $𝜏$ being the characteristic temporal scale to approach stationarity. Interestingly, for $0 < 𝑏 / 𝐷 ≪ 1$ one obtains the continuum approximation of the celebrated Fisher log series (), the RSA distribution of a metacommunity.

One can analytically solve Eq. at any time with arbitrary initial conditions, so that the evolution of the population can be traced even far from stationary conditions. The time-dependent solution with reflecting boundary conditions ( $𝑏 > 0$ ) and initial population $𝑥 0$ is

𝑃 𝑟 (𝑥, 𝑡 | 𝑥 0, 0) = (1 𝐷 𝜏) 𝑏 / 𝐷 𝑥 𝑏 / 𝐷 - 1 𝑒 - 𝑥 / 𝐷 𝜏 [(1 / 𝐷 𝜏) 2 𝑥 0 𝑥 𝑒 - 𝑡 / 𝜏] 1 / 2 - 𝑏 / 2 𝐷 1 - 𝑒 - 𝑡 / 𝜏 \times e x p [- (1 / 𝐷 𝜏) (𝑥 + 𝑥 0) 𝑒 - 𝑡 / 𝜏 1 - 𝑒 - 𝑡 / 𝜏] 𝐼 𝑏 / 𝐷 - 1 [(2 / 𝐷 𝜏) \sqrt 𝑥 0 𝑥 𝑒 - 𝑡 / 𝜏 1 - 𝑒 - 𝑡 / 𝜏],

(43)

where $𝐼 𝜈 (𝑧)$ is the modified Bessel function of the first kind (). Because $𝐼 𝜈 (𝑧) ≃ (𝑧 / 2) 𝜈 / Γ (𝜈 + 1)$ when $𝑧 \to 0 +$ , for large time intervals Eq. converges to $𝑃 0 (𝑥)$ in Eq. . By setting $𝑥 = 𝐷 𝑦 2 / 4$ ( $𝑦 > 0$ ), when $𝑏 / 𝐷 = 1 / 2$ the process can be mapped into the simpler Ornstein-Uhlenbeck process () with a boundary (reflecting or absorbing) at $𝑥 = 0$ , and when $𝑏 / 𝐷 = 1$ it can be mapped into the (reflecting) Rayleigh process ().

From the Fokker-Planck equation , it is easy to derive the evolution of the mean population per species. It is simply

⟨ 𝑥 (𝑡) ⟩ = 𝑏 𝜏 + (𝑥 0 - 𝑏 𝜏) 𝑒 - 𝑡 / 𝜏,

(44)

where $⟨ 𝑥 (0) ⟩ = 𝑥 0$ is the initial population per species in the community. The mean number of individuals per species converges to $𝑏 𝜏$ at stationarity, with standard deviation $𝜏 \sqrt 𝑏 𝐷$ . The autocovariance $𝜅 (𝑡) = ⟨ 𝑥 (𝑡) 𝑥 (0) ⟩ - ⟨ 𝑥 (𝑡) ⟩ ⟨ 𝑥 (0) ⟩$ under stationary conditions is given by

𝜅 (𝑡) = \int \infty 0 [𝑥 𝑥 0 - ⟨ 𝑥 (𝑡) ⟩ ⟨ 𝑥 (0) ⟩] 𝑃 𝑟 (𝑥, 𝑡 | 𝑥 0, 0) 𝑃 0 (𝑥 0) 𝑑 𝑥 𝑑 𝑥 0,

(45)

where $𝑃 𝑟 (𝑥, 𝑡 | 𝑥 0, 0)$ is the reflecting solution in Eq. and $𝑃 0 (𝑥 0)$ is the stationary solution in Eq. . However, $𝜅 (𝑡)$ can be more rapidly obtained from Eq. . In fact, since in the Itô prescription $𝜉 (𝑡)$ and $𝑥 (𝑡)$ are uncorrelated, one gets $⟨ \sqrt 𝑥 (𝑡) 𝜉 (𝑡) ⟩ = ⟨ \sqrt 𝑥 (𝑡) ⟩ ⟨ 𝜉 (𝑡) ⟩ = 0$ because we have assumed that $⟨ 𝜉 (𝑡) ⟩ = 0$ . Therefore the evolution equation for $⟨ 𝑥 (𝑡) 𝑥 (0) ⟩$ becomes $\partial 𝑡 ⟨ 𝑥 (𝑡) 𝑥 (0) ⟩ = 𝑏 ⟨ 𝑥 (0) ⟩ - ⟨ 𝑥 (𝑡) 𝑥 (0) ⟩ / 𝜏$ , whose solution is $⟨ 𝑥 (𝑡) 𝑥 (0) ⟩ = (𝑏 𝜏) 2 + 𝑏 𝐷 𝜏 2 𝑒 - 𝑡 / 𝜏$ . Hence, under stationary conditions, we have $⟨ 𝑥 ⟩ = 𝑏 𝜏$ , $⟨ 𝑥 2 ⟩ = ⟨ 𝑥 ⟩ 2 + 𝑏 𝐷 𝜏 2$ , and then $𝜅 (𝑡) = 𝑏 𝐷 𝜏 2 𝑒 - 𝑡 / 𝜏$ for $𝑡 > 0$ . Therefore, the time correlation function is a simple exponential. Because the equation has only one temporal scale, $𝜏$ is also the characteristic relaxation time of the process, i.e., the temporal scale upon which the system at stationarity recovers from a small perturbation.

D. The species turnover distribution

According to the NT, the turnover of ecological communities reflects their continuous reassembly through immigration and emigration and local extinction and speciation. Species’ histories overlap by chance due to stochasticity, yet their lifetimes are finite and distributed according to the underlying governing process. Nontrivial stationary communities are reached because old species are continuously replaced by new species, bringing about a turnover of species that can be studied and modeled within our framework.

To measure species turnover one usually considers the population of a given species at different times, then studying the temporal evolution of their ratios. For an ecosystem close to stationarity, one can look at the distribution $𝒫 (𝜆, 𝑡)$ , i.e., the probability that at time $𝑡$ the ratio $𝑥 (𝑡) / 𝑥 (0)$ is equal to $𝜆$ , where $𝑥 (𝑡)$ and $𝑥 (0)$ are the population of a given species at time $𝑡 > 0$ and $𝑡 = 0$ , respectively. Thus, the STD $𝒫 (𝜆, 𝑡)$ by definition is

𝒫 (𝜆, 𝑡) = ⟨ 𝛿 (𝜆 - 𝑥 / 𝑥 0) ⟩ = \int \infty 0 𝑑 𝑥 0 \int \infty 0 𝑑 𝑥 𝑃 (𝑥, 𝑡 | 𝑥 0, 0) 𝑃 0 (𝑥 0) 𝛿 (𝜆 - 𝑥 / 𝑥 0) .

(46)

Here $𝑃 (𝑥, 𝑡 | 𝑥 0, 0)$ is either the reflecting or absorbing solution defined in Eqs. and , and $𝑃 0 (𝑥)$ is given in Eq. , and $𝜆 > 0$ and $𝑡 > 0$ .

In Eq. , one can use the time-dependent reflecting or absorbing solution according to different ecological dynamics. We should use the reflecting solution when not concerned with the extinction of the species present at the initial time point and especially when accounting for any new species introduced through immigration or speciation. The expression for the reflecting STD can be found in .

One can show that it has the following power-law asymptotic behavior for a fixed $𝑡$ :

𝒫 (𝜆, 𝑡) \approx {𝑘 1 (𝑡) 𝜆 𝑏 / 𝐷 + 1 for 𝜆 ≫ 1, 𝑘 2 (𝑡) 𝜆 𝑏 / 𝐷 - 1 for 𝜆 ≪ 1,

(47)

where the functions $𝑘 1$ and $𝑘 2$ are independent of $𝜆$ . Customarily, rather than the random variable $𝜆$ , ecologists study $𝑟 = l o g (𝜆)$ that is distributed according to $𝑔 (𝑟, 𝑡) = 𝑒 𝑟 𝒫 (𝑒 𝑟, 𝑡)$ and that can be compared to empirical data in Fig. . We obtain an estimate of $𝜏$ , the characteristic time scale for the BCI forest, which is around $3500 \pm 1000 y r s$ [for trees with $> 10 c m$ of stem diameter at breast height (dbh)] and $2900 \pm 1100 y r s$ (for trees with $> 1 c m$ dbh), where the broad uncertainty is due to the fact that the data are sampled over relatively short time intervals.

FIG. 4.

STD for the interval 1990–1995 in the BCI forest. The main panel shows the results for individuals of more than 10 cm dbh, and the inset the results for individuals of more than 1 cm dbh (Center for Tropical Forest Science website). The black line represents the analytical solution given by Eq. . From .

These fits not only provide direct information about the time scale of evolution but also they underline the importance of rare species in the STD. $𝑏 / 𝐷$ is closely tied to the distribution of rare species and in fact $𝑃 0 (𝑥) \sim 𝑥 𝑏 / 𝐷 - 1$ for $𝑥 ≪ 𝐷 𝜏$ . Alternatively, the dependence of the STD on $𝑏 / 𝐷$ and $𝜏$ suggests that at any fixed time $𝑡 > 0$ rare species are responsible for the shape of the STD.

E. Persistence or lifetime distributions

A theoretical framework to study and analyze persistence or extinctions of species in ecosystems allows one to understand the link between environmental changes [like habitat destruction or climate change (; ; ; ; )] and the increasing number of threatened species.

The persistence or lifetime $𝜏$ of a species is defined as the time interval between its emergence and its local extinction (see Fig. ) within a given geographic region (; ). In statistical physics, this is known as the distribution of first-passage times to zero of the stochastic processes describing the species-abundance dynamics. According to the neutral theory of biodiversity, species can span very different lifetime intervals and thus, at a local scale, persistence times are largely controlled by ecological processes like random drift, dispersal, and immigration.

FIG. 5.

Schematic representation of persistence time (or lifetime) of a species $𝜏$ and survival times $𝜏 𝑠$ , defined as the time to local extinction of a species randomly sampled among the observed assemblages at a certain time $𝑇$ . From .

1. Discrete population dynamics

The simplest baseline model to study the persistence-time distribution is a random walk in the species abundance $𝐧$ , i.e., ME with $𝑏 𝑛 = 𝑑 𝑛 = 𝑐$ and absorbing boundary condition in $𝑛 = 0$ . According to this scheme, local extinction is equivalent to a random walker’s first passage to zero and thus, the resulting persistence-time distribution has a power-law decay with exponent $3 / 2$ (; ).

A further step in modeling lifetime distributions can be made by taking into account birth, death, and speciation (; ; ; ) through a mean-field scheme of the voter model with speciation (; ) (see description in Sec. ), i.e., ME with birth and death rates given by Eqs. and in the large $𝐽 𝑀$ limit, $𝑏 (𝑛) = (1 - 𝜈) 𝑛 / 𝐽 𝑀$ and $𝑑 (𝑛) = 𝑛 / 𝐽 𝑀$ for $𝑛 \geq 0$ . The corresponding persistence or lifetime distribution is given by $𝑝 𝜏 (𝑡) = - 𝑑 𝑃 (0, 𝑡) / 𝑑 𝑡$ , where $𝑃 (0, 𝑡)$ is the probability that a species has a zero population at time $𝑡$ . The asymptotic behavior of the resulting persistence-time distribution [i.e., $𝑝 𝜏 (𝑡)$ ] exhibits a power-law scaling modified by an exponential cutoff ():

𝑝 𝜏 (𝑡) \propto 𝑡 - 2 𝑒 - 𝜈 𝑡,

(48)

for $𝑡$ greater than some lower cutoff value. Here $𝑡$ is measured in units of $𝐽 𝑀$ and $𝜈$ is now a speciation rate rather than a probability [Eq. is derived in Appendix ]. Note that, in this context, $𝑝 𝜏 (𝑡)$ has a characteristic time scale $1 / 𝜈$ for local extinctions determined by the speciation or migration rate. The general case when none of the coefficients is zero and are given by Eqs. can also be solved. In this case, $𝑝 𝜏$ displays a crossover from the $𝑡 - 3 / 2$ to the $𝑡 - 2$ behavior at a certain characteristic time and finally, an exponential decay beyond yet another characteristic time scale (). It has been shown numerically (; ) that for the spatial voter model in dimensions $𝑑 = 1$ , 2, 3, the exponent of the power law in Eq. changes to $𝑡 - 𝛼$ , with $𝛼 < 2$ depending on the topological structure underlying the voter model. In particular, $𝛼 = 3 / 2$ in $𝑑 = 1$ and $𝛼 = 2$ for any $𝑑 \geq 3$ , whereas in $𝑑 = 2$ one gets $𝑝 𝜏 (𝑡) \propto 𝑡 - 2 𝑒 - 𝜈 𝑡 l n 𝑡$ as shown in and references therein.

We have seen in Sec. that the RSA pattern does not depend on the biological details of the ecosystem under analysis. Thus, one may wonder if the persistence-time distribution is also a universal macroecological pattern. Indeed, it has been shown (; ) that the power law with an exponential cutoff shape predicted for the persistence-time distribution by the NT (see Appendix ) is common to very different types of ecosystems (see Fig. ). Other exact formulas for species ages and species lifetimes have been proposed in neutral () and non-neutral models (; ).

FIG. 6.

Comparison between persistence empirical distributions for North American breeding birds, Kansas grasslands, New Jersey BSS forest, and an estuarine fish community and the corresponding theoretical species persistence-times pdfs. The circles and solid lines show the observational distributions and fits, respectively. The finiteness of the time window $Δ 𝑇 𝑤$ imposes a cutoff in the maximum observable persistence time and thus, only lifetimes where $𝜏 < Δ 𝑇 𝑤$ have been considered and the theoretical predictions have been adjusted appropriately (see Appendix ). Adapted from .

Another interesting and related quantity is the survival distribution $𝑃 𝜏 𝑠$ defined as the probability that a species randomly sampled from the community at stationarity is still present in the community after a time $𝑡$ . This quantity depends on the initial conditions as $𝑃 𝜏 𝑠 (𝑡) = \int 𝑝 𝜏 (𝑡 | 𝑛 0) 𝑃 0 (𝑛 0) 𝑑 𝑛 0$ , where $𝑝 𝜏 (𝑡 | 𝑛 0)$ is the lifetime distribution for a species that initially has a population $𝑛 0$ . Assuming that the stationary distribution of population abundances is the Fisher log series, then $𝑃 𝜏 𝑠 (𝑡) \sim 𝑡 - 1$ when $𝑡 ≪ 𝑡 *$ , whereas $𝑃 𝜏 𝑠 (𝑡) \sim 𝑒 - 𝑡 / 𝑡 *$ when $𝑡 ≫ 𝑡 *$ (; ). It can be shown that this asymptotic behavior of the survivor distribution is valid regardless of the functional shape of the birth and death rates $𝑏 (𝑛)$ and $𝑑 (𝑛)$ involved in the ME driving the evolution of $𝑃 (𝑛, 𝑡)$ ().

2. Continuum limit

In the continuum limit, a crucial distribution for the analysis of species’ extinction is the time-dependent solution of Eq. with absorbing boundaries at $𝑥 = 0$ . We refer to for its complete derivation. This probability distribution exists only when $𝑏 < 𝐷$ and is given by

𝑃 𝑎 (𝑥, 𝑡 | 𝑥 0, 0) = (𝐷 𝜏) - 1 1 - 𝑒 - 𝑡 / 𝜏 e x p [- (1 / 𝐷 𝜏) (𝑥 + 𝑥 0 𝑒 - 𝑡 / 𝜏) 1 - 𝑒 - 𝑡 / 𝜏] \times (𝑥 𝑥 0 𝑒 𝑡 / 𝜏) 𝑏 / 2 𝐷 - 1 / 2 𝐼 1 - 𝑏 / 𝐷 [(2 / 𝐷 𝜏) \sqrt 𝑥 0 𝑥 𝑒 𝑡 / 𝜏 𝑒 𝑡 / 𝜏 - 1] .

(49)

Note that Eq. is finite at $𝑥 = 0$ but $l i m 𝑥 \to 0 + 𝑃 𝑎 (𝑥, 𝑡 | 𝑥 0, 0) \neq 0$ .

The lifetime distribution can be calculated analytically using Eq. , i.e., $𝑝 𝜏 (𝑡) = - (𝑑 / 𝑑 𝑡) \int \infty 0 𝑃 𝑎 (𝑥, 𝑡 | 1, 0) 𝑑 𝑥$ . It can be shown that the lifetime distribution calculated in this way displays the same asymptotic behavior as Eq. , although the functional form is different (). In the continuum limit, one can also obtain the mean extinction time ():

⟨ 𝑡 ⟩ = 𝜏 \int 10 𝑥 - 𝑏 / 𝐷 - 1 1 - 𝑥 𝑑 𝑥 = - 𝜏 [𝛾 + 𝜓 (1 - 𝑏 / 𝐷)],

(50)

which depends only on $0 < 𝑏 / 𝐷 < 1$ and $𝜏$ . Here $𝛾 = 0.577 \dots$ is the Euler constant and $𝜓 (𝑧) = Γ' (𝑧) / Γ (𝑧)$ is the logarithmic derivative of the gamma function (). Equation can be used to fit the RSA of various tropical forests yielding $⟨ 𝑡 ⟩ / 𝜏 = 1.94$ , 1.67, 0.67, 0.95, and 1.38 for Yasuni, Lambir, Sinharaja, Korup, and Pasoh, respectively (). These time scales are in accordance with the estimates of extinction times presented elsewhere () and it is quite interesting that $⟨ 𝑡 ⟩ / 𝜏$ depends on $𝑏 / 𝐷$ only, which can be calculated from the steady-state RSA without the need for dynamic data. The values of $𝑏 / 𝐷$ obtained from various tropical forests () suggest that $⟨ 𝑡 ⟩ ≃ 𝜏$ , although this is not built into the model. In general, if $𝜏 ≫ ⟨ 𝑡 ⟩$ , extinction would be much faster than recovery and the ecosystem will not reach a steady state. However, if $𝜏 ≪ ⟨ 𝑡 ⟩$ the ecosystem would recover from external disturbances very rapidly with respect to the extinction time and, therefore, it would be very robust. This would leave little room for the action of evolution. Therefore, the fact that $⟨ 𝑡 ⟩ ≃ 𝜏$ suggests that ecosystems at stationarity might be marginally stable—not so stable that they are frozen in time and not so fragile that they are prone to extinction. From estimates of the model parameters $𝑏 / 𝐷$ , $𝜏$ and thus predictions of $⟨ 𝑡 ⟩$ , many biological and ecological features of the ecosystem may be understood.

IV. Neutral Spatial Models and Environmental Fragmentation

So far we have considered models, which assume that all individuals, at a given time, experience essentially the same conditions. They live in well-mixed habitats in which environmental heterogeneity and spatial distance are not important. These models are conceptually and mathematically simple and this is the basic reason why they have been widely used. However, the spatial structure of ecosystems is able to control the shape of many patterns and it is also a critical factor for understanding species’ coexistence (). In fact, models with well-mixed populations predict that the best competitor is able to displace a pool of species competing for the same limiting resource. There is empirical and theoretical evidence that space plays a crucial role in maintaining species diversity in communities with a single limiting resource, as demonstrated in the grasslands of the Cedar Creek Natural History Area (USA) ().

Space is therefore an essential element for understanding the organization of an ecosystem and most empirical observations are spatial. Often, the dynamics and composition of a community cannot be disentangled from its spatial aspects. Unfortunately, it is very difficult to derive analytical predictions for spatial (stochastic) models. The main difficulty is due to the fact that these are out-of-equilibrium models ().

Spatial effects can be incorporated into the theoretical framework, with the ME of Eq. remaining formally the same by considering the index $𝑖$ in $𝑛 𝑖$ as a composite index $𝑖 = (𝛼, 𝐫)$ , where $𝛼$ identifies the species and $𝐫$ identifies its spatial location. The set of all spatial locations $𝐫$ will be denoted by $Λ$ . Now the transition rates should take into account dispersal from nearby locations. To simplify, we will use the notation $𝑛 𝑖 (𝐫)$ , where $𝑖$ indicates one of the $𝑆$ species and $𝐫$ a spatial location.

At present, no coherent spatial neutral theory exists but rather there is a collection of models and techniques that can explain some spatial patterns. In this section, we review some of those approaches.

The relationship between the number of species and the area sampled is probably one of the oldest quantities studied in ecology (). referred to it as “one of community ecology’s few genuine laws.” The SAR is defined as the average number of species $⟨ 𝑆 (𝐴) ⟩$ sampled in an area $𝐴$ .

postulated a power-law relationship $⟨ 𝑆 (𝐴) ⟩ = 𝑐 𝐴 𝑧$ . Empirical curves show an inverted S shape () (see Fig. ), with a linear behavior at small and large areas, and a power law with an exponent $𝑧$ at intermediate scales. This behavior seems to be pervasive and has been reported for distinct ecosystems.

FIG. 7.

Triphasic shape of the species-area relationship. The left panel shows the three behaviors on different scales. At a local scale the relationship is linear, becoming a power-law relationship at the regional scale and returning to linear at very large intercontinental scales. The right panel shows empirical data of species diversity at the regional scale (from the BCI forest). From .

Despite some notable exceptions (), the value of the exponent $𝑧$ has attracted the most attention in the studies of SAR. This value of $𝑧$ is far from universal, ranging from 0.1 to 0.5 () and showing dependence on latitude (), body mass, taxa, and general environmental conditions (; ). The exponent $𝑧$ is interpreted as a measure of biodiversity.

Several models tried to reproduce the empirical behavior and a simple but useful assumption involves considering different species as independent realizations of the same process (). It is important to note that this assumption is stronger than neutrality, because neutrality does not imply independence. Under these assumptions the SAR is given by Eq. .

The EAR is defined as the number of species that are completely contained (i.e., endemic) in a given area (see Table ). It is not generally simply related to the SAR (). If the species are considered as independent realizations of a unique process, then

⟨ 𝐸 (𝐴) ⟩ = 𝑆 t o t 𝑃 0 (𝐴 𝑐),

(51)

where $𝑆 t o t$ is the total number of species in the system, $𝑃 0 (𝐴 𝑐)$ is the probability that a species is not present in $𝐴 𝑐$ , the complement of $𝐴$ , and that it is therefore completely contained in $𝐴$ .

The $𝛽$ diversity is a spatially explicit measure of biodiversity. A simple and useful measure of $𝛽$ diversity is the similarity index, i.e., the fraction of common species shared between different locations. It can also be defined as the probability $𝐹 (𝑟)$ that two individuals at a distance $𝑟$ are conspecific (). This quantity may be related to the two-point correlation function $𝐺 𝑖 𝑗 (𝑟)$ . Under the assumption of translational invariance, this latter is defined as

𝐺 𝑖 𝑗 (𝑟) ≔ ⟨ 𝑛 𝑖 (𝑥) 𝑛 𝑗 (𝑦 + 𝑟) ⟩ = 1 𝑆 [1 𝑉 \sum 𝑥 \sum 𝑦 𝑛 𝑖 (𝑥) 𝑛 𝑗 (𝑦) 𝛿 𝑖 𝑗 𝛿 (∥ 𝑥 - 𝑦 ∥ - 𝑟)],

(52)

where $𝑛 𝑖 (𝑥)$ is the number of individuals of the species $𝑖$ in the location $𝑥$ , $∥ 𝑥 - 𝑦 ∥$ is the distance between $𝑥$ and $𝑦$ , $𝑉$ is the number of site locations, and $𝛿 (∥ 𝑥 - 𝑦 ∥ - 𝑟)$ selects only pairs at a distance $𝑟$ . $𝐹 (𝑟)$ is the probability that two individuals at a distance $𝑟$ are conspecific, i.e., it is the ratio between the number of pairs of individuals belonging to the same species at a distance $𝑟$ and the total number of pairs of individuals at a distance $𝑟$ . Therefore we obtain

𝐹 (𝑟) = \sum 𝑖 \sum 𝑥 \sum 𝑦 𝑛 𝑖 (𝑥) 𝑛 𝑖 (𝑦) 𝛿 (∥ 𝑥 - 𝑦 ∥ - 𝑟) \sum 𝑖 \sum 𝑗 \sum 𝑥 \sum 𝑦 𝑛 𝑖 (𝑥) 𝑛 𝑗 (𝑦) 𝛿 (∥ 𝑥 - 𝑦 ∥ - 𝑟)

(53)

that can be rewritten as

𝐹 (𝑟) = \sum 𝑖 ⟨ 𝑛 𝑖 (𝑥) 𝑛 𝑖 (𝑥 + 𝑟) ⟩ \sum 𝑖 \sum 𝑗 ⟨ 𝑛 𝑖 (𝑥) 𝑛 𝑗 (𝑥 + 𝑟) ⟩ = \sum 𝑖 𝐺 𝑖 𝑖 (𝑟) \sum 𝑖 \sum 𝑗 𝐺 𝑖 𝑗 (𝑟) .

(54)

If the spatial positions of different species are independent, then $𝐺 𝑖 𝑗 (𝑟) = ⟨ 𝑛 (𝑥) ⟩ 2 ≔ 𝜌 2$ , and by defining $𝐺 𝑖 𝑖 (𝑟) = 𝐺 2 (𝑟)$ , the $𝛽$ diversity reads

𝐹 (𝑟) = 𝐺 2 (𝑟) (𝑆 - 1) 𝜌 2 + 𝐺 2 (𝑟) .

(55)

Most theoretical work consists of attempts to relate these ecological quantities with other spatial and nonspatial observable factors. For instance, a typical problem is to calculate the SAR knowing the $𝛽$ diversity and the RSA over a global scale. One of the future challenges will be to relate and predict spatial patterns on different scales (upscaling and downscaling), having only local information on one or more patterns.

A. Phenomenological models

Phenomenological models do not assume any microscopic dynamics but they are rather based on a given phenomenological distribution of individuals in space. The simplest assumption is to consider individuals at random positions in space (; ).

This null model, usually known as a “random placement” model, can be used to obtain predictions for the SAR and EAR having the RSA or the SAD as an input. Even though this assumption is not realistic, the random placement model turns out to be very useful to capture the relevant aspects of the relationship between RSA and SAR, and it also gives reasonable predictions that can be benchmarked against empirical data. In addition, this assumption also allows direct connections between SAR and EAR to be formulated ().

Consider $𝑆 (𝐴 0)$ species in a region of total area $𝐴 0$ . Species $𝑖$ has an abundance $𝑛 𝑖$ and the $𝑁 = \sum 𝑖 𝑛 𝑖$ individuals are uniformly distributed at random in the area $𝐴 0$ . If we observe a subregion of area $𝐴$ , the probability of observing a particular individual is simply $𝐴 / 𝐴 0$ , while the probability of not observing it is $1 - 𝐴 / 𝐴 0$ . The probability of not observing species $𝑖$ will then be $(1 - 𝐴 / 𝐴 0) 𝑛 𝑖$ , given the fact that the positions of individuals are independent. We can then obtain the average number of species observed in an area $𝐴$ as

⟨ 𝑆 (𝐴) ⟩ = 𝑆 (𝐴 0) - 𝑆 \sum 𝑖 = 1 (1 - 𝐴 𝐴 0) 𝑛 𝑖,

(56)

where $𝑆 (𝐴 0) = 𝑆 t o t$ is the total number of species in the system.

The simple framework of the random placement model also allows the EAR to be calculated. Using Eq. we obtain

⟨ 𝐸 (𝐴) ⟩ = 𝑆 \sum 𝑖 = 1 (𝐴 𝐴 0) 𝑛 𝑖 .

(57)

Despite the simplicity of the approximation, the EAR evaluated using random placement captures the quantitative behavior of several observed ecosystems ().

Under random placement assumptions, one can obtain the EAR from the SAR and vice versa. To calculate the EAR in Eq. , we calculated the number of species with zero individuals in the area complementary to that of interest. This number is equal to the difference between the number of species in the whole area and the number of species in the complementary area. The complementary area has a nontrivial shape and under general assumptions this quantity is not easy to calculate. Under the random placement assumption, the number of species in the complementary area is the SAR of the complementary area. The EAR is therefore

⟨ 𝐸 (𝐴) ⟩ = 𝑆 (𝐴 0) - 𝑆 (𝐴 0 - 𝐴) .

(58)

presented a careful analysis of the reliability of this extrapolation to predict the empirical EAR and showed that the random placement approximation describes the empirical data well.

One can obtain a closed form expression for the EAR and SAR by starting with a RSA distribution. In the case of a Fisher log series [see Eq. ], the SAR reads

⟨ 𝑆 (𝐴) ⟩ = 𝑆 (𝐴 0) - \infty \sum 𝑛 = 1 𝜃 𝑟 𝑛 𝑛 (1 - 𝐴 𝐴 0) 𝑛 = 𝜃 l o g (1 + 𝑟 1 - 𝑟 𝐴 𝐴 0),

(59)

which follows from the observation that $𝑆 (𝐴 0)$ is equal to $- 𝜃 l o g (1 - 𝑟)$ of Eq. . The same calculation can be performed for the EAR, obtaining

⟨ 𝐸 (𝐴) ⟩ = - 𝜃 l o g (1 - 𝑟 𝐴 𝐴 0) .

(60)

In real ecosystems, individuals are not distributed uniformly in space but rather, due to dispersal limitation, individuals of the same species tend to be clustered. This is confirmed by empirical $𝛽$ -diversity plots (). The phenomenological approach of random placement can be generalized to a nonuniform distribution of individuals in space, taking into account the empirically observed spatial clustering of conspecific individuals (). Individuals are distributed in space via a Poisson cluster process (PCP). In a PCP, the centers of the cluster of points are uniformly distributed in space. A random number of points are distributed around each center according to a given spatial kernel. The process depends on two distributions: the number of points in each cluster and the spatial kernel. It is possible to show () that these two distributions may be related to the RSA and the $𝛽$ diversity. An analytical formula for the SAR and the EAR can therefore be obtained given the RSA and the $𝛽$ diversity. This approach reproduces the S-inverted shape observed in empirical systems, showing that this shape can be explained simply in terms of spatial correlations of conspecific individuals.

B. Spatial stochastic processes

There are two possible ways to include space in a neutral stochastic model, either implicitly or explicitly.

1. Spatially implicit model

Spatially implicit models are based on the observation that one can relate the sample area $𝐴$ to the total number of individuals $𝐽$ (), given that $𝐽 = 𝜌 𝐴$ , where $𝜌$ is the density. One can therefore obtain species-area curves in nonspatial models, looking at the scaling of the number of species $𝑆$ with the number of individuals $𝐽$ . Spatial implicit models can thus be thought of as mean-field “well-mixed” models, where one can neglect dispersal limitation.

In a metacommunity the number of species $⟨ 𝜙 (𝑛) ⟩$ of a population $𝑛$ is a Fisher log series (see Sec. )

⟨ 𝜙 (𝑛) ⟩ = 𝜃 𝑟 𝑛 𝑛,

(61)

where $𝑟$ is the ratio of the birth to death rate. The average number of individuals in the metacommunity is $𝐽 𝑀 = 𝜃 𝑟 / (1 - 𝑟)$ and the expected number of species can be easily computed

⟨ 𝑆 (𝐽 𝑀) ⟩ = 𝜃 \infty \sum 𝑛 = 1 𝑟 𝑛 𝑛 = - 𝜃 l o g (1 + 𝑟) = 𝜃 l o g (1 + 𝐽 𝑀 𝜃) .

(62)

This result corresponds exactly to that found using the random placement in Eq. . At small sample sizes, the number of species $𝐽 𝑀 ≪ 𝜃$ is equal to the number of individuals $⟨ 𝑆 (𝐽 𝑀) ⟩ = 𝐽 𝑀$ . In other words when small areas are sampled the individuals belong to different species and the number of species grows along with the number of individuals. With larger sample sizes, the number of species $𝐽 𝑀 ≫ 𝜃$ grows logarithmically with the number of individuals. This approach allows one to calculate the SAR directly from the RSA distribution and it is clearly applicable to any RSA distribution.

Real ecosystems are of course spatially explicit, but one might wonder how spatially implicit models or, more generally, models that do not consider space explicitly, are predictive and how their parameters are related to spatially explicit ones. A way to assess this is to measure the efficacy of nonspatial models in predicting the behavior of spatially explicit models (). As expected, nonspatial models have a good predictive power when the dispersal lengths are sufficiently large and they are particularly good in predicting nonspatial patterns such as the RSA.

2. Spatially explicit model

Spatially explicit models are typically defined as birth-death-diffusion processes. A model is fully specified given a ME and can be obtained in several ways, i.e., it is possible to write several different MEs that include space in a neutral model. The ME is not tractable analytically and one has to introduce approximations in order to get analytical results. Spatially explicit models are particularly difficult to solve because of the lack of detailed balance (see Sec. ).

The voter model () was originally introduced to describe opinion formation, whereby voters are located in a network and each one has one opinion among $𝑞$ possibilities. In ecological applications, voters become individuals and their opinion corresponds to the species they belong to (). An individual (voter) is chosen at random and is replaced by a copy of one of its neighbors (see Fig. ). This process has an absorbing state, because when a species disappears there is no way to introduce it again. In the long run, the system is populated by only one species. In order to overcome this problem, one can introduce the possibility of new species entering the system (; ), a model we will refer to as the MVM with speciation.

FIG. 8.

Microscopic moves of the multispecies voter model in a 2D lattice and with a nonregular network (). In the 2D structure, each site is occupied by one and only one individual, whose color represents the species it belongs to. At each time step, one random individual is replaced by a daughter of one of its neighbors with probability $1 - 𝜈$ (black lines). The probability that a speciation event occurs is $𝜈$ , wherein the individual is replaced by an individual of a new species (red lines). In the case of a nonregular network (right panel), the number of neighbors depends on the site considered.

Consider a lattice of $𝑑$ dimensions (for a typical ecological landscape $𝑑 = 2$ ), where $𝑎$ is the lattice spacing with exactly one individual at each site ( $𝑎 𝑑$ is the average volume occupied by one individual), and a total number of $𝑁$ individuals. At each time step, one individual chosen at random is removed and is replaced with a copy of one of its neighbors with a probability $1 - 𝜈$ , or with an individual of a species not already present in the system with a probability $𝜈$ . The case $𝜈 = 1$ is trivial, whereas the case $𝜈 = 0$ has an absorbing state, i.e., a state where a single species is present. The MVM is clearly a neutral model, and the microscopic dynamics is the same for all species. It is also a zero-sum process, because the total number of individuals $𝐽$ is conserved. The mean-field version of the voter model, where the neighbors of each node are all the other nodes, is governed by the ME , with the transition rates given by Eqs. and .

We want to write an equation for $𝐹 (𝑟)$ (the probability that two individuals picked at random at a distance $𝑟$ belong to the same species). Assuming translational invariance, one obtains

𝐹 (𝑟_) 𝑡 + 1 = 𝐹 (𝑟_) 𝑡 (1 - 2 𝑁) + 1 - 𝜈 𝑑 𝑁 𝑑 \sum 𝜇 = 1 [𝐹 (𝑟_+^𝑒 𝜇) 𝑡 + 𝐹 (𝑟_-^𝑒 𝜇) 𝑡],

(63)

where $𝐹 (𝑟_) 𝑡$ is the probability that two individuals separated by $𝑟_$ at time $𝑡$ belong to the same species. The solution is obtained with the boundary condition $𝐹 (0_) 𝑡 = 1$ .

The stationary solution of Eq. can be obtained from a Fourier series in the continuum limit, by taking the limit of $𝑎 \to 0$ , $𝑁 \to \infty$ , and $𝜈 \to 0$ , and constraining $𝛾 2 = 2 𝑑 𝜈 / 𝑎 2$ to be a constant. In this way one obtains () the following second order differential equation:

\nabla 2 𝐹 (r_) - 𝛾 2 𝐹 (r_) + 𝑐 𝛿 𝑑 (r_) = 0,

(64)

where $𝛿 𝑑$ is the $𝑑$ -dimensional Dirac delta, and whose solution is

𝐹 (𝑟) = 𝑐 𝛾 𝑑 - 2 (2 𝜋) 𝑑 / 2 (𝛾 𝑟) (2 - 𝑑) / 2 𝐾 (2 - 𝑑) / 2 (𝛾 𝑟),

(65)

where $𝐾 𝑧 (𝑟)$ is the modified Bessel function and $𝑐$ is a constant such that $\int 𝑟 < 𝑎 𝑑 𝑑 𝑟 𝐹 (𝑟_) = 1$ with $𝑎$ being the average distance between neighboring trees. By using the properties of $𝐾 𝑧 (𝑟)$ , one finds that in one dimension the solution is $𝐹 (𝑟) = 𝑐 𝜉 e x p (- 𝑟 / 𝜉) / 2$ , with $𝜉 = 𝛾 - 1$ . Interestingly, the exact solution in the (1D) discrete lattice is $𝐹 𝑥 = e x p (- 𝑥 / 𝜉)$ , where $𝑥 = 0, 𝑎, 2 𝑎, \dots$ and

𝜉 𝑎 = (l o g 1 - 𝜈 1 - \sqrt 𝜈 (2 - 𝜈)) - 1,

(66)

with $𝑎$ being the lattice spacing. Equation converges to $𝜉 = 𝛾 - 1$ in the continuum limit. The correlation length depends only on the speciation rate $𝜈$ and it grows as $1 / \sqrt 𝜈$ when $𝜈 \to 0$ , a limit that is distinct from the special case $𝜈 \equiv 0$ , in which at stationarity a single species is present across the whole finite lattice (monodominance). The $𝛽$ diversity obtained in Eq. was shown to be in good agreement with empirical data from different forest censa (). The prediction could be improved taking into account the Janzen-Connell effect ().

It is possible to obtain interesting results on the SAR via extensive numerical simulations. In the case of NT, one can take advantage of the coalescent approach (; ). Instead of simulating the stochastic dynamics directly, one can reconstruct the genealogy of the individuals in the sample area by regressing in time. The main value of this method is that there is no need to wait for any transient state to decay and, therefore, this approach is much faster than forward dynamics () as well as allowing infinite landscapes to be simulated.

The coalescent approach has been applied to the MVM with a different dispersal kernel (). Instead of a nearest-neighbor diffusion, once an individual is removed it is replaced by the offspring of another individual with a probability that depends on its distance from the individual removed. In an infinite landscape, the SAR shows the characteristic inverted-S shape. The model depends only on two parameters (up to a choice of the functional form of the dispersal kernel): the speciation rate $𝜈$ and the dispersal length $𝜉$ . The SAR scales as

⟨ 𝑆 (𝐴, 𝜉, 𝜈) ⟩ = 𝜉 𝑟 𝑆 (𝐴 𝜉 - 𝑟, 𝜈),

(67)

where the exponent $𝑟$ is independent of $𝜈$ and very close to 2 (as expected by dimensional analysis). The exponent $𝑧$ of the power-law regime can be calculated as the derivative of the curve evaluated at the inflection point of the SAR (in log-log scale), i.e., the minimum value of $𝑑 l o g (𝑆) / 𝑑 l o g (𝐴)$ . Additional simulations (), obtained with a wider spectrum of speciation rates, suggests a logarithmic relationship

𝑧 = 1 𝑞 + 𝑚 l o g (𝜈),

(68)

where $𝑞$ and $𝑚$ are two real parameters, which confirms the original prediction of . This inverse logarithmic trend seems very robust and it has also been observed in other spatial models (). Indeed, it was shown that power-law dispersal kernels better fit data than other short-ranged kernels (), whereby the estimation of the speciation rate corresponding to a given value of the exponent gives much smaller values.

An attempt to connect spatial and temporal patterns can be found (), where the persistence-time distribution was studied in MVM (see Sec. ). It was shown that the empirical persistence-time distribution is consistent with that predicted by MVM and

𝑝 (𝑡 | 𝐴) \sim 𝑡 - 𝛼 𝑒 - 𝑡 / 𝜏 (𝐴),

(69)

where $𝜏$ is the average time of persistence in an area $𝐴$ . The exponent $𝛼$ is universal and depends only on the dimension of the system, while the time scale $𝜏$ is a function of the sampled area. Empirically, the time scale $𝜏$ scales with the area as $𝐴 𝛽$ . In the MVM $𝜏 = 1 / 𝜈$ . Using this fact and that the rate of appearance of a new species is $𝜆 = 𝜈 𝑁 \sim 𝜈 𝐴$ , one can relate these quantities with the SAR. Indeed, the SAR is the product of the rate of appearance of new species and their average persistence time

𝑆 (𝐴) = 𝜆 ⟨ 𝑡 ⟩ \sim 𝐴 1 - 𝛽 (𝛼 - 1) .

(70)

Here we assumed that $⟨ 𝑡 ⟩ \sim 𝜏 2 - 𝛼$ to obtain a scaling relationship that connects the exponent $𝑧$ and the other exponents $𝑧 = 1 - 𝛽 (𝛼 - 1)$ . This allows spatial patterns to be connected with temporal ones in a manner consistent with the empirical patterns shown in Fig. .

C. Environmental fragmentation and habitat loss

Resources are not equally distributed and, even over small scales, their distribution in space is far from uniform. This spatial heterogeneity affects the distribution of individuals and species in space and has clear implications for the conservation of ecosystems. The space within which ecosystems are embedded is often fragmented and a species may be present in only part of the landscape. Moreover, different patches are not independent but they are rather connected via immigration.

In the absence of speciation, the voter model () predicts monodominance in dimension $𝑑 \leq 2$ . Spatial heterogeneity can be modeled as quenched disorder () and when two species are considered, it is postulated that different locations on a lattice prefer one species over the other. At each site $𝑖$ , a binary variable $𝜎 𝑖$ resides that takes values of $\pm 1$ , depending on which species is present. Spatial heterogeneity can be modeled as a quenched external field $𝜏 𝑖$ , which also takes values $\pm 1$ , and the dynamics are fully specified by the transition rates

𝑊 [𝜎 𝑖 \to - 𝜎 𝑖] = 1 - 𝜀 𝜏 𝑖 𝜎 𝑖 2 𝑧 \sum 𝑗 \in \partial 𝑖 (1 - 𝜎 𝑖 𝜎 𝑗),

(71)

where $\partial 𝑖$ is the set of nearest neighbors of $𝑖$ , and $𝑧$ is the size of this set. The quantity $𝜀$ measures the strength of the preference. The main result from this model () was that this randomness enhances the coexistence of species. Indeed, if $𝜀$ is larger than $𝜀 𝑐 = \sqrt 2 / (2 + 𝑁)$ , species coexist in any dimension in the limit of a large system. One can introduce the quantity $𝜙 = (1 / 𝑁) \sum 𝑖 𝜎 𝑖$ and if coexistence occurs $⟨ 𝜙 2 ⟩ < 1$ , whereas this quantity tends to 1 otherwise. It was shown that () for small $𝜀$

𝑃 (𝜙) \propto e x p {- (𝑁 / 2) [𝜀 2 / (1 - 𝜀 2)] 𝜙 2} (1 - 𝜀 2) (1 - 𝜙 2) + 2 𝜈,

(72)

where a small mutation rate [ $𝜈 ≪ 2 / (2 + 𝑁)$ ] was introduced to regularize the solution. If $𝜀 > 𝜀 𝑐$ , then the average $⟨ 𝜙 2 ⟩ < 1$ , i.e., coexistence is stable.

NT has also been applied to predict the extinction rate of species after habitat loss (). Habitat loss corresponds to a reduction of the area available and therefore to the total number of individuals. When this area is destroyed, the endemic species suddenly disappear and what follows is a delayed series of extinctions due to habitat loss. The community was modeled as a neutral assembly of species and a typical time scale of extinctions was obtained, along with its dependence on the number of species and the habitat destroyed. The predictions obtained in this way reproduced the available data of avifaunal extinctions well ().

V. Beyond Neutrality

A. Reconciling neutral and niche theory

The concept of a niche is central in classical ecology (; ; ). An ecological niche is “the requirements of a species for existence in a given environment and its impacts on that environment” () and it describes how an organism or a population responds to changes in resources, competitors, and predators. A possible mathematical realization of the concept of a niche is the Hutchinsonian niche, which involves a $𝑛$ -dimensional hypervolume, where the axes are environmental conditions or resources. A position in this space represents a set of behaviors and traits characterizing a species or a group of individuals. In niche theory much relevance is given to the specific traits of species and their interdependence. A central concept in niche theory is competitive exclusion, which states that two species cannot occupy the same niche, as two identical, yet distinct species, cannot coexist for an indefinite time.

The mathematical representation of an ecological community that includes niche aspects typically coincides with the Lotka-Volterra equations. In this case, the focus is on the properties of the fixed points (or other dynamical attractors) of these systems of equations and the typical problem that is analyzed is their stability, in relation to the parameters and the species present in the system.

The main difference between neutral and niche theory therefore depends on which mechanism plays the main role in shaping ecosystems (). Neutral theory assumes that random processes, such as dispersal, demographic stochasticity, speciation, and ecological drift, have a stronger impact on many of the observed patterns than niche differences. Niche theory assumes the opposite, that the quantitatively important processes are related to differences in species and their interdependence.

As we might expect in a real ecosystem both stochasticity and niche differences play a role, and it is natural to try to quantify how neutral behavior emerges from a niche model. In many cases, niche-based and neutral models yield compatible fits of biodiversity patterns (; ; ; ), and it is impossible to distinguish between the mechanisms by looking at those patterns. As pointed out by neutrality emerges when species have the same or very similar fitness.

Neutrality has often been proposed to emerge under some conditions from models considering niche differences (; ; ), and neutrality and niche theory were proposed to be the extremes of a more general model (; ). In both cases, a Lotka-Volterra equation is introduced to describe community dynamics, while population dynamics is modeled by also taking into account demographic stochasticity and immigration. By considering different values of parameters, one can move from a scenario where species differences matter a lot and the stable configuration is very close to the solution of the deterministic Lotka-Volterra equation to a scenario where demographic stochasticity is more important and the community behaves like a neutral community. A slightly different approach has also been considered (), analyzing a stochastic version of Lotka-Volterra dynamics and quantifying, when stochasticity is varied, the difference between the prediction of the neutral model with the full model. In this case, instead of a continuum of strategies, neutral and niche regimes are two macroscopic phases separated by a phase transition. Others have also incorporated neutral and non-neutral features within the same framework (; ): species with comparable sizes were considered functionally equivalent, thereby entailing a neutral dynamics; and parameters such as speciation, dispersal, and populations of organisms of distinct sizes were obtained from allometric scaling laws. With this approach, were able to explain why species richness reaches a maximum at intermediate body size.

A different approach was based on effectively considering niche theory as a model where per-capita death and birth rates are species dependent (). In this case, one might expect the difference between species abundances to reflect the differences between these parameters. More precisely, in a neutral scenario, all the species fluctuate around a given abundance value, while when niche characteristics play a role, each species fluctuates around its own distinct value of abundance. A third scenario was proposed (), wherein the per-capita birth and death rates do not have a monotonic effect on species abundance. The symmetry between species, due to the neutrality of the process, can yet be spontaneously broken. In this case, the stable states are not symmetric in terms of species abundance.

So far macroevolutionary patterns such as phylogenetic trees have not been a major focus in NT. However, there is evidence that NT is not able to satisfactorily capture phylogenetic diversity (). Multiple patterns of evolutionary history in bacterial communities seem to deviate from the predictions of NT (). However, recent models which incorporate either neutral population dynamics into the cladogenesis () or a mild selection into an otherwise neutral model () are able to slow down the diversification process as well as match phylogenies observed in nature. These represent promising developments, which can improve and broaden the spectrum of predictions of NT.

B. Emergent neutrality

Ecologists have highly criticized NT because of its unrealistic assumptions. The patterns that we have studied so far can in fact be explained without introducing species differences, and this has led some ecologists to oppose NT because it assumes that nature is actually governed by neutral processes, whereas it is not. Clearly, no one believes that nature is truly neutral. The patterns of community ecology are actually generated by a cocktail of processes, and it is both inappropriate and dangerous to consider processes in isolation from a macroecological pattern or empirical data set.

However, it is informative to study whether, how, and which ecological processes can drive a community toward a state in which demographic stochasticity and immigration play a crucial role in the face of strong species’ differences that are dictated by classical competitive exclusion. Such a state should allow similar species to emerge in the niche space with the ability to coexist for sufficient time. Indeed, when this problem was studied, it was shown that species can evolve into groups of relatively more similar species that coexist for very long times (). A large number of species were placed at random along a hypothetical niche axis, which represents a specific trait, assuming that interspecific competition can be calculated through niche overlap. Running a classical Lotka-Volterra competition model and studying evolution, groups of multiple species were evident that aggregated around similar values in the niche axis and they could coexist for many generations before the majority of them head toward an inexorable extinction. Eventually, only one species survives from each group, producing the expected pattern of single species equally spaced in the niche space.

In other words, the niche similarity of species prevents competitive exclusion from swiftly selecting the best competitor among a group of similar species, allowing their coexistence for very long times even though only the superior species will ultimately persist.

This model may be considered as one of the possible steps toward a reconciliation of niche and neutral theories. Species that are initially ecologically nonequivalent, and that therefore behave in a non-neutral fashion, are driven by community and evolutionary processes toward states in which the dynamics may well be better approximated by neutral models over appropriate spatial and temporal scales. More recently, further support for this approach came from showing that the model can produce multimodal RSAs (). Immigration may also be an important component in the neutral-like behavior of communities (; ). Parasitoids competing for a common species () have been used to show that clusters of species separated by gaps emerge along the niche axis, confirming—using a quite different approach—that processes exist that can lead community dynamics to be effectively neutral ().

The basic Lotka-Volterra model has been extended to investigate the possibility that some processes decrease the risk of competitive exclusion so that species lumps are not only transient, but ultimately permanent (). Density-dependent regulation was introduced that stabilizes the coexistence of species within a group. This approach unfortunately has a drawback that the mechanism introduces a discontinuity in the competition strength of the species (), which means that unmodeled species differences may be responsible for coexistence in the community. However, it has been shown that more realistic density-dependent terms, or even other mechanisms [e.g., migration ()], can eliminate this problem, making the approach more robust ().

C. Maximum entropy models

The maximum entropy principle is a useful method to obtain the least biased information from empirical measurements (). This powerful tool, borrowed from statistical mechanics, has a wide range of applications (), including ecology (; ). In its ecological application, the Max Ent principle is an inference method () used to evaluate the effective strength of interactions among species based on either species-abundance data () or simply the presence or absence of the species (). This methodology provides a way to systematically incorporate the most important species interactions into the development of a theory beyond the purely noninteracting case. In addition, the Max Ent principle was implemented as a method to predict biodiversity patterns across different spatial scales using only the information on local interactions (; ). Here we discuss the extension of Max Ent to the study of spatial biodiversity patterns. Consider an ecosystem in which $𝑆$ species (belonging to the same trophic level, see Table ) live within a given area $𝐴$ , divided into $𝑁$ adjacent sites of equal area. Let us assume that there are empirical records of the species contained within each site. From these data one can calculate, for example, the average presence of any species in the ecosystem and the co-occurrence of any pair of species in neighboring sites. When applying Max Ent, we can consider these mean occurrences and co-occurrences as given constraints. We introduce the binary random variable $𝜎 𝛼 𝑖$ , which records the occurrence of species $𝛼 = 1, \dots, 𝑆$ at each site $𝑖 = 1, \dots, 𝑁$ . If species $𝛼$ is present in plot $𝑖$ , then $𝜎 𝛼 𝑖 = 1$ , otherwise $𝜎 𝛼 𝑖 = 0$ . Therefore the “state” of any species $𝛼$ can be characterized by the random vector $𝝈 𝛼$ . There is empirical evidence that species belonging to the same trophic level interact weakly (; ; ) and therefore, in a first approximation, one may assume that species occur in a given geographical location independently of one another. Because of such independence, the probability of finding the system in the configuration $𝝈 = (𝝈 1, 𝝈 2, . ., 𝝈 𝑆)$ is thus $𝑃 (𝝈) = \prod 𝑆 𝛼 = 1 𝑝 𝛼 (𝝈 𝛼)$ , where $𝑝 𝛼$ gives the probability distribution of finding a species $𝛼$ in the configuration $𝝈 𝛼$ . In order to build the Max Ent model, we maximize Shannon’s entropy $𝐻 = - \sum 𝝈 1 \dots \sum 𝝈 𝑆 𝑃 (𝝈) l n 𝑃 (𝝈) = - \sum 𝛼 \sum 𝝈 𝛼 𝑝 𝛼 (𝝈 𝛼) l n 𝑝 𝛼 (𝝈 𝛼)$ , imposing the average occurrence constraint, i.e., $⟨ 𝑀 ⟩ = \sum 𝛼 \sum 𝝈 𝛼 𝑝 𝛼 (𝝈 𝛼) 𝑀 𝛼 (𝝈 𝛼)$ with $𝑀 𝛼 (𝝈 𝛼) = \sum 𝑖 𝜎 𝛼 𝑖$ , and the average co-occurrence constraint, i.e., $⟨ 𝐸 ⟩ = \sum 𝛼 \sum 𝝈 𝛼 𝑝 𝛼 (𝝈 𝛼) 𝐸 𝛼 (𝝈 𝛼)$ with $𝐸 𝛼 (𝝈 𝛼) = \sum (𝑖, 𝑗) 𝜎 𝛼 𝑖 𝜎 𝛼 𝑗$ , where $(𝑖, 𝑗)$ indicates two nearest neighbor locations. Both $⟨ 𝑀 ⟩$ and $⟨ 𝐸 ⟩$ are meant to match the corresponding empirical averages, as calculated from the empirical records of species contained in the region. Maximization () provides an expression for the probability $𝑝 𝛼$ of finding a species $𝛼$ in the configuration $𝝈 𝛼$ . Thus,

𝑝 𝛼 (𝝈 𝛼) = 𝑒 𝐽 𝛼 𝐸 𝛼 (𝝈 𝛼) + ℎ 𝛼 𝑀 𝛼 (𝝈 𝛼) 𝑍 𝛼,

(73)

where $𝑍 𝛼 = 𝑍 𝛼 (ℎ 𝛼, 𝐽 𝛼)$ is the partition function. From Eq. one can characterize the spatial biodiversity patterns of the ecosystem, i.e., calculate the SAR or the EAR (). This type of approach may well be suited to infer biodiversity properties of communities over larger spatial scales by upscaling the model results at local scales. However, only a few studies have tackled this problem (see Sec. ).

VI. Open Problems and Future Challenges

In this section, we consider how one might advance theory in terms of adding essential details without making the system unnecessarily complex or having an explosion of species-specific parameters. What is the relative importance of species traits, their interactions, and spatial and environmental effects? Throughout this review, we have treated the total number of species and the total number of individuals as parameters which are provided as input. How would one determine and predict these parameters? A mechanistic explanation is needed for why one region may be more biodiverse than another, what sustains biodiversity, and how evolutionary pressures sculpt ecological communities.

A. Analytical spatially explicit models

At the moment, we do not have a spatially explicit theory that can analytically predict the most important ecological patterns. This would be important, because it would allow us to understand what drives biodiversity across spatial and temporal scales. Indeed, it is likely that biological processes are not equally important across scales. As evident in particle physics, we might eventually find that ecosystems will need to be described by different effective theories according to the range of scales in which we are interested. In a pioneering paper, which received great interest, a ME was proposed that seemed to permit analytical calculations of the SAR (). Indeed, an analytical expression for the SAR was obtained. However, the promised solution was neither correct, nor was it an approximation of the actual SAR (). The main technical difficulty preventing an analytical solution was that the model was inherently out of equilibrium. Detailed balance, which is the condition that makes the calculation of stationary probabilities possible, did not hold (). The time is ripe for explorations of this kind. The challenge is to develop new and powerful techniques in nonequilibrium statistical mechanics. Others have applied NT to study biodiversity patterns on river networks. By exploiting the network connectivity of such dendritic landscapes (), they were able to explain several empirical large-scale spatial biodiversity patterns (, ; ) and show how changes in riverine ecosystems may impact the spread of species and local species richness (; ).

B. Linking different macroecological patterns

Over the last few decades, ecologists have come to appreciate the importance of spatial patterns and processes, and the explicit introduction of space has the potential to revolutionize what we know about natural populations and communities (; ). It has become apparent that key ecological patterns, such as SAR, RSA, and spatial patterns of species distributions and turnover, are intimately intertwined and scale dependent (see Fig. ). However, despite a plethora of models that address spatial patterns, only a few practical methods have been proposed to link them across different scales. Specifically, spatial approaches (; ) lack the needed analytical machinery, whereas most theoretical approaches are not spatially explicit or sufficiently flexible (; ). Even neutral theory was conceived to reflect the idealized behavior of natural systems at equilibrium, rather than to reflect nonpristine landscapes produced by environmental change or management (). Therefore, a general methodology is required to predict and link these ecological patterns across scales that is sufficiently robust and flexible to allow its application to a range of natural or managed systems. One possible way to tackle this challenging problem is through a theoretical framework inspired by ideas coming from phenomenological renormalization (). The fundamental assumption at the core of this theoretical setting is that the functional form of the RSA remains the same across all spatial scales, even though the parameters of the curve are likely to vary. Because of this assumption, the spatial dependence of the abundance distribution can be obtained by making the RSA parameters suitable functions of scale. Together with the functional shape of the RSA, the other model input is the spatial PCF, which describes the correlation in species’ abundances between pairs of samples as a function of the distance between them (; ). If populations were randomly distributed in space, distinct communities would on average share the same fraction of species regardless of their spatial separation, and therefore the PCF would not depend on distance. In contrast, in highly aggregated communities correlations in abundance would fall off steeply with increasing distance. The PCF not only measures the rate of turnover in species composition but it also reflects the variation of population clustering across scales, given that the variance in species abundances at any particular scale can be calculated directly from the PCF (). Therefore, the PCF is related to the spatial species-abundance distribution. Thus, the PCF can link the effects of aggregation, similarity decay, species richness, and species abundances across scales. Building on the intrinsic relationship among these patterns, while accounting for spatial correlations with fidelity, is critical for predicting the biodiversity profiles across scales when information on a limited number of fine-scale scattered samples is available (; ). The explosion of publicly available large-scale biodiversity data for paradigmatic ecosystems such as the Amazonian forest () and ocean plankton () makes this problem one of the most exciting scientific challenges in this field.

C. Environmental noise

Most neutral (symmetric) models are solely governed by the underlying stochastic birth-death process and do not take into account any environmental stochasticity (). Recently it was shown that environmentally induced variations of the demographic rates dominate the long-term dynamics and have an important impact on some dynamic properties [such as age-size relationships and species extinction time ()], while not affecting the already good accord of neutral models with ecological static patterns (such as the RSA) (). However, analytical results on neutral (symmetric) models for ecological communities subject to correlated environmental noise are mainly missing [although some recent results have been obtained for white noise (; )], incentivizing physicists to explore these paths.

D. Multitrophic ecosystems: Ecological networks

Although we are progressing in understanding the suitability and limits of noninteracting and nonspatial models, most neutral models still assume that species interact randomly with each other. However, a network approach to modeling ecological systems provides a powerful representation of the interactions among species (; ; ). Ecological networks may be viewed as a set of different species (nodes) and connections or links (edges) that represent interspecific interactions (e.g., competition, predation, parasitism, and mutualism). The architecture of ecological interaction networks has become a bubbling area of research, and it seems to be a critical feature in shaping and regulating community dynamics and structure diversity patterns (; ). An important step relevant for multitrophic systems will be to obtain a general framework within which a network of preferences or disfavors modifies the birth and death rates of different species and can be superposed on neutral models (like the voter model presented in this review). Recent studies of ecological networks have considered the exciting task of anticipating critical transitions in such systems and to design structures that are less vulnerable to collapse (; ). Connecting stochastic quasineutral models, ecological networks, and critical transitions within a unified theoretical framework is an important challenge as it will enhance our capacity to understand and thus manage the crucial interplay between ecological dynamics and species interactions.

VII. Conclusions

In this paper, we have attempted to describe some of the theoretical frameworks that can be used to understand key issues related to biodiversity and that will serve to address important questions. These frameworks are necessarily elementary and incomplete, yet they have the advantage of being tractable and related to the central issues. Unlike standard approaches to more traditional physics, here the Hamiltonian function or the interactions among the components are completely unknown and even identifying the state variables may sometimes be a nontrivial task. We have a fleeting picture of what an ecosystem is and it is not necessarily in equilibrium. We have little knowledge of the myriad degrees of freedom and their interactions. The real challenge is to discern the most essential degrees of freedom and to develop a framework to understand and predict the emergent ecological behavior.

The basic message of this review is that to resolve these challenging problems, ideas and techniques must be recruited from different disciplines. We are still at the beginning of this adventure. Moving forward is not only important but it is also urgent. The pressures of habitat destruction, pollution, and climate change are having highly undesirable consequences on the health of ecological communities. To address practical issues related to conservation biology, we need models that can be used across scales in order to extrapolate information on biodiversity from accessible regions to inaccessible yet important scales. There is plenty of room for ideas that matter, and community ecology can greatly benefit from the contribution of other disciplines, including physics.

Acknowledgments

We have greatly benefited from discussions with Stephen Hubbell, Stefano Allesina, David Alonso, Enrico Bertuzzo, Ryan Chisholm, Stephen J. Cornell, Paolo D’Odorico, William F. Fagan, Andrea Giometto, William E. Kunin, Egbert Leigh, Alan J. McKane, Miguel A. Muñoz, Andrew Noble, James P. O’Dwyer, Andrea Rinaldo, Ignacio Rodriguez-Iturbe, James L. Rosindell, and Nadav Shnerb. S. S. acknowledges the University of Padova Physics and Astronomy Department Senior Grant No. 129/2013 Prot. 1634.

Appendices

Appendix A: Density Dependence and Assembling of Local Communities

In this Appendix, we present an alternative method to introduce density dependence and how ecosystems emerge by assembling local communities.

First, one can set $𝑏 𝑘 (𝑛) = 𝑏 (𝑛 + ϒ 𝑘)$ (for the $𝑘$ th species) and $𝑑 (𝑛) = 𝑑 𝑛$ , where $ϒ 𝑘$ incorporates the effects of both intraspecific interactions, such as those giving rise to density dependence, and the immigration occurring from a metacommunity. This approach can be applied to two distinct ecosystems: coral reefs and tropical forests.

The steady-state solution of the ME for $𝑃 𝑘 (𝑛)$ yields a negative binomial distribution:

𝑃 𝑘 (𝑛) = (1 - 𝑟 𝑘) ϒ 𝑘 Γ (𝑛 + ϒ 𝑘) Γ (ϒ 𝑘) 𝑟 𝑛 𝑘 𝑛!

(A1)

with a mean $⟨ 𝑛 𝑘 ⟩ = 𝑟 𝑘 ϒ 𝑘 / (1 - 𝑟 𝑘)$ , where $Γ$ is the gamma function.

The number of species containing $𝑛$ individuals is given by $𝜙 𝑛 = \sum 𝑆 𝑘 = 1 𝐼 𝑛, 𝑘$ , where $𝐼 𝑛, 𝑘$ is a random variable that is 1 with a probability $𝑃 𝑘 (𝑛)$ and 0 with a probability $1 - 𝑃 𝑘 (𝑛)$ . Thus, the RSA is given by

⟨ 𝜙 𝑛 ⟩ = 𝑆 \sum 𝑘 = 1 𝐼 𝑛, 𝑘 = 𝑆 \sum 𝑘 = 1 𝑃 𝑘 (𝑛) = 𝜃 𝑟 𝑛 𝑛! Γ (𝑛 + ϒ),

(A2)

where $𝜃 = 𝑆 / [(1 - 𝑟) - ϒ - 1] Γ [ϒ]$ is the biodiversity parameter (), and we dropped the $𝑘$ dependence because of the symmetric hypothesis. For a small $ϒ$ , the RSA for the communities resembles the Fisher log series and it does not have an interior mode.

Note that a nontrivial $𝑘$ dependence might arise even under the neutral hypothesis. For example, one can set $ϒ 𝑘 = ˜ 𝑚 𝑝 𝑘$ , where $˜ 𝑚$ is a measure of the immigration rate from the metacommunity, in units of the birth rate $𝑏$ and $𝑝 𝑘$ is the fraction of individuals in the surrounding metacommunity belonging to the $𝑘$ th species. As a result, one can obtain the following RSA for the community of tropical forests:

⟨ 𝜙 𝑛 ⟩ = 𝜃 𝑥 𝑛 𝑛! \int \infty 0 Γ (𝑦 + 𝑛) Γ (𝑦 + 1) 𝑒 - 𝜔 𝑦 𝑑 𝑦 \equiv 𝜃 𝑥 𝑛 𝑛! 𝑓 (𝑛, 𝜔),

(A3)

where $𝜔 = 𝜃 / ˜ 𝑚 - l n (1 - 𝑥)$ . This approach provides a virtually indistinguishable fit to the empirical data as

⟨ 𝜙 𝑛 ⟩ = 𝜃 𝑥 𝑛 𝑛 + 𝑐

considered earlier with the advantage of having ecologically meaningful parameters.

The average number of species observed in the local community is

𝑆 o b s = ⟨ 𝑆 ⟩ = 𝑆 - ⟨ 𝜙 0 ⟩ = 𝑆 \sum 𝑘 = 1 (1 - 𝑟) - ϒ = 𝑆 [1 - (1 - 𝑟) ϒ] .

(A4)

If the sample considered has $𝐽 𝐿$ individuals, and thus the community dynamics obeys the zero-sum rule (i.e., it has a fixed total population), then the multivariate probability distribution is

𝑃 (𝐧 | 𝐽 𝐿) = 𝒩 𝑆 \prod 𝑘 = 1 𝑃 𝑘 (𝑛) 𝛿 (𝐽 𝐿 - 𝑛 1 - 𝑛 2 - \dots - 𝑛 𝑆) = (𝐽 𝐿 + \sum 𝑘 ϒ 𝑘 - 1 𝐽 𝐿) - 1 𝑆 \prod 𝑘 = 1 (𝑛 𝑘 + ϒ 𝑘 - 1 𝑛 𝑘) \times 𝛿 (𝐽 𝐿 - 𝑛 1 - 𝑛 2 - \dots - 𝑛 𝑆)

(A5)

that is a compound multinomial Dirichlet distribution where $𝒩$ is the normalization constant.

Metacommunity composition: Now let us gradually assemble the metacommunity of coral reefs by considering it as an assemblage of local communities. Let us start by considering the joint RSA of two local communities $𝐴$ and $𝐵$ with $𝑛 𝐴$ and $𝑛 𝐵$ individuals, respectively. The probability that the species has $𝑛$ individuals in the metacommunity formed by $𝐴$ and $𝐵$ is ()

𝑃 (𝑛 = 𝑛 𝐴 + 𝑛 𝐵) = \sum 𝑛 𝐴 + 𝑛 𝐵 = 𝑛 𝑃 (𝑛 𝐴) 𝑃 (𝑛 𝐵) \propto 𝑟 𝑛 𝑛! Γ (𝑛 + 2 ϒ),

(A6)

where the actual spatial locations of the local reef communities have been neglected (all local communities are well mixed in the metacommunity, i.e., a mean-field approximation). The elegant result of the $ϒ$ ’s adding to each other follows ecologically from their interpretation as immigration rates. For the metacommunity, we can introduce speciation with a rate $𝜈 ≪ 1$ as we have seen it has a crucial role when $𝑛 = 0$ (under the assumption of neutrality, the species label of the new species is of no consequence), i.e., $ϒ 𝑘 = 𝜈$ and $𝑃 (𝑛 𝑘) = 𝜈 𝑟 𝑛 / 𝑛 + 𝑂 (𝜈 2)$ .

Extending the calculation of the joint RSA distribution to more and more local communities, it can be shown that the RSA of the metacommunity is characterized by an effective immigration parameter $𝐿 ϒ$ , where $𝐿$ is the total number of local communities comprising the metacommunity, and it becomes log normal like if $𝐿 ≫ 1$ , in agreement with the available data ().

A Fisher log series is observed in two limiting cases—in the metacommunity in which there are no immigration events and in the very small local community that has a high immigration rate from the metacommunity characterized by a Fisher log-series RSA.

Appendix B: Persistence or Lifetime Distributions

In this Appendix, we present the derivation for the persistence (or lifetime) distribution that leads to the empirical choice made in Eq. . The master equation we want to solve is

\partial 𝑃 (𝑛, 𝑡) \partial 𝑡 = 𝑏 (𝑛 - 1) 𝑃 (𝑛 - 1, 𝑡) + 𝑑 (𝑛 + 1) 𝑃 (𝑛 + 1, 𝑡) - [𝑏 (𝑛) + 𝑑 (𝑛)] 𝑃 (𝑛, 𝑡),

(B1)

where the birth rate is $𝑏 (𝑛) = (1 - 𝜈) 𝑛 / 𝐽 𝑀$ with $𝑏 (- 1) = 0$ , and the death rate is $𝑑 (𝑛) = 𝑛 / 𝐽 𝑀$ and $𝑛 \geq 0$ . We can redefine the time scale $𝑡 \to 𝐽 𝑀 𝑡$ so that the factor $1 / 𝐽 𝑀$ disappears from the birth and death rates. As the initial condition we choose $𝑃 (𝑛, 𝑡 = 0) = 𝛿 𝑛, 𝑛 0$ with the initial population as $𝑛 0$ . Our goal is to calculate the survival probability defined as

𝒫 (𝑡 | 𝜈) = \sum 𝑛 \geq 1 𝑃 (𝑛, 𝑡) = 1 - 𝑃 (𝑛 = 0, 𝑡) .

(B2)

Thus, we introduce the generating function

𝐺 (𝑧, 𝑡) = \sum 𝑛 \geq 0 𝑃 (𝑛, 𝑡) 𝑧 𝑛

(B3)

whose radius of convergence is $\geq 1$ . Note that $𝐺 (𝑧 = 1, 𝑡) = 1$ and $𝐺 (𝑧 = 0, 𝑡) = 𝑃 (𝑛 = 0, 𝑡) = 1 - 𝒫 (𝑡 | 𝜈)$ . Using Eq. the time evolution of the generating function is derived immediately

\partial 𝐺 (𝑧, 𝑡) \partial 𝑡 = [(1 - 𝜈) 𝑧 - 1] (𝑧 - 1) \partial 𝐺 (𝑧, 𝑡) \partial 𝑧

(B4)

with the initial condition $𝐺 (𝑧, 𝑡 = 0) = 𝑧 𝑛 0$ . The previous equation is a linear partial differential equation and it can be solved by standard methods. One introduces a function $𝑍 (𝜏)$ satisfying the time evolution equation

𝑑 𝑍 (𝜏) 𝑑 𝜏 = - [(1 - 𝜈) 𝑍 (𝜏) - 1] [𝑍 (𝜏) - 1],

(B5)

with a “final” condition $𝑍 (𝜏 = 𝑡) = 𝑧$ . Using Eq. , one is led to

𝑑 𝐺 (𝑍 (𝜏), 𝑡) 𝑑 𝜏 = 0, \forall 𝜏,

(B6)

implying that

𝐺 (𝑧, 𝑡) = 𝐺 (𝑍 (𝑡), 𝑡) = 𝐺 (𝑍 (0), 0) = 𝑍 (0) 𝑛 0 .

(B7)

The solution of Eq. , with the chosen final condition gives

𝑍 (𝜏 = 0) = 1 - 𝐴 (𝑧, 𝑡) 1 - (1 - 𝜈) 𝐴 (𝑧, 𝑡) with 𝐴 (𝑧, 𝑡) = 1 - 𝑧 1 - (1 - 𝜈) 𝑧 𝑒 - 𝜈 𝑡

(B8)

and thus

𝐺 (𝑧, 𝑡) = (1 - 𝐴 (𝑧, 𝑡) 1 - (1 - 𝜈) 𝐴 (𝑧, 𝑡)) 𝑛 0 .

(B9)

Finally, we get the survival probability

𝒫 (𝑡 | 𝜈) = 1 - 𝐺 (𝑧 = 0, 𝑡) = 1 - [1 + 𝜈 (𝑒 𝜈 𝑡 - 1) - 1] - 𝑛 0 .

(B10)

In the scaling, if we consider the limit of fixed $𝜈 𝑡$ as $𝑡$ becomes large, we get the scaling form

𝒫 (𝑡 | 𝜈) = 1 𝑡 𝐹 (𝜈 𝑡),

(B11)

with

𝐹 (𝑥) = 𝑛 0 𝑥 𝑒 𝑥 - 1 .

(B12)

The lifetime distribution is simply given by

𝑝 𝜏 (𝑡) = - 𝑑 𝒫 (𝑡 | 𝜈) / 𝑑 𝑡 = 1 𝑡 2 𝑓 (𝜈 𝑡),

(B13)

where the last equality holds in the scaling limit as

𝑓 (𝑥) = 𝑛 0 𝑒 𝑥 (𝑥 𝑒 𝑥 - 1) 2,

(B14)

which for $𝑥 \to 0$ tends to a nonzero constant. As $𝑥 \to \infty$ , $𝑓 (𝑥)$ decays exponentially, leading to the empirical scaling form of the kind given by Eq. .

Appendix C: The Species Turnover Distribution with Absorbing Boundary Conditions

In this Appendix, we present another temporal pattern predicted by the model defined by Eq. or . For an ecosystem at or near stationarity, the model can provide the exact expression for the STD when using the time-dependent absorbing solution of Eq. : the so-called species turnover distribution with absorbing boundary conditions. Unlike the reflecting species turnover distribution, this distribution corresponds to a new kind of measure that only accounts for the species present at the initial time, and it does not take into account any new species introduced by immigration or speciation or any old species that reappear after their apparent extinction until $𝑡 > 0$ . This amounts to the selection of a particular sample and the study of the temporal behavior of those selected individuals. The species turnover distribution with absorbing boundary conditions $𝒫 𝑎 𝑏 𝑠 (𝜆, 𝑡)$ can be obtained through Eq. where the conditional probability must now be $𝑝 𝑎 (𝑥, 𝑡 | 𝑥 0, 0)$ defined as in Eq. . The final expression is

𝒫 𝑎 𝑏 𝑠 (𝜆, 𝑡) = s i n (𝜋 𝑏 / 𝐷) 𝜋 (1 - 𝑏 / 𝐷) 𝑒 - 𝑡 / 𝜏 (𝑒 𝑡 / 𝜏 - 1) 𝑏 / 𝐷 (𝜆 + 1) 2 \times 2 𝐹 1 (1, 32, 2 - 𝑏 / 𝐷; 4 𝜆 𝑒 - 𝑡 / 𝜏 (𝜆 + 1) 2),

(C1)

where $𝑏 / 𝐷 < 1$ . Note that $𝒫 𝑎 𝑏 𝑠 (𝜆, 𝑡)$ decays exponentially to zero at a rate $(1 - 𝛽) / 𝜏$ , regardless of $𝜆$ .

It is noteworthy that the STD’s with absorbing and reflecting boundaries are indistinguishable whenever $𝑡 ≪ 𝜏$ and $𝜆$ are not too small (for $𝑏 / 𝐷 < 1$ ). Since the BCI data are sampled over relatively short times intervals (at most 10 years), the distributions are almost the same within the interval $1 / 2 < 𝜆 < 3 / 2$ .