Abstract

Ecosystems display a complex spatial organization. Ecologists have
long tried to characterize them by looking at how different measures
of biodiversity change across spatial scales. Ecological neutral
theory has provided simple predictions accounting for general
empirical patterns in communities of competing species. However,
while neutral theory in well-mixed ecosystems is mathematically well
understood, spatial models still present several open problems,
limiting the quantitative understanding of spatial biodiversity. In
this review, we discuss the state of the art in spatial neutral
theory. We emphasize the connection between spatial ecological
models and the physics of non-equilibrium phase transitions and
how concepts developed in statistical physics translate in
population dynamics, and vice versa. We focus on non-trivial scaling
laws arising at the critical dimension D=2 of spatial neutral
models, and their relevance for biological populations inhabiting
two-dimensional environments. We conclude by discussing models
incorporating non-neutral effects in the form of spatial and
temporal disorder, and analyze how their predictions deviate from
those of purely neutral theories.

Community ecology aims at shedding light on how competing species
assemble and coexist in their habitats Cody and Diamond (1975). This
has proven to be a formidable challenge. A main reason is that
ecological dynamics span a wide range of spatial and temporal scales,
from those typical of individuals to those characterizing large
populations or communities. Ecologists have
empirically characterized biodiversity at the different spatial
scales; for example, counting the average number of species hosted in
a given area – species area relationship (SAR)
Preston (1960); Rosenzweig (1995)–, or the distribution of their
abundances – species abundance distribution (SAD)
MacArthur (1960); Tokeshi (1993). Often, the ecological
forces determining these patterns act at a given spatio-temporal
scale but can affect others as well. The inverse problem,
i.e. linking observed patterns with the causes originating
them at different scales, is arguably the central problem in ecology
Levin (1992).

This kind of problem sounds familiar to experts in statistical
physics, where large-scale emergent behavior results from interactions
among simple local units. Tools of statistical physics are indeed very
useful to make progress on the aforementioned crucial issues in
ecology. In particular, a natural approach to such complex problems is
to radically simplify them. To this aim, we consider ecosystems made
up of competing non-motile species, such as trees, or having a
motility range much smaller than the typical linear size of the
population, such as communities of
microorganisms. Further possible simplifications are
that all emergent phenomena originate at the single-individual scale
and, more drastically, that differences among individuals, possibly
belonging to different species, can be neglected. These assumptions
constitute the basis of the ecological neutral theory proposed
by Hubbell Hubbell (2001).

Ecological neutral theory Hubbell (2001) was built upon
theoretical ideas of Kimura’s neutral theory of population genetics
Kimura (1983). Both theories underscore the
role of stochastic demographic fluctuations in determining the fate of
populations and completely neglect deterministic effects stemming
from fitness differences. The assumption of ecological neutrality has elicited
heated controversies, as it hinted that classical ecological
concepts, such as niches, might play a marginal role in
structuring communities of competing species. Despite these
contentions, neutral theory had a considerable impact on
ecological thinking, owing to its ability to quantitatively predict
non-trivial patterns of biodiversity with simple models
characterized by very few
adjustable parameters
Alonso et al. (2006); Rosindell et al. (2011); Azaele et al. (2016).

The focus of this review is on spatially-explicit neutral and
near-neutral population models. Explicitly describing space is crucial
to address the fundamental ecological questions sketched at the
beginning of the introduction. However, spatially-explicit models –
that are often
variants of familiar models in non-equilibrium statistical physics
Durrett and Levin (1994) – are
still poorly understood, especially if compared with their well-mixed
counterparts Etienne and Rosindell (2011). One of the most studied
neutral model is the voter model with speciation, or
multi-species voter model
Durrett and Levin (1996); Rosindell and Cornell (2007); Pigolotti and Cencini (2009), which
generalizes the more common two-species voter model
Liggett (1985). The stepping-stone model
Kimura (1953); Kimura and Weiss (1964); Korolev et al. (2010); Cencini et al. (2012)
and the contact process Liggett (1985); Marro and Dickman (1999); Hinrichsen (2000)
are other examples of spatial models that have been studied in
both the physics and population biology literature. We shall
discuss how these analogies can be used to advance our understanding
of spatial ecology and the main open problems. This review heavily
relies on extensive numerical computations of lattice models based on
previous works by the authors. This might have biased the choice of
some topics and we apologize if some relevant works are not properly
discussed.

The review is organized as follows. In Sect. II we introduce
the multispecies voter model on a lattice and its dual representation
in terms of coalescing random walkers. We then discuss its predictions
of macroecological patterns: the SAR, and the SAD. For the latter, we
compare two recent analytical approaches
Zillio et al. (2008); Shem-Tov et al. (2017); Danino et al. (2016a) with novel computational
results. We mainly discuss the two-dimensional case
due to its ecological relevance, but also briefly present the
one-dimensional case for comparison. We conclude the section by
presenting new results on an important dynamical property: the distribution of
species persistence-times. In Sect. III we
discuss other neutral models, where, at variance with
the voter model, lattice sites are not necessarily occupied by exactly
one individual at all times. In particular, we consider the stepping
stone model, where each lattice site hosts a local community of
individuals. This generalization is relevant for
modeling microorganisms and their macroecological patterns. We then
consider a multispecies variant of the contact process, where lattice
sites can be either empty of occupied by single individual. In
Section IV we introduce non-neutral effects on a
simplified two-species competition model, where adjusting a single
parameter one can tune the departure from neutrality, here modeled as
a specific habitat preference. Physically, this habitat preference can
be thought as a form of quenched disorder. We discuss how this
disorder generically favors species coexistence using the language of
statistical mechanics, and also discuss other forms of disorder such
as temporal heterogeneity. Finally, Sect. V is
devoted to perspectives and conclusions.

ii.1 Description of the model

A paradigmatic example of spatial neutral model is the voter
model with speciation, Durrett and Levin (1996), which is is a
multispecies generalization of the voter model
Liggett (1985). The latter is a widely studied model
that has been applied in diverse contexts, from
population genetics to spatial conflicts Clifford and Sudbury (1973), spreading of
epidemic diseases Pinto and Munoz (2011), opinion dynamics
Castellano et al. (2009) and linguistics Croft (2002).

The voter model with speciation is defined on a lattice, where each
site hosts one individual belonging to some species. At each discrete
time step, a lattice site is chosen at random and the residing
individual is removed (death event). Then,
as illustrated Fig. 1, the dead individual is replaced:

With probability ν, by an
individual of a new species not present in the system (speciation
event). Notice that, because of speciation, the total number of
species is not fixed. In population genetics, this type of event is
interpreted as a mutation within the same species
Moran (1958); Kimura and Weiss (1964).

With complementary probability (1−ν), by a new individual of
an existing species (reproduction event). In this case, the newborn
belongs to the same species of a parent individual
chosen at random in the neighborhood of the vacant site. In the
simplest case, the nearest-neighbors (NN) are chosen with uniform
probability. More generally, the parent individual is selected
according to a probability distribution P(→r) (the dispersal kernel) over the neighbors within a
distance →r.

Figure 1: Examples of transitions in the 2D voter model
with speciation.

Most of this section will be devoted to the ecologically relevant case
where the system is a two-dimensional (2D) square lattice,
although we will briefly present some results in 1D
for comparison.

ii.2 Duality

The voter model with speciation is dual to a system of
coalescing random walkers with an annihilation rate
Holley and Liggett (1975); Bramson et al. (1996); Durrett and Levin (1996). In this
context, “duality” means that each trajectory of one system can be
mapped in one of the other system having equal probability
Holley and Liggett (1975). The dual process is constructed as follows. We start by placing on each lattice
site a random walker. The dynamic of the dual process proceeds
backward in time. At each discrete (backward) time step, with
probability 1−ν, a randomly chosen walker is moved to a new site,
where the dispersal kernel P(→r) here plays the role of the
distribution of possible displacements. If the site is occupied, the
two walkers coalesce, i.e. one of the two is removed keeping trace of
the coalescing partner. With complementary probability ν a
randomly chosen random walker is annihilated, i.e. removed from the
system. This event corresponds to a speciation event in the forward
dynamics. The whole tree of coalescing random walkers, before
annihilation, represents the entire genealogical tree of a species up
to the speciation event that originated it.

The standard forward in time evolution of the voter-model with
speciation and its dual dynamics are sketched, for
the one-dimensional case, in Fig. 2a and 2b,
respectively.

Figure 2: a) Example of space-time dynamics of the 1D voter model
with speciation. b) Corresponding dual dynamics: coalescing and
annihilating random walkers. c) Snapshot of a configuration of
the 2D voter model simulated with the dual dynamics, with
ν=510−7 and nearest-neighbor (NN) dispersal. d) Same as c)
but with a longer dispersal range (uniformly distributed in a square of
side K) with K=7. Each color labels a different
species.

Duality is a very useful property to understand the physics of the
voter model. For example, it immediately stems from duality that the
ν→0 limit is fundamentally different in D≤2 and
D>2. As a matter of fact, in D≤2 the random walk is recurrent,
meaning that the probability of two randomly chosen individuals to
belong to the same species approaches one as ν→0. In
other words, in the absence of speciation, one has
monodominance of one species in the long term. The same property does
not hold in D>2, where random walkers are not recurrent and,
in an infinite system, multiple species coexist on the long term even in the
limit ν→0. Interestingly, the ecologically most relevant case,
D=2, is the critical dimension of this model. We shall see that this
fact is a source of non-trivial behaviors of
ecologically relevant quantities.

Duality is also an extremely powerful tool for computational analyses
Rosindell and Cornell (2007); Pigolotti and Cencini (2009). If one is interested in the
static, long-term, properties of the voter model with speciation, it
is numerically much more efficient to simulate the dual dynamics than
the forward one. In a dual simulation, after all walkers coalesced or
were annihilated, species can be assigned to the start site of each
walker, obtaining a stationary configuration of the voter
model. Beside computational speed, this approach has also the
advantage of eliminating finite-size effects induced by the boundary
conditions, as the coalescing random walkers can be simulated in a
virtually infinite system. For illustrative purposes, in
Fig. 2c and 2d we show two configurations of
the 2D voter model obtained with the dual dynamics for two different
dispersal kernels.

ii.3 β−diversity

The first ecological pattern we consider is the
β-diversity, which is a measure of how the species composition
in an ecosystem varies with the distance. We define the
β-diversity as the probability F(→r), that two
randomly chosen individuals at a distance →r are conspecific,
i.e. belong to the same species. We remark that,
although this is the natural definition in this context, other
definitions have been used
in the ecological literature Tuomisto (2010). Mathematically, F(→r)
can be expressed in terms of the two-point correlation function
Gsi,sj(→r)=⟨nsi(→x)nsj(→x+→r)⟩, where
nsi(→x) denotes the number of individuals of species si
at location →x

F(→r)=∑iGsi,si(→r)∑i,jGsi,sj(→r),

(1)

where the sums extend over all species in the ecosystem
Azaele et al. (2016). Eq. (1) can be used
to estimate the β-diversity as the ratio between the
couples of conspecific over the total number of couples in a sample.

Let us now study the evolution equation of F(→r,t) for the
voter model with speciation and
NN dispersal. Although we shall focus on the 2D case,
it is useful to present the general calculation
in D dimensions. Following
Chave and Leigh (2002); Zillio et al. (2005); Azaele et al. (2016) we write

F(→r,t+1)

=

(1−2N)F(→r,t)+

(2)

+

1−νDND∑k=1[F(→r+→ek,t)+F(→r−→ek,t)].

The first term in the r.h.s. of Eq. (2) represents the
fact that F does not change if two generic individuals at distance
→r are not removed in a given time step and
therefore survive. The second term represents the events in which
one of the two individuals dies (with prob. 2/N), no speciation
occurs (with prob. 1−ν) and the dead individual is replaced by a
conspecific from the 2D neighbor sites. Taking the continuous limit
N→∞ with the lattice spacing a→0, the speciation
probability ν→0, and a finite value of κ2=2Dν/a2, one obtains
at stationarity the differential equation

1rD−1ddrrD−1dFdr−κ2F(r)+cδD(r)=0

(3)

where δD is the D−dimensional Dirac delta, and because of
isotropy the β-diversity F(r) is now function of r=|→r|
only. The solution of
Eq.(3) is Azaele et al. (2016)

F(r)=cκD−2(2π)D/2(κr)(2−D)/2K(2−D)/2(κr),

(4)

where Kz is the modified Bessel function of the second kind of
order z and the constant c is fixed by the condition
∫r<adDrF(→r)=1. We recall that Eq. (4) is a
continuous expression, valid for distances much larger than the
lattice spacing Chave and Leigh (2002). Although we derived
Eq. (4) for NN dispersal, the same results
hold for a general dispersal kernel for distances larger than the
kernel range, provided that the kernel range is finite.

For D=2, Eq. (4) implies that F(r)∝K0(κr), which is characterized by a slow logarithmic decay, ∼−ln(rκ), up to distances of order 1/κ∼1/√ν,
followed by a faster, exponential falloff. Remarkably, the
β-diversity empirically measured in several tropical forests in
Central and South America is consistent with a logarithmic decay for
large distances Condit et al. (2002). We remark that this
logarithmic decay is the signature that D=2 is the critical
dimension for the voter model. In contrast, in D=1,
Eq. (4) becomes F(r)∝√rκK1/2(κr)∼exp(−rκ).
We mention for later convenience that, in D=1 with
NN dispersal, Eq. (2) can be solved
without using the continuous approximation,
giving Zillio et al. (2005)

F(r)=exp(−α(ν)r),withα(ν)=ln[(1−ν)(1−√(ν(2−ν))],

(5)

where α(ν)≈√(2ν) for ν→0.

Although the β-diversity decays exponentially on scales
1/κ∼1/√ν both in 1D and 2D, there are
important differences. Because 2D is the critical dimension, a large
biodiversity (i.e. a large average number of species) can be sustained
by very low values of the speciation rate ν. This implies that in
2D there are many species living on scales much smaller than
1/κ, where the correlations decay logarithmically. Conversely,
in 1D to maintain biodiversity one needs a large value of ν, so
that 1/κ is the only characteristic scale and there is no
additional structure on scales smaller than 1/κ. This
crucial point will be further elucidated in the rest of the section,
where we will discuss other observables in 2D (subsections
II.4 and II.5) and compare them with their
1D counterparts.

ii.4 Species-Area Relationships

We now focus on the SAR, defined as the average number of species, S
of a given taxonomic level occupying a given area of size A. SARs
are widely studied as a measure of spatial biodiversity
and quantify how larger habitats support more species than
smaller ones Rosenzweig (1995). Empirical measures of SARs at
multiple scales often reveal three different regimes
Preston (1960); Rosenzweig (1995); Hubbell (2001). At small areas, the
number of species increases rather steeply, nearly linearly, with the
sampled areas. A similar steep increase is observed at very large,
continental scales. Instead, at intermediate scales, a slower,
sublinear growth is often found. Such a growth is well
approximated by a power law S∼Az, z<1, over a wide range of
taxa Arrhenius (1921), though a logarithmic behavior
S≈ClnA has also been proposed. An extensive meta-study by
Drakare et al. Drakare et al. (2006) reconsidered a large body of
SAR studies from the literature, revealing that the power law provides
a better fit in about half of the cases. This study
also observed that the exponent z correlates positively with the
body size of the considered group of species, so that small
microorganisms typically display very shallow SAR curves as compared
with larger organisms (see also Horner-Devine et al. (2004) and
Sect. III.1).

Simulations of the (dual) voter model with speciation yields SARs
qualitatively similar to those obtained from field data, see
Fig. 3a. In the voter model, the steep initial regime
is mostly determined by the dispersal range K. For areas
significantly larger than K2, a sublinear growth is observed (see
Fig. 3b. In this regime, the growth
becomes progressively more shallow as the speciation rate ν is
decreased. For larger scales, the logarithmic slope
of the SAR curves become steeper again. The area at which this final
crossover occurs increases as ν is decreased.

An interesting question is whether the sublinear growth regime in the
voter model can be characterized by a power-law S∼Az and, in
this case, what is the value of the exponent z as a function of
ν. To address this question, we begin by
reviewing a classic estimate of z by Durrett and Levin
Durrett and Levin (1996) relying on duality (see
Sect. II.2). The speciation rate ν sets a time scale
1/ν which also corresponds to a characteristic length scale
ξ=1/√ν because of the diffusive behavior of random walkers
in the dual model. Walkers with an initial separation much
larger than ξ are likely to be annihilated before
coalescence occurs. This observation alone explains the linear scaling
of S(A) for areas A≫ξ2=ν−1. At these scales, species
are uncorrelated, as can also be inferred from the analysis of the
β-diversity in the previous section. For a system of coalescing
random walkers in 2D, the density of occupied sites ρ(t) decays
asymptotically as Bramson and Lebowitz (1991); Peliti (1986)

ρ(t)∼lntπt.

(6)

Figure 3: Species Area Relationships (SAR) and their scaling behavior
in the voter model with speciation. a) Number of species S as a
function of the sampled area A for different speciation rates as
in the caption. The triphasic shape is evident for larger
speciations rates. Simulations from Pigolotti and Cencini (2009) were
performed with a square dispersal kernel, i.e. P(→r) is a
uniform distribution on a square of side K centered on the empty
site, with K=7. b) Local slopes, dlnS/dlnA for the curves
shown in panel a. c) Dependence of the exponent z on ν as
obtained from the local slopes for both the square kernel with
K=7 and NN dispersal. The exponent is estimated at the
inflection point of the SADs, i.e. at the minimum of the local
slopes. Also shown is the prediction of Eq.(8) (black
solid line) where the black triangles correspond to the values
provided in Durrett and Levin (1996). d) Plot of 1/z vs
ln(ν) of the same data of panel c to highlight the
logarithmic behavior of Eq. (9)

The characteristic logarithmic coarsening of clusters observed in the 2D
voter model without speciation can be related to the logarithm
appearing in Eq. (6) Dornic et al. (2001). Assuming
ν≪1, the annihilation rate at time t in an area
ξ×ξ can be approximated as the annihilation rate per
walker ν times the number of walkers in the absence of
annihilations ξ2ρ(t). Integrating over time, we find that the
total number of annihilations, i.e. the total number of species, is
Bramson et al. (1996)

S(ξ2)

≈

νξ2∫1/ν=ξ2t0dtρ(t)=ln2(ξ2)−ln2(t0)2π≈

(7)

≈

2π(lnξ)2,

where t0 is the time at which the asymptotic expression
(6) starts to be valid. The upper temporal cut-off is set to
1/ν (with 1/ν=ξ2) because the number of killing events
occurring after a time ∼1/ν is bounded by the number of
walkers in the system, which is ξ2ρ(1/ν)∼lnξBramson et al. (1996). Finally, combining Eq. (7), the
fact that S(1)=1 and matching a power law behavior S=Az in the
range of scales from A=1 to A=ξ2, one finds
Durrett and Levin (1996)

z=ln[S(A)]ln(A)=2ln[ln(1/√ν)]+ln(2/π)ln(1/ν).

(8)

Also in this case, the logarithmic dependence of the exponent z
on ν derives from the fact that D=2 is the critical dimension for
the voter model.

More recent results disputed the validity of Eq. (8).
Scaling arguments hinted that z should approach a finite value
z≈0.2 in the limit of vanishing ν (see Zillio et al. (2005)
and Sec. II.5.1), while numerical simulations suggested a
power law dependence, z∼ν0.15Rosindell and Cornell (2007). Finally, further numerical simulations, based on
the dual representation of the voter model with speciation (see
Sect. II.2) and spanning a very wide range of speciation
rates from 10−3 to 10−11 confirmed the logarithmic behavior
predicted by Eq. (8) Pigolotti and Cencini (2009). The exponents
measured in such simulations, shown in Fig. 3c, are
well fitted by a phenomenological expression of the form

z=1q+mln(ν)

(9)

which is consistent with Eq. (8) up to order lnlnν,
see also Fig. 3d. However, fitted values of the prefactors
q and m are not consistent with Eq. (8). This discrepancy is
probably due to pre-asymptotic effects as well
as to the approximation of assuming a power-law range between A=1
and A=ln(1/ν).

Let us briefly discuss the role of
the dispersal kernel. As illustrated in Figs. 3c
and 3d, a comparison between NN dispersal and a
square dispersal kernel of range K=7 demonstrates that the exponent
z depends to some extent on the dispersal kernel. However,
numerical evidence Rosindell and Cornell (2007); Pigolotti and Cencini (2009) suggests that
when the dispersal kernel range is large enough (approximately
K≥5) the exponents are very weakly dependent on K. Moreover,
SARs obtained with different values of K can be rescaled onto a
universal function of A and ν via the
transformation S=f(A,ν,K)=Kχϕ(A/Kχ,ν) with a fitted
value of χ≈1.97. To the best of our knowledge, a formal
derivation of this scaling law and of the exponent χ is currently
an open problem.

The non-trivial area dependence of the SAR results is
a special feature of the critical dimension D=2. To highlight this
point, we now discuss the D=1 case as comparison. This
case is also relevant to describe quasi one-dimensional ecosystems,
such as river basins Suweis et al. (2012a). For simplicity, we limit
ourselves to the case of NN dispersal.

To the best of our knowledge, also in D=1, an exact
expression for the average number of species, S(L), in a segment
of length L is unknown. Nevertheless, it is possible to provide a
lower and upper bound for S(L). In D=1, the density of walkers
behaves as ρ(t)∼1/√t, to be contrasted with
eq. (6) valid in the 2D
case. Dimensional arguments then suggest that the
average number of species must a function of L√ν only,
i.e. S(L;ν)=Ψ(L√ν). Computational results
(Fig. 4a and inset) support well this simple argument. As
shown in the figure, the non-trivial power-law regime characteristic
of 2D SARs is absent in D=1. Indeed, the function Ψ is
linear for large arguments, with a coefficient around 1.2 and it
is nearly constant for L√ν≪1.

Figure 4: Species Area Relationshipt for the voter model in D=1. a)
Average number of species S versus the system size L for
different ν as labeled. Inset: same curves plotted vs
L√ν, notice the excellent collapse. b) SAR for
ν=10−5 compared with the theoretical upper
(10) and lower (11) bounds.

We can derive an upper bound to S(L) using that,
in D=1, individuals are organized in Ns(L;ν) segments of
conspecific individuals, so that S≤Ns, with the
equality holding if no species is present in more than one segment.
We compute Ns from the probability Pi−1,i≡F(|i−j|),
with F(r) given by Eq. (5), that two sites i
and j are occupied by conspecific individuals
Derrida and
Jung-Muller (1999)

S≤Ns

=

L−L−1∑i=1Pi−1,i=L−(L−1)F(1)=

(10)

=

L−(L−1)e−α(ν),

which for ν→0 can be approximated as Ns≈1+√2ν(L−1).

The lower bound follows from Jensen’s
inequality (see also Derrida and
Jung-Muller (1999))
applied to the frequency of species represented by the individual in
site i∈[0,L−1], here denoted ni(L), which yields

S=∑i⟨1ni(L)⟩≥∑i1⟨ni(L)⟩,

(11)

where ⟨ni⟩=∑jPi,j and Pi,j=F(|i−j|)
is again given by Eq. (5) and can be easily summed
numerically.

In Fig. 4b we compare the numerically
obtained SAR with the upper (10) and lower
(11) bounds. Notice that the upper bound is very close to the
actual SAR, implying that most species are organized in single
segments.

ii.5 Species-Abundance Distributions

We now discuss Species-Abundance Distributions (SADs), P(n;A), that
measure the relative abundance of species in a given area A. More
precisely, denoting S(A) the total number of species sampled in an
area A, each composed by ni (i=1,…,S(A)) individuals,
P(n;A)dn is the probability that a randomly picked species has an
abundance between n and n+dn. While the expression of P(n;A) for
well-mixed neutral models is known Volkov et al. (2003), computing it for
spatially explicit models, such as the voter model with speciation,
has proven to be a rather hard problem. We first discuss in section
II.5.1 an approach based on standard
finite-size scaling, and
underline its limitations. In Sec. II.5.2, we discuss how
this approach can be generalized at the critical dimension, present
numerical results, and discuss a recent attempt to compute P(n;A)
exploiting duality. Although we focus ond
comparing the scaling theory with results from the voter model with
speciation, we remark that the theoretical approach presented
in this section is more general and can be applied to a vast class
of models at the critical dimension.

Power-law scaling relation

In the voter model with speciation, the SAD is not only a function of the system size
A, but also of the speciation rate ν. Although we are mainly
interested in 2D, it is instructive to consider the general case in
which A=LD, where L is the linear size of the sample. Following
Zillio et al. (2005); Azaele et al. (2016), we assume a standard scaling
form for the SAD

P(n;A,ν)=n−βΨ(nνα,AνD/2)

(12)

where the exponents α and β remain unspecified for the
time being, whereas the exponent D/2 stems from the diffusive nature
of neutral models ν∼t−1∼L−2∼A−2/D. Note
that in models with long-range, non-diffusive dispersal
Rosindell and Cornell (2009) the scaling form might
differ. Equation (12) describes a power-law dependence of
P on n, holding up to a scale determined by the scaling function
Ψ, that depends on dimensionless combinations of the population
size n, the speciation rate ν, and the system size A. To the
best of our knowledge, there is no available analytical prediction for
the exponent β. The exponent α can be estimated in the
dual formulation of the voter model with speciation, where the
population size n is the number of coalescences that occur before an
annihilation (see Sec. II.2). This implies that α
is the same exponent characterizing the temporal decay of the
density of coalescing random walkers, ρ(t)∼t−α.
However, ρ(t) decays as ρ(t)∼t−min(1,D/2) for
D≠2 and ρ(t)∼log(t)/t in D=2, see
eq. (6) and Bramson and Lebowitz (1988); Peliti (1986). Consequently, one should expect
the power-law scaling of Eq. (12) to hold in D=1 and
D≥3, but not at the critical dimension D=2, where logarithmic
corrections should appear.

Generalized scaling relation

The dependence on ν was omitted as the above scaling law was
applied to observational data for which the speciation rate is unknown
and assumed to be fixed. The key aspect of Eq. (13) is that
f and g, are general functions and not necessarily power-laws as
in conventional scaling, allowing for the possibility to include
logarithms or other functional dependencies. The scaling function
Ψ(x) is still assumed to be a power law

Ψ(x)∼x−Δ

(14)

for small values of x, where Δ is an exponent to be
determined. Thus, also Eq. (13) postulates a power-law
dependence on n, but with a more general cut-off for large
areas. After specifying the functions f and g, Eq. 13 can
be tested by plotting P(n;A)/g(A) versus x=n/f(A) for a set of
different areas and assessing the quality of the data collapse onto a
single curve, Ψ(x).

To determine the functions f and g, we impose that P(n;A)
has to be normalized,
∫∞n0dng(A)Ψ(n/f(A))=1, and that its average
value has to be
⟨n⟩=∫∞1dnng(A)Ψ(n/f(A)).
Substituting the scaling form (14) into these two equations,
it is possible to derive conditions that the functions f and g
must obey, depending on the value of Δ. In particular, the case
Δ=1 is marginal and needs to be treated with care (other values
Δ≠1 are analyzed in the Appendix). Approaching such a
limit as Δ=1−ϵ with ϵ≪1,
Eq.(14) becomes

Ψ(x)=x−1+ϵ∼1x[exp(ϵ)ln(x)]∼1x[1+ϵln(x)]

(15)

up to first order in ϵ. At the same order in ϵ, the
two conditions for P(n;A) become
1∼g(A)f(A)ln(f(A))[1+ϵ2ln(f(A))]
and ⟨n⟩∼g(A)f(A)2, respectively, from which we
finally obtain

f(A)

=

⟨n⟩ln⟨n⟩[1+ϵ2ln⟨n⟩]

g(A)

=

1⟨n⟩ln2⟨n⟩[1+ϵ2ln⟨n⟩]2

(16)

up to first order in ϵ. Notice that both functions f and g include logarithmic
corrections. By means of a similar calculation, one can estimate the
k-th moment ⟨nk⟩, and verify that all the moment
ratios ⟨nk⟩/⟨nk−1⟩ scale in the
same way, up to a multiplicative constant

⟨nk⟩⟨nk−1⟩=∫dnnkP(n;A)∫dnnk−1P(n;A)∝f(A)k≥1.

(17)

revealing a highly anomalous scaling.

Zillio et al.Zillio et al. (2008), showed that this scaling form
provides a much better collapse of empirical data from the Barro
Colorado tropical forest than a power-law scaling relation such as
Eq. (12). This supports the idea that Δ is close
to its marginal value 1 in tropical forests.

Figure 5: SAD and
data collapse. Results are presented for different linear
system sizes and different speciation rates ν, keeping the
product Aν=200 constant. a) SADs for different linear sizes
from L=400 to L=2500. b) Collapse of SADs by means of
Eqs.(13) and (16). The fitted parameter in the
functions f and g is ϵ=0.08. c) Naive collapse
without logarithmic corrections, where deviation for perfect
collapse are evident. d) Collapse with the scaling form of
Eqs.(13) and (16), but setting
ϵ=0. Also in this case the discrepancy is evident.

We tested computationally whether Eqs. (13) and (16)
provide a good collapse of SADs obtained from the voter model with speciation
and whether the relationship between the moments,
Eq.(17), holds. In simulations, an additional
parameter is the speciation rate ν. As discussed above, ν
appears in scaling relationships via the dimensionless combination
AνD/2, that in 2D equals Aν. Thus, although
Eqs. (13) and (16) do not include speciation
explicitly, we expect these relationships to hold if Aν is kept
constant. We therefore performed computational analyses fixing
Aν=200, although the conclusions are robust against this choice.
Results are summarized in Figure 5 which shows plots of
the SAD, for systems with different linear size, L and different
speciation rates ν (with L2ν=Aν=200). Observe in
Fig. 5a that the smaller the size (or the larger the
speciation rate) the smaller the maximal
abundance. Figure 5b show the data collapse as given by
Eqs. (13) and (16), where ⟨n⟩ is the
average number of individuals measured in each area A and ϵ
is a free parameter that we fitted obtaining ϵ=0.08 and
a remarkable collapse of the different
curves. The small value of ϵ, is consistent with the assumed
small deviation from Δ=1. A similar collapse for Aν=20
leads to an even smaller value ϵ≈0.069 (not shown).
We verified that either removing all logarithmic corrections (thus
plotting results as a function of ⟨n⟩) or simply
fixing ϵ=0 in Eq. (13) and (16) leads to less
convincing collapses, as shown in Fig. 5c and
5d, respectively. Clearly, these deviations can pass
unnoticed in the presence of statistical fluctuations. Probably, this
is the reason why in Rosindell and Cornell (2013) a simple scaling law was claimed to
hold for the 2D voter model with
speciation. Finally, we also verified that moment ratios scale as
f(A), as predicted by Eq.(17) and illustrated in
Fig.6.

Figure 6: Moment ratios for different values of
k. As predicted by Eq.(17), in the case Δ≈1 all moment
ratios ⟨nk⟩/⟨nk−1⟩ scale in the
same way with f(A) up to a multiplicative constant. As in
Fig. 5, the fitted value is ϵ=0.08.

In summary, a non-standard scaling form, including logarithmic
corrections, provides an excellent collapse both for empirical data
and for numerical simulations of the 2D voter
model. We remark that the scaling theory is
phenomenological, and the small parameter ϵ controlling the
importance of logarithmic corrections is, at this level, a
non-universal free parameter. These results are in sharp contrast
with the one-dimensional case, where logarithmic corrections are not
expected. Indeed, Fig. (7) shows that the naive
scaling form P(n;A)⟨n⟩ vs. n/⟨n⟩
(derived in Appendix A for the case Δ≠1) yields a perfect
collapse for SADs in one-dimensional systems.

It is interesting to remark that the data collapsed in
Zillio et al. (2008) were obtained from tropical forests of different
areas A. It is reasonable to assume that the speciation rate ν
do not vary much among these forests. Therefore, the product Aν
is not fixed, as in our computational analyses. A possible
explanation is that, although the collapse achieved in this way is
not perfect, the deviations from perfect scaling are too small to be
appreciated in observational data due to the limited sample size. We
have verified in simulations (not shown) that
keeping ν constant (rather than Aν constant) small
deviations from perfect collapse are observed.

We conclude this section mentioning that a heuristic expression for
the SAD has been recently derived for the voter model with speciation
following a completely different approach
Danino et al. (2016a); Shem-Tov et al. (2017). Let us define P(x,t) as the
distribution of the number of individual of a given species at time
t. If we approximate x as a continuous quantity, we can
heuristically write a Fokker-Planck equation for the evolution of
P(x,t)

∂tP(x,t)=ν∂x[xP(x,t)]+∂2x[I(x)P(x,t)]

(18)

where the first term in the right hand side is the negative drift due
to speciation, and the second is the fluctuation in population size,
where I(x) is the average number of interfaces of a species of size
x. The crucial underlying approximation is to neglect fluctuations
of I(x), which is appropriate if the distribution
of the number of interfaces at fixed value of x is a very peaked
function. In this simple framework, all the dependence on the spatial
dimension of the voter model is recap into the function I(x). The
steady-state solution of Eq. (18) is

Pst(x)=e−ν∫dxxI(x)I(x).

(19)

From duality considerations Danino et al. (2016a); Shem-Tov et al. (2017), the average
number of interfaces must scale in 2D as I(x)=x/(1+clnx) where
c is a non-universal constant. Notice how the expression of I(x)
includes familiar logarithmic terms and that the
constant c plays the role of the exponent ϵ in the
scaling theory. Substituting this expression into
Eq. (19) leads to an explicit expression for the SAD,
which obeys a scaling law with logarithmic corrections similar to
Eq. (16), though not identical. A more detailed comparison between
this result and the previous scaling form is an interesting issue, but
beyond the scope of this review.

ii.6 Species persistence-times

So far, we have considered neutral predictions of
static ecological observables. However, neutral theory can also be
used to predict time-dependent properties. A chief example is the
distribution of survival times. The survival time τ (also
called ”persistence time”) within a geographic region is defined as
the time occurring between the speciation event originating a given
species and its local extinction
Pigolotti et al. (2005). Recent empirical work on north-american
birds and herbaceous plants revealed that the probability of
observing a persistence time τ decays as as power laws
P(τ)∼τ−1.83 and P(τ)∼τ−1.78
respectively, with area-dependent exponential cut-offs
Bertuzzo et al. (2011); Suweis et al. (2012b).

In the voter model with speciation, the survival
probability as a function of time can be computed analytically.
Also in this case, the calculation relies on duality
Bramson and Lebowitz (1991); Peliti (1986); Lee (1994). In 2D and in the limit of vanishing
ν one obtains

P(τ)∼lnττ2

(20)

while standard power-law scaling P(τ)∼τ−1/2 is
expected in 1D. For non-negligible values of ν, these scaling
forms are cut-off by a ν-dependent exponential factor
exp(−ντ) in either dimension. Also in this case, diffusive
scaling relates the characteristic time scale 1/ν with a length
scale ξ via ξ∼√ν. This explains the aforementioned
area-dependent cut-offs observed in empirical data Bertuzzo et al. (2011).

Figure 8: Species persistence times. (a) Probability distribution
function of species persistence times for different values of the
speciation rate ν as in label. (b) and (c) show the pdf
rescaled with the logarithmic correction, P(τ)τ2/lnτ,
and with a power law, P(τ)τ1.9,
respectively.

Species persistence times in simulations of the 2D voter model with
speciation are shown in Fig. 8a. Panels (b) and
(c) show compensated plots of the simulation results. The simulations
support the prediction of eq. (20) (panel b), and also
illustrate that a power law with an exponent close to 2 (1.9 in
this case) provides a good approximation of the scaling predicted by
Eq. (20) in a broad range of scales (panel c), consistently
with the empirical findings in Bertuzzo et al. (2011); Suweis et al. (2012b).

In the voter model with speciation, the habitat is saturated and each
site is always occupied by an individual. In this section, we study
neutral spatial models where the number of individuals that can
inhabit a site is varied. We consider three
variants: the stepping-stone model with speciation,
where each site can host many individuals but the landscape
remains saturated; the contact process with speciation, where
occupancy is limited to a maximum of one individual per site, but
sites can also be empty; and the O’Dwyer-Green model, where occupancy
is unbounded.

iii.1 Stepping-Stone Model with speciation

In the voter model, each lattice site hosts a single
individual. This assumption is appropriate for big sessile
species, such as trees, where each individual occupies a
well-defined area and exploits its local resources. On the other
side of the spectrum, microorganisms, such as small eukaryotes or
bacteria, are often present in very large numbers on tiny spatial
scales, where all individuals share the same resources. For these
species, it is more appropriate to think of the habitat as subdivided
into small patches, connected by migration and each hosting a large
number of individuals directly competing with each other
Fenchel and Finlay (2004). To model such ecological cases, in
this section we consider the stepping-stone model
Kimura and Weiss (1964); Korolev et al. (2010) with speciation, which
generalizes the voter model with speciation to the case in which each
site hosts a fixed number M of individuals.

Similar to the voter model with speciation, at each time step an
individual is randomly chosen and killed. With probability ν, it
is replaced by an individual of a novel species. With complementary
probability (1−ν), a reproduction event occurs. The parent of the
new individual is selected with probability (1−μ) among the
surviving M−1 individuals present at the same site, and with
probability μ among the M individuals in a randomly chosen
neighboring patch (according to a probability distribution on the
neighbors P(→r), similar to the case of the voter model). The
particular case of M=1 reduces to the voter model with speciation up
to a time rescaling t→μt. Like the voter model, the stepping-stone
model admits a dual representation in terms of coalescing random
walkers with annihilation, which can be exploited for efficient
numerical simulations. The main difference with respect to the
dual of the voter model is that, in the dual stepping-stone model,
at each step a random walker can either move to another site or stay
in the site of origin. Coalescence can happen in both circumstances,
corresponding to reproduction of an individual from neighboring sites
or from the same site. For full details on the implementation we refer
to Cencini et al. (2012).

As revealed by numerical simulations of the stepping-stone model based
on the dual representation, SARs are qualitatively similar to those of
the voter model, although the exponents z are, in general, smaller
than in the voter model Cencini et al. (2012). In particular,
the exponent depends not only on ν, but also on the combination of
parameters Mμ, which determines the regimes of the model. For
Mμ≪1, each local site is likely to contain only one species. In
this limit, each site behaves as one individual up to a time
rescaling, so that one should expect the same exponents as in
the voter model with speciation. In the opposite limit Mμ≫1,
there is a large diversity of species at each site. An analytical
argument suggests that, in this latter limit, the exponent should be a
factor two smaller than in the former limit
Cencini et al. (2012). Let us study the limit
Mμ≫1 in the dual representation. Since random walkers in the
same site have a low probability of coalescence, they will wander for
a long time before coalescing. Therefore, we can assume that,
asymptotically, they will behave as in the well-mixed case. This
implies that their density in an area smaller or equal than ξ2
approximately decays according to the mean-field formula

ρ(t)∼1t.

(21)

Observe that in this case the characteristic length is
ξ=√μ/ν, as random walks diffuse with
probability μ at each time step. Proceeding as in Eq. (7), the
average number of species in an area ξ2 can be estimated as

S(ξ2)∼νMξ2∫μ/ν=ξ2t0dtt=Mμln[ξ2t0].

(22)

To compute z, we also need an estimate for S(1), that in this case
is not trivially equal to one. As the population is assumed to be
well-mixed in an area equal to ξ2 or smaller, the composition of
a single site can be thought as a sample of M individuals from this
well-mixed population. The probability distribution of the abundance
in such a sample is given by Ewens’ sampling formula
Ewens (1972). Substituting its expression yields

S(1)=M−1∑j=0MμMμ+j≈Mμln(1+μ−1).

(23)

Combining Eqs. (22) and (23) and assuming a
power law in the range from A=1 to A=ξ2, we find an exponent

z∼ln(ξ2)lnln(ξ2)=lnln(ν/μ)ln(ν/μ)

(24)

which, to the leading order, is a factor 2 smaller than the
corresponding estimate for the voter model (8). The decrease
of the exponent z with the combination of parameters Mμ is
confirmed in numerical simulation, see Fig. 9, although
the asymptotic reduction is less than the factor two predicted by the
approximate estimate of eq. (24).

Figure 9: Species-area exponents for the Stepping Stone Model at
fixed ν=10−6, different local population size M and
dispersal rate μ, with NN dispersal. The
numerical estimate of the exponent z in the voter model for
NN-dispersal and the same value of the speciation rate is shown
for comparison.

Summarizing, the stepping-stone model at large local community size
M yields smaller values of the species-area exponent z than the
voter model Cencini et al. (2012). This fact is consistent with
the ecological observation that microbial communities, characterized
by very large local community sizes, typically display very shallow
species-area relations, and that in general there seems to be a
positive correlation between the exponent z and the body size of a
taxonomic group Horner-Devine et al. (2004). In the stepping-stone model, a
decrease in the SAR exponent is observed in the regime
Mμ≫1 where each site hosts a large number of species and
therefore provides a buffer for biodiversity
Cencini et al. (2012). This interpretation is also consistent
with the “cosmopolitan” nature of many microbial species, i.e. the
fact that relatively small communities of microbes host a biological
diversity comparable with that observed in the whole planet
Finlay and Fenchel (2004); Fenchel and Finlay (2004). This feature has sometimes been
explained invoking the fact that microbes have the possibility of
long-range dispersal Finlay and Fenchel (2004). However, numerical
simulations show that, in the voter-model with speciation, long-range
dispersal leads to steeper, rather than shallower SARs Rosindell and Cornell (2009).

iii.2 Contact Process with speciation

In the voter model, every dead individual is instantly replaced by a
newborn, leading to a constantly saturated environment. The implicit
underlying assumption is that the birth rate is infinite, so that
death events are the rate-limiting steps. Such assumption
constitutes a good approximation in resource-rich ecosystems. In
less rich ecosystems, where the birth rate is finite, the environment
is not always saturated and empty gaps can exist
Loreau (2000).

To explore this latter case, we study here the contact process with
speciation, which is the multi-species
variant of the well-known contact process
Griffeath (1981); Holley and Liggett (1975); Durrett and Levin (1994); Marro and Dickman (1999); Hinrichsen (2000). As
usual, we consider the model on a 2D square lattice. Sites of the
lattice can be occupied by individuals belonging to different species
or empty. The model is defined in continuous time; each individual
dies at a rate d and reproduces at a rate b. In case of a death,
the site is simply left vacant. A reproduction event is considered
successful only if the individual has at least one vacant neighboring
site. In such a case, one of the vacant neighboring sites is
chosen at random. With probability ν, the site is occupied by an
individual of a new species (speciation event); with complementary
probability, (1−ν), the newborn is of the same species as the
parent.

As in the standard contact process
Griffeath (1981); Holley and Liggett (1975), the parameter determining
the asymptotic density of occupied sites ρ is the dimensionless
birth-to-death ratio η=b/d. For
η<ηc≈1.649 the absorbing state in which all sites are
empty is stable. A non-equilibrium phase transition at η=ηc
separates this region from a stable active phase (η>ηc)
characterized by a non-vanishing value of ρ that depends on
ηMarro and Dickman (1999); Hinrichsen (2000). For η→∞
one has ρ→1 and the model is equivalent to the voter
model with speciation Durrett and Levin (1994).

Figure 10: (top) Snapshots of configurations of the contact process with
speciation at different values of the birth-to-death rate ratio
η and ν=10−4. In each panel, each color represents a
different species. (bottom) SARs at different
values of the birth-to-death rate η (shown in the figure
legend) and ν=10−5. (inset) Red dots: estimated exponent z
as a function of ν for the contact process with speciation at
η=1.68. Red dashed line is a linear fit; black dashed line is
the corresponding fit for the voter model with speciation for
comparison. We have chosen a NN dispersal kernel in
all panels.

The CP is a self-dual model. Therefore, duality cannot be
exploited in numerical simulations as in the case of the voter model. Forward
simulations show that the SAR and the corresponding exponents are
remarkably similar to the voter model with speciation even at small
values of η, corresponding to very fragmented ecosystems as shown
in Fig. 10. For values of η very close to ηc (but
within the active phase) and small values of ν, SAR exponents tend
to be smaller than in the voter model, see inset of Fig. 10.

In principle, in a very fragmented ecosystem it would not make sense
to sample empty areas, or areas with very few individuals. With this
idea in mind, an alternative to the standard definition of SAR used so
far is to weigh the sample of a given area with its
number of individuals, i.e. of occupied sites. Adopting this
definition one finds qualitatively different SARs for small values of
ηCencini et al. (2012). In particular, these SARs do not
seem to be characterized by a clear power-law range. We refer the
reader to Ref. Cencini et al. (2012) for a broader discussion
of this issue.

iii.3 O’Dwyer-Green model

We have seen that finding exact results for neutral spatial models
constitutes a formidable problem, and even in the simple case of the
voter model only asymptotic results are known.

To make progress in this direction, O’Dwyer and Green proposed a
spatial neutral model in which individuals do not compete, i.e. the
site occupancy is not bounded O’Dwyer and Green (2010). In their model, each
individual can reproduce at a rate b, giving rise to a newborn
located according to a dispersal distribution, die at a rate d, or
speciate at a rate ν, giving rise to a newborn of a new
species. The model was studied at the critical point b+ν=d. The
lack of interaction considerably simplifies the mathematical
treatment: the model can be mapped into a field theory from which the
authors of O’Dwyer and Green (2010) obtained an analytical expression for the
species-area law and the dependence of z on ν. In particular,
the solution was derived by writing an equation for the distribution
of a generic species, which was solved by imposing detailed balance.
However, Grilli and coworkers Grilli et al. (2012) pointed out a
flaw in this procedure. In this model all species are transient, as
the birth rate of each species is always smaller
than the death rate because of speciation. This implies
that all species eventually go extinct, so that the detailed balance
(i.e. equilibrium) assumption is not valid.

An often overlooked aspect of the O’Dwyer and Green model is the lack
of a carrying capacity. Although well-mixed neutral
models commonly do not have a carrying capacity (beside that of the
entire ecosystem), a local carrying
capacity, i.e. a maximum occupancy of each lattice site, is
a standard ingredient in spatial neutral theory, shared by all models we discussed so far.
In the O’Dwyer and Green model, since the dynamics of the entire ecosystem is
a critical branching process, the population at each site undergoes
huge fluctuations. This fact implies as a drawback that numerically
simulating the steady-state of the model and sampling its
configurations is extremely difficult. While the authors of
Grilli et al. (2012) clearly pointed out that the detailed balance
solution leads to several inconsistencies and is therefore not valid,
to the best of our knowledge there have been no attempt of comparing
this solution with numerical simulations to see if detailed balance
can provide a reasonable approximation of the dynamics in some
particular regimes or limits.

Currently, the research of spatial neutral models that can be solved
analytically is still open O’Dwyer and Cornell (2017). In this
direction, although this review focuses on lattice models, we mention
a recent phenomenological attempt based on a spatial Fokker-Planck
equation where both space and population sizes are continuous
variables Azaele and Peruzzo (2016).

In the previous sections, we focused on neutral
ecological models. However, in real ecosystems the neutral
assumption is (at best) a crude approximation. It is thus
interesting to examine some of the main effects of non-neutral
forces, also because many biodiversity patterns
that are well predicted by neutral models are also found in richer,
non-neutral models McGill (2003); Tilman (2004); Gilbert and Lechowicz (2004).
A main difficulty in comparing neutral and non-neutral models is the
large number of possible ecological effects (and corresponding
parameters) that typically enter the latter. In
this section, with the aim of understanding basic non-neutral
effects in a simple setting, we present a minimal model introduced
in Pigolotti and Cencini (2010), where one can continuously move from a neutral
to a non-neutral scenario by varying a single parameter, tuning
the amount of spatial disorder. We then discuss generalizations to
other types of spatio-temporal disorder.

iv.1 Habitat-preference model

We consider a variant of the voter model where
different sites are preferred habitats for each one of the competing
species. For the sake of simplicity, we limit ourselves to the case
of two species A and B with NA and NB individuals,
respectively. We assume habitat saturation, so that the total
population is N=NA+NB=L2 where the system is a square lattice of
size L with periodic boundary conditions. Individuals of type A
and B can also migrate to the system from an infinite reservoir
where they are equally
represented. Each lattice site can be of type a or b,
i.e. being a preferred habitat for colonization by species A or B,
respectively. After colonization, mortality and dispersal do not
depend on being on a preference site. Ecologically, this means that
the fitness advantage belongs to the seeds and not to the individuals
themselves (see Chesson and Warner (1981) for a different choice). The a
vs b character of each site is chosen randomly at the beginning and
it remains fixed over time – quenched disorder. To maintain the model
globally symmetric, we assume equal proportions of a and b sites
and that intensity of the two biases (a favoring A and b
favoring B) are identical. The dynamics proceeds as follows. At
each discrete time step, a lattice site is randomly chosen with
uniform probability and the residing individual is killed. The
individual is replaced either by an immigrant from the reservoir (with
probability μ) or by an offspring of an individual residing in one
of the four neighboring sites (with probability (1−μ)). In both
cases, the colonization probability is biased by an additional factor
γ for the individuals that have preference for the empty site.
In formulas, the probability of colonization of a site x={a,b} by
an individual X={A,B} (Y={B,A}) having (not having) preference
for that site is

respectively, where nX (nY) denotes the number of individuals of
species X={A,B} (Y={B,A}) in the neighborhood of the
considered site. Similar models have been proposed also in the context
of heterogeneous catalysis Frachebourg et al. (1995) and social dynamics
Masuda et al. (2010). For γ=0 and μ=0, the standard (neutral)
voter model with two species is recovered. For γ=0 but
μ≠0, it corresponds to the noisy voter model
Kirman (1993); Granovsky and Madras (1995).

Also in this model, the results can depend on the choice
of the dispersal kernel P(r). Here we focus on the NN
dispersal and global dispersal (GD), i.e. a mean-field version of
(25). The GD case can be thought as a variant of the
two islands model Moran (1962) of population genetics, where each
island host N/2 individuals and is favorable to one of the two
species. In the mean-field version, the state of the system is
univocally determined by the numbers of individuals NAa and
NBb residing on their island of preference. The numbers of
individuals outside their island of preference are
NBa=N/2−NAa and NAb=N/2−NBb. The dynamics is then
fully specified by the probabilities per elementary steps that
NXx (with X={A,B} and x={a,b}) increases or decreases by
a unit:

WNXx→NXx+1

=

(12−NXxN)WxX(NA,NB)

WNXx→NXx−1

=

NXxNWxY(NA,NB)

(26)

where WxY and WxX are given by eqs. (25) with
nX and nY replaced by NX=NXx+NXy and
NY=NYy+NYx, respectively.

iv.2 Extinction times

In the absence of immigration (μ=0) and for finite populations
N<∞, persistent coexistence of the two species is not possible:
demographic stochasticity eventually drives one of the species to
extinction (the absorbing state) with the fixation (in the
jargon of population genetics) of the other species. In this
case, information on the system can be obtained by studying the
dynamics toward extinction Chesson and Warner (1981). Of particular interest
is the average extinction time, ⟨Text⟩, and its
dependence on system properties, such as the deviation from neutrality
and the population size.

In the neutral case (γ=0), as discussed, the system recovers
the voter model with NN dispersal and the Moran model Moran (1958)
in the version with global dispersal. In this limit, the extinction
time is set by the population size. In particular, for large N we
have ⟨Text⟩∼NlnN for NN-dispersal
Krapivsky (1992) and ⟨Text⟩∼N for global
dispersal Moran (1958); Gillespie (2004). To inquire the effect of habitat
preferences we performed simulations of the model (25)
with an initial condition NA=NB=N/2 until the
extinction of one of the two species.

Figure 11: Extinction times for the model with NN
dispersal without immigration (ν=0). Mean extinction time
⟨Text⟩ as a function of N for different values
of γ as in label. The blue curve approximating the neutral
γ=0 data points corresponds to the neutral expectation
⟨Text⟩∝NlnN, the black curves over the
symbols for γ≠0 correspond to exponential fits of the
form ⟨Text⟩∝exp(C(γ)N). The inset
shows (symbols) C(γ) vs γ, while the black solid line
display the best fit C(γ)=Aγβ with
β≈1.63. The average extinction time is
obtained by an annealed average, i.e. by randomizing the preference
sites at each realization. Each point represents an average over 103
realizations.

Figure 11 shows the average extinction time,
measured in generations, i.e. N elementary steps of
eqs. (25), as a function of the population size N for
different values of γ. For
γ=0 we observe the NlnN behavior expected in the neutral
case. Habitat preference (γ>0) leads to a dramatic increase of
the average extinction time, which becomes exponential in N

⟨Text⟩∝exp(C(γ)N),

(27)

for large enough N. The dependence of the constant C(γ)
on γ, shown in the inset, is well-fitted by a power-law with
exponent ≈1.63. The mean-field version of the model presents
similar qualitative features with the only difference that
⟨Text⟩∝N for γ=0 and with some
differences in the γ dependence of C(γ), as shown in
Pigolotti and Cencini (2010).

The exponential dependence of the average extinction times on N
indicates that habitat preference has a stabilizing impact on the population
dynamics. Indeed, when N is large enough, the two species coexist on
any realistic time scale. The stabilizing effect of habitat
preference reflects also in the probability of fixation Pfix,
i.e. the probability that a species, say A, gets fixated when
initially present as a fraction x=NA/N of the population. In the
neutral case, standard computation Gillespie (2004) shows that
Pfix(x)=x. As shown in Pigolotti and Cencini (2010), when γ is
increased, Pfix(x) develops a much steeper dependence on x and
quickly reaches values ≈1/2 even for small x, provided
that γ is large enough. In other words, the stabilization
due to habitat preference tends to compensate any initial
disproportion between the population of the two species.

iv.3 Coexistence

In the presence of immigration (μ>0), a locally extinct
species can recolonize, leading to a dynamical coexistence between the
two species. However, if the typical recolonization time 1/μ is
large compared to the average extinction time ⟨Text⟩, such recovery from extinction is slow and
unlikely. Therefore, most of the time the ecosystem is
dominated by one of the two species.
Therefore, the
distribution of the population size of any of the two species, P(X)
(X=A,B) is peaked at 0 and at the population size N,
corresponding to dominance of either of the two species. We denote
this regime as monodominance, see
Fig. 12a. In the opposite limit
⟨Text⟩≫1/μ, temporary extinctions are very
unlikely and the distribution is peaked at NA=NB=N/2 leading to
pure coexistence of the two species
(Fig. 12c). For intermediate values of
μ, temporary extinctions are still possible though the
replenishment due to immigration will tend to equilibrate the two
populations. In this case of mixed coexistence, the
distribution is characterized by three local maxima at NX=0,N/2,N
(Fig. 12b).

Figure 12: Different regimes of coexistence for the case with
NN dispersal and immigration for the model with
habitat preference. Top panels show the stationary distribution
P(NA) for γ=0.3 and (a) N=50 with μ=10−3, (b)
N=300 with μ=2×10−3, and (c) N=100 with
μ=10−3, corresponding to a typical distribution in the cases
of monodominance, mixed regime and pure coexistence, see
text. Bottom panels show how the three regimes partition the
N,μ-parameter space for different values of γ: (d)
γ=0 corresponding to the neutral case, (e) γ=0.3 and
(f) γ=1. The three points in (e) correspond to the
distributions displayed in the top panels, as labelled by the color
coding.

Figs. 12d,e,f show the three
regimes of coexistence in the N−μ parameter space for
the model with NN-dispersal for different habitat preference
strength γ (increasing from left to right). In the mean-field
model, we find the same qualitative features, except that for
γ=0 the mixed regime is absent, so that one has a direct
transition from monodominance to pure coexistence Pigolotti and Cencini (2010).

The main emerging feature is that increasing habitat preference
expands the region of parameter space corresponding to mixed
coexistence at the expenses of monodominance. Surprisingly, the pure
coexistence regime seems to be insensitive to the degree of habitat
preference. In particular, the critical line μc(N) separating it
from the mixed regime seems to be the same that separates
coexistence from monodominance in the neutral model (γ=0) with
global dispersal, which is given by the expression
μc(N)=2/(2+N). This result can be obtained in the
following way. For γ=0, the transition rates
(IV.1) can be expressed in terms of the rates for NA to
increase/decrease by one

WNA→NA±1

=

N2±(N2−NA)N×

(28)

×

⎡⎢
⎢⎣(1−μ)N2∓(N2−NA)N+μ2⎤⎥
⎥⎦.

Then, the equilibrium distribution P(NA) can be computed imposing
the detailed-balance condition

P(NA+1)P(NA)=WNA→NA+1WNA→NA−1,

(29)

which must hold at stationarity since the process is one dimensional
Gardiner (2009). To determine μc(N) for the transition from
monodominance to coexistence, it is sufficient to determine
whether, for small NA, P(NA) is an increasing or a
decreasing function. Using (29) with (28)
and imposing P(NA+1)>P(NA) one obtains, after some algebra, the
inequality [(2+N)μ−2](N−2NA−1)>0, which is verified whenever
μ>2/(2+N). Notice that, in the case of global dispersal, the
distribution is uniform along this line, i.e. for μ=μc one
finds P(NA)=1/N.

iv.4 Generalizations of the habitat-preference model

To gain physical insight into the different regimes shown in
Fig. 12, a variant of the habitat
preference model was introduced and analyzed for the global dispersal
case in Borile et al. (2013). By considering the first two terms of a
system-size expansion of the master equation, results in the
infinite-size limit and finite-size corrections were derived. In the
infinite-size limit, i.e. neglecting the effect of fluctuations, the
introduction of a non-vanishing local preference generates a
deterministic force, which can be described as an effective potential
V(δ) for the relative difference of densities
δ=(NA−NB)/N. This potential has a minimum at the coexistence
state, δ=0, corresponding to a maximum in the probability
distribution at NA=NB=N/2. In other words, species coexistence
emerges for infinitely large sizes. On the other hand, for finite
systems, when fluctuations are considered, the only possible true
steady states are the absorbing states δ=±1, where the
effective potential V(δ) is singular. The minimum at δ
is separated from the negative singularities by two potential
barriers. As strength of the local preference and/or N increase, the
basin of attraction of the coexistence state becomes larger and deeper
and the two symmetric barriers become closer to the absorbing states
and higher. Consequently the time needed to escape the coexistence
state becomes much longer, therefore unaccessible in computer
simulations. Thus, three different regimes can thus be identified:
the absorbing, intermediate (quasi-active) and active phase (much as
in Fig. 12). In the absorbing phase,
symmetry is broken and one of the two species reaches extinction with
certainty. This regime is equivalent to the monodominance regime in
Fig. 12. The active phase is
characterized by a coexistence of both species, and survives
fluctuations only in the infinite-size limit. This corresponds to the
coexistence phase of Fig. 12. Finally,
the intermediate state is a mixture of the two previous ones: the
absorbing states and the coexistence state are locally stable, thus,
the system is tri-stable, and the steady state depends on the initial
conditions. This is the mixed state of
Fig. 12. These results provide a nice
analytical example of how noise can effectively change the shape of a
deterministic potential. Still, the presence of absorbing states –
with the associated singularities in the steady state distribution –
prevent true phase transitions from occurring: the only possible
steady state for any finite system is an absorbing one. Only in the
infinite-size limit, noise vanishes and the coexistence state becomes
truly stable Borile et al. (2013).

Another study scrutinized the case in which there are local
habitat preference only at some specific locations in space, while
all other sites are neutral Borile et al. (2015). An interesting
example which has been analyzed in details is that of a square lattice
where only the left (resp. right) boundary has a preference for
species A (resp. B), (Borile et al. (2015), see also
Mobilia (2003); Mobilia et al. (2007)). The conclusion is that even mild biases
at a small fraction of locations induce robust and durable species
coexistence, also in regions arbitrarily far apart from the biased
locations. As carefully discussed in Borile et al. (2015) this result
stems from the long-range nature of the underlying critical bulk
dynamics of the neutral voter model, and is robust to the introduction
of non-symmetrical biases –i.e. stronger for one of the species–
except for the fact that the state of coexistence is no longer
symmetric. These conclusions have a number of potentially important
consequences, for example, in conservation ecology as it suggests that
constructing local “sanctuaries” for different competing species can
result in global increase of stability of their populations, and
thus an enhancement of biodiversity, even in regions
arbitrarily distant from the protected zones Borile et al. (2015).

iv.5 Temporally-dependent habitat preferences

We have seen that spatial quenched disorder generically fosters
species coexistence. Another important question is what happens when
the preference for a species are time-dependent, i.e. if neutrality is
temporarily broken in favor of one of the coexisting species, while
the ecosystem remains neutral on average. This question has a long
tradition in ecology. Several theoretical studies have looked at the
impact of environmental fluctuations on population growth and
ecosystem stability Chesson and Warner (1981); Ridolfi et al. (2011). On one hand,
environmental stochasticity enhances fluctuations and extinction
rates, that can have a destabilizing effect on the ecological
community. On the other hand, it can also foster stability, as the
temporal alternance of species can effectively reduce the strength of
interspecific competition.

Similarly to the case of spatial disorder, one can design
quasi-neutral models where habitat-preferences for different species
are time-dependent, i.e. where in each time window there is a
preference for a randomly chosen species. Different works have
recently analyzed this type of models, showing that time-dependent
habitat preference greatly improves predictions of empirical
ecological patterns with respect to purely neutral theories
Danino et al. (2016b); Shem-Tov et al. (2017); Hidalgo et al. (2017); Spanio et al. (2017). In particular, it has been
claimed that these models provides more realistic estimates of
dynamical quantities, such as average species persistence times and
distributions of species turnover Azaele et al. (2006), compared
with their neutral counterparts.

iv.6 Models with density dependence

In ecology, one speaks of density-dependence or Allee effect
when the fitness of an individual depends on the abundance of the
species it belongs to. The underlying mechanisms can be very diverse,
from cooperative defense/feeding to spreading of parasites among
conspecific. An interesting scenario is that of negative
density-dependence, i.e. when individuals belonging to more abundant
species have lower fitness. It is established that, in well mixed
systems, negative density-dependence significantly favors species
coexistence Chesson (2000). Versions of the voter model
implementing a negative density-dependence have been studied in the
literature Molofsky et al. (1999); Schweitzer and Behera (2009). In these models, the
reproduction probability of an individual depends on the number of
conspecific individuals in a given local neighborhood. Strictly
speaking, these models are not neutral: the neutral hypothesis is
defined at the level of individuals Hubbell (2001), and here
individuals belonging to species of different abundance do not have
the same fitness. However, these models, as the other models
considered in this Section, are still symmetric, since all species are
treated on equal footing. Interesting phenomena like the possibility
of spontaneous breakdown of such a symmetry –thus leading to
asymmetric species coexistence– have been recently uncovered at the
mean field level Borile et al. (2012).

The range of ecological problems discussed in this review is by
no means exhaustive, and we believe there are many directions that
still need to be explored or fully understood.

A prominent example is the role of different speciation mechanisms on
spatial biodiversity. In the models discussed in this review,
speciation events involve a single individual (point speciation
mode, in the language of evolutionary ecology). This assumption is
convenient from the modeling perspective, but leads to fitted values
of the speciation rate that are incompatible with independent
estimates Ricklefs (2006). This assumption also tends to
generate too many young species which last for a short time and
overweights rare species. To address these issues, recently, another
mechanism called protracted speciation has been proposed in the
context of neutral models Rosindell et al. (2010). In
protracted speciation, the speciation event does not occur at a single
generation, but is a gradual event lasting for some generations.
Introducing protracted speciation partially solves some of the
aforementioned problems Rosindell et al. (2010). In real
ecosystems, even more speciation mechanisms are at play
Coyne and Orr (2004). For example, in parapatric speciation, two
spatially-separated population of the same species can diverge and
give rise to two different species. This would correspond to a
speciation event involving a group of individuals
rather than a single one. The role of different speciation modes in
maintaining biodiversity and in patterning the spatial organization of
species is still under discussion and modeling results can provide
very useful contributions to this debate.

As mentioned in the Introduction, ecological neutral theory elicited a
heated debate which is far from being solved as, in many cases,
non-neutral models based on the concept of niche and neutral models
yield similar fits of biodiversity patterns
Chave et al. (2002); McGill (2003); Tilman (2004). In recent years a new
view on this debate has been emerging. In Chase and Leibold’s words:
“niche and neutral models are in reality two ends of a continuum with
the truth most likely in the middle”
Chase and Leibold (2003). Indeed, the ecological forces underlying
niche and neutral models are not mutually exclusive, and demographic
stochasticity plays an important role also in non-neutral
settings. However, it has been difficult to clarify the importance of
different neutral and non-neutral mechanisms, as most non-neutral
model are characterized by a large number of parameters. Some progress
in this direction has been obtained in simplified settings which,
similarly to the model presented in Sect. IV, allow
for a controlled departure from neutrality. For instance, Haegeman and
Loreau Haegeman and Loreau (2011) added the main ingredients of neutral
theory, demographic stochasticity and immigration, to a Lotka-Volterra
competition model. Similar problems have been studied in
Refs. Noble et al. (2011); Noble and Fagan (2011); Pigolotti and Cencini (2013). An interesting
future direction would be to study similar models in a spatial
context.

In many ecological communities, in particular of microbial organisms,
ecological and evolutionary timescales are not separated.
Eco-evolutionary models describing both processes are becoming more
and more important Villa Martin et al. (2016). Neutral theory has provided a simple
framework to describe patterns in these communities, for example in
gut microbiota Jeraldo et al. (2012). These systems call
for new theoretical efforts and new observables, such as
generalizations of the β-diversity taking into account genetic
differences among individuals Houchmandzadeh (2017).

We have seen throughout this review how some observables measured by
ecologists corresponds to well known quantities in statistical
physics: for instance, the β-diversity is closely related with a
two-point correlation function. Other observables, such as SARs and
SADs, are less common in statistical physics. A potentially fruitful
future direction is to consider other observables which are common in
statistical mechanics, such as multi-point correlation functions, and
measure them in ecosystems. In this direction, it is very interesting
the study of species clustering in Plotkin et al. (2002) based on
the theory of continuum percolation.

In summary, we presented an overview of different
stochastic spatial models in population ecology. We have seen that
even very simple models are a source of challenging problems in
statistical physics. In particular, because of speciation, each
species is bound to extinction and is therefore ultimately
transient. This feature is in contrast with traditional classical spin
system defined on a lattice where, even when in out-of-equilibrium
conditions, the number of spin components is fixed from the beginning.
Further, ecosystems are typically two-dimensional and, due to the
underlying diffusive behavior, D=2 is the critical dimension for
these models. We have shown that this fact often leads to logarithmic
corrections to scaling laws, which have been difficult to analyze both
analytically and numerically. Despite these difficulties, remarkable
progress has been made in recent years. We believe that
cross-fertilization between statistical physics and ecology will
be more and more important in the future to deepen our quantitative
understanding of how ecosystems are organized.

Appendix: General scaling relationships

In this brief Appendix, we discuss general condition imposed on the
functions f and g by the properties of the function Ψ,
depending on the exponent Δ, see eq.(13),
eq. (14) and Zillio et al. (2008).
Let us write the normalization condition for P(n;A)

∑nP(n;A)≈g(A)f(A)∫Λn0/f(A)dxx−Δ=1.

(30)

The infrared cutoff Λ is related to the fact that the function
ψ(x) is a power-law for small x only and rapidly decays for
larger arguments, see e.g. Fig. 5.
The integral is singular for small x and Δ>1 and thus

1∼g(A)f(A)f(A)Δ−1=g(A)f(A)Δ.

(31)

On the other hand, if Δ<1, the integral is weakly dependent on f(A),
so that

1∼g(A)f(A).

(32)

Similarly, the first moment of Ψ is

⟨n⟩∼g(A)f2(A)f(A)Δ−1=g(A)f(A)Δ+1

(33)

if 1<Δ<2 and

⟨n⟩∼g(A)f2(A)

(34)

for Δ>2. Combining the expressions above,
different regimes emerge as a function of Δ: if Δ<1, f(A)=⟨n⟩, while for 1<Δ<2, f(A)=⟨n⟩1/(2−Δ), while no specific prediction for f(A)
can be made in the case Δ≥2. In particular, for Δ<1 one has a simple scaling form f(A)=⟨n⟩ and
g(A)=1/⟨n⟩ which applies, for example, to the 1D
case as described in the main text. The marginal case Δ=1 is
treated in detail in Sec. II.5.

Acknowledgements.

MAM is grateful to the Spanish-MINECO for financial
support (under grant FIS2013-43201-P; FEDER funds), as well as to
J. Hidalgo, S. Suweis, A. Maritan, C. Borile for a long term
collaboration on topics related to the content of this paper.