Cross-scale neutral ecology and the maintenance of biodiversity

One of the first successes of neutral ecology was to predict realistically-broad distributions of rare and abundant species. However, it has remained an outstanding theoretical challenge to describe how this distribution of abundances changes with spatial scale, and this gap has hampered attempts to use observed species abundances as a way to quantify what non-neutral processes are needed to fully explain observed patterns. To address this, we introduce a new formulation of spatial neutral biodiversity theory and derive analytical predictions for the way abundance distributions change with scale. For tropical forest data where neutrality has been extensively tested before now, we apply this approach and identify an incompatibility between neutral fits at regional and local scales. We use this approach derive a sharp quantification of what remains to be explained by non-neutral processes at the local scale, setting a quantitative target for more general models for the maintenance of biodiversity.


Backward Equation Derivation
We consider sessile individuals in d-dimensional space that give birth to conspecifics at a rate (b − ν) and die at rate b. Offspring are distributed relative to their parents according to a dispersal kernel B, i.e. the probability density that an offspring of an individual at position r is located at a position r is B(r − r ). We assume that there are no density dependent interactions, so that there is no zero-sum rule and that species (and, more generally, lineages) do not interact. The difference between the birth and death rates represents the production of offspring that speciated (to a novel species) at birth.
We define P (k, A, r, t) as the probability that a single individual at position r at time t 0 has exactly k conspecific descendants in a region A at time t + t 0 (the dynamics are translationally invariant in time, so this probability does not depend on t 0 ). We consider a single individual at time time 0, and enumerate the possibilities for its lineage during a small ensuing interval ∆t, and the consequent value of P (k, A, r, t):

Process
Probability Probability of k descendants in A at time t Nothing 1 − (2b − ν)∆t P (k, A, r, t − ∆t) Death b∆t δ k,0 Birth, offspring at r (b − ν)∆t B(r − r )dr k m=0 P (m, A, r, t − ∆t)P (k − m, A, r, t − ∆t) If the individual dies, then there will be zero descendants (δ ij [= 1 if i = j, 0 otherwise] is the Kroneker delta). If the individual gives birth, then there are two statistically independent lineages, and the total number of descendants in A is the sum of the number of descendants from the two lineages.
We assume that the dispersal kernel B decays to zero rapidly enough that P (k − m, A, r, t) varies much more slowly in space, and perform a Taylor expansion for P (m, A, r , t) so that P (m, A, r , t)B(r − r )dr = B(r − r) 1 + (r − r) · ∇ + ((r − r) · ∇) 2 2 + . . . P (m, A, r, t)dr ≈ P (m, A, r, t) + K∇ 2 P (m, A, r, t), where K = 1 2 (r − r) 2 B(r − r)dr and we have assumed that B is normalised so that B(r − r)dr = 1, and has rotational symmetry so that (r − r)B(r − r)dr = 0. The backward equation can then be written as This equation is to be solved for t > 0, with initial condition P (k, A, r, t) = δ k,1 if r is inside the region A, and zero otherwise, representing the fact that if t = 0 the lineage consists of the founding individual only.
Finally, we define the generating function of P as: so that where σ 2 is twice the variance of the dispersal kernel. The boundary condition is that G(z, r, 0, A) is = z if r is within the sampling region defined by A (e.g. a line segment of length L in the one dimensional problem, and an area of a particular geometry in the 2-dimensional case).
A similar backward equation can be derived for species that, rather than being sessile and dispersing at birth, performed random walks. This equation turns out to be identical to eqn. (5) except that the factor G in front of the Laplacian is replaced by 1. As discussed later, this factor does not affect the Species Area Curve or Species Abundance Distribution in the biologically relevant limit ν b.

Species Area Curve
We now consider an approximation method to find solutions of Eq. (5). Considering spatial dimension 1, first we define: so that φ satisfies: with an initial condition φ(x, t = 0, L) = R(x, L), where R(x, L) is a rectangular function, equal to zero for x < −L/2 and x > L/2, and equal to one for −L/2 < x < L/2.
The function φ(x, t) is the probability that an individual appearing in a speciation event at a time t in the past, and at location x, will have one or more descendants in the focal region between −L/2 and +L/2 in the present day. In order to derive the Species-Area relationship from this probability distribution (where 'area' indicates the one-dimensional length of the focal region, L), we need to integrate over all speciation events, which in the neutral model occur at a rate νρ per unit time per unit area, where ρ is the equilibrium density of individuals in space. Hence, our goal is to derive a solution for: Due to the nonlinearity in Eq. (7), solving for φ(x, t, L) exactly does not seem tractable. But the initial and final conditions for φ suggest a linear approximation, if we treat φ 2 (x, t, L) φ(x, t, L)R(x, L). While true at t = 0, and true at late times when φ → 0, this does not hold for general t, and with this approximation for S(L) we would underestimate the number of species at large values of L.
The problem is clear-at intermediate times, φ(x, t, L) will be non-zero outside of the focal region, and will interpolate between one and zero in side the focal region.
We therefore handle this discrepancy by approximating φ 2 (x, t, L) as φ 2 (x, t, L) = φ(x, t, L) while x is within the focal region, and outside of the focal region we set φ 2 (x, t, L) ∝ φ(x, t, L) with a (we expect small) constant of proportionality to be determined. In addition, we replace the factor (1 − φ) in front of the Laplacian by 1, because when ν b the dominant contribution to the integral that determined S comes from large values of t, by which time φ 1. This leads to an equation of the form: So in fact, we have two equations to solve, as the approximation we have used leads to different equations inside and outside of the focal region defined by L: We can now integrate over time, before solving these equations as a function of space, to obtain: , respectively. Note that we have also introduced a new, effective rate of introduction of new species per unit space and time, bν eff ρ, instead of νρ, for consistency at small values of L with the term ν eff φ out in Eq. (11). It may seem like we have introduced a free parameter or parameters by allowing for ν eff , but in fact this effective rate is fixed by the large scale behavior, i.e. as L → ∞. In this limit, which leads to a solution So in order to match the standard neutral result for a well-mixed community, at large scales we have that with no free parameters.
We can solve the pair of equations (11) by imposing that there is no singular behavior, that Φ out asymptotes to zero for large x, and that at the boundaries ±L/2 both Φ out and Φ in and their first derivatives match. The result is: where we have defined the inverse length-scale m = b/D = 1/σ which is related to the standard deviation of the dispersal distance. We now integrate over space to obtain our approximate prediction for the one-dimensional Species-Area relationship.
In the following subsections we will show that despite our approximation, this provides an extremely accurate prediction of the relationship across a broad range of areas.
We now provide the corresponding result in two spatial dimensions. We apply exactly the same approximation, but where now we interpret R(x, y, L) as the 'top-hat' function, which is equal to one inside a circular region of radius L, and is equal to zero outside. We then solve for functions inside and outside of this circular region: This leads to solutions which depend only on a radial coordinate, r (distance from the origin), and not on the corresponding polar coordinate: and integrating over all space we find: where again ν eff = − ν 1−ν log(ν/b). In terms of sample area A = πL 2 we can rewrite this as: In the main text we replaced the inverse length-scale m with a length-scale σ = 1/m, but both can be related directly to the parameter D we introduced in the formulation of this problem.

Species Abundance Distribution
We now consider the same kind of approximation method to find solutions of Eq. (5), but instead of considering just total species richness, we define (again first considering the one-dimensional case): so that φ again satisfies: with an initial condition φ(x, z, t = 0, is the same rectangular function as above. To obtain the Species Abundance Distribution, we need to integrate over all speciation events, which in the neutral model occur at a rate νρ per unit time per unit area, where ρ is the equilibrium density of individuals in space. Hence, our goal is to derive a solution for: where using this definition, Ψ(z, L) is related to the Species Abundance Distribution in a sample taken from the region of size L by We now extend our previous approximation for the species area curve. We follow the same principle to set φ 2 (x, z, t, L) (1 − z)φ(x, z, t, L) when x is within the focal region, while outside of the focal region, we set φ 2 (x, z, t, L) = g(z)φ(x, z, t, L) for a function g(z) to be determined by the requirement that we match the known behavior at large values of L. We again replace the factor (1 − φ) in front of the Laplacian by 1, as this makes a vanishing difference when ν → 0. This leads to an equation of the form: Again, we have two equations to solve, as this approximation leads to different equations inside and outside of the focal region defined by L: We can now integrate over time, before solving these equations as a function of space, to obtain: , respectively. Note that we have also introduced a new, effective rate of introduction of new species per unit space and time, bg(z)ρ, instead of νρ, for consistency at small values of L with the term bg(z)φ out in Eq. (27). The function g(z) is then fixed by the large scale behavior. In this limit, which leads to a large-scale solution So in order to match the standard neutral result for a well-mixed, non-zero-sum community [1], which is we need to set We note also that g(0) = ν eff , and so this approximation simply reduces to the approximation we used to derive the Species-area Curve above; our solution above for φ(z, t, L) is equal to φ(x, z = 0, t, L) here, so that our approximation for the Species-area curve satisfies these same equations but with z = 0 (as it should do).
The result of solving this pair of equations is essentially the same as above. The dependence on z of the parameter does not affect the solution as a function of x, it just changes the parametrization in that solution: where again m = b/D = 1/σ in the main text, and for ease of notation we introduce as in the main text, and such that g(z), f (z) and h(z) are related by g(z) = h(z)f (z)/(1 − z). We now integrate over space to obtain our approximate prediction for the one-dimensional Species Abundance Distribution.
The same approach in 2 spatial dimensions for a circular region of radius L and area A = πL 2 leads to: and

Universal behaviour as ν → 0
Here, we show that the species abundance distributions given by our model exhibit the same universality property found in simulations by Rosindell and Cornell [2] i.e. that the species abundance distributions form a family of curves, parametrised by the single paramemeter Aν bσ 2 .
First, we show that the exact solution to the backward equation for φ has this scaling property. If we define then eqn. ((7)) becomes and the initial condition becomes Q(X, . Therefore, Q only depends on the parameters through the combinations i.e. Q =Q(X, T, Z, Y ) for some functionQ. The generating function for the abundance distribution in 2D is then given by (see eqn. ((23))) for some functionΨ.
While we do not have an expression for the exact solution Ψ, we can verify that our approximate solution Ψ 2d has the same scaling behaviour.
which is of the same form as eqn. ((37)).
We have thus shown that the generating function of the abundance distribution is a function of two parameter combinations only. We will now show that this is equivalent to the observation that the species abundance distribution is a one-parameter family of curves [2]: S(k, A) = νS νk, Aν bσ 2 .
(note that the expression in ref. [2] describes the scaling of logarithmic Preston classes of abundance, and hence is missing the prefactor ν). To see how this is related to our scaling expression for Ψ(z, A), we write where m = νk. We need to proceed with caution in caseS has a non-integrable singlarity in its first argument. Abundance in Preston classes of low order appears from Fig 2 in the paper to approach a finite limit, so we assume thatS ∼ s(Y ) νk at small (νk). Without loss of generality, we write dy is the exponential integral and f is a (finite) function of two arguments. The component of this expression that depends on Z takes the same scaling form as found above for the backward equation model described above. The term that is independent of Z does not contribute to any particular S(k, A), but does contribute to the total species richness, and indeed we can identify

Biological Interpretation of the Approximation Method
By approximating the non-linear term in the defining backward Equation (7) by a heterogeneous linear term, we found the pair of equations (11) to solve for the species area relationship: (For simplicity we work in one spatial dimension but the interpretations are identical in 2d.) We could equally well interpret these not just as an approximation to Eq. (7), but as a biological model in their own right. If we do so, can we reinterpret these equations and understand biologically why this linear approximation works? Eqs. (11) constitute a system where there is only mortality (driving loss of species from the focal, sample area), dispersal, and input from speciation. This might be expected, since species are only removed from the focal region when there is a mortality event, and only added when there is a speciation event landing in the focal region, or dispersal in from outside. However, in these equations the effective rates of species loss are different for species which originated outside the focal region (rate ν eff ), versus those that originated inside the focal region (rate b), and that is what we must explain.
Remembering that in the original derivation above of Eq. (7), the per capita birth rate was b − ν and mortality rate was b, this interpretation of rate of species loss in the equation for Φ in becomes clear: species are lost from the focal region at the rate at which a single individual dies. I.e. we are approximating that species which originate within −L/2 < x < L/2, will only not be found at time t in this sample region if it goes extinct, and this rate of loss is approximated by the rate of loss b of a single individual. For species outside the sample region, the rate of loss from mortality is bν eff = bν b−ν log(b/ν). What is this number? In fact, it is equal to b/ n , where n is the expected population size of an extant species (i.e. total number of individuals divided by total number of extant species). On average, a species originating outside the focal region at any point in the past, is lost from the focal region at a effective rate, b/ n .

Comparison with Field Theory/Forward-in-time Equations
In an earlier paper, one of the authors derived a forward-in-time approach to these same spatial neutral models [3]. In that approach, we also began with the case of a spatially-discrete landscape. Here we will recap the basic features and approximation we made in that paper, and where they break down relative to our current approach. We work here with individuals that diffuse across the landscape, rather than sessile individuals, but as discussed above this does not change the patterns of abundance and occupancy when the speciation rate is small.
We first describe the state of the discrete system using the probability distribution P (. . . , n i . . . , t) that there are n i individuals at each spatial location, i, at time t, belonging to a focal species. Individuals in this spatially-discrete model die with a per capita mortality rate, d, produce new offspring at a per capita birth rate, b, and may transfer to nearest neighbour cells at a rateD. In addition, there is a speciation process modeled as immigration from outside the system at a ratẽ k from 0 to 1 individual: We could also remove this last term, introducing speciation, and thus allow each species to reach permanent extinction. We would then sum the contributions to the present day state from all species that originated at some point in the past, assuming a uniform speciation rates across time and space. Before taking the limit of continuous space, we rewrite the dynamics of our discrete community in terms of a moment generating function. This generating function is defined by a sum over all spatial configurations of individuals: Rewriting Eq.(38) in terms of this generating function, we find a new defining equation:

Taking a Continuum Limit
We denote the lattice spacing by ∆, and define the continuum limit as follows: Finally, to define the continuum limit for the sum over nearest neighbours, we consider a square, d-dimensional lattice: With these identifications, the multivariate generating function Eq. (39) becomes a functional of the source H(x): where n(x) , n(x 1 )n(x 2 ) etc are expectation values of the number densities and correlations of individuals as a function of spatial location (express this better). This generating functional satisfies the continuum limit of Eq. (40), the following functional differential equation: where we have performed an integration by parts and assumed that the source H(x) vanishes at infinity. Finally, we make a change of variables for the source,

Equal Time Correlation Functions
The n-point spatial correlation functions for this model, taken at equal times, t, satisfy a set of partial differential equations. These equations are obtained by expanding Z[J(x), t] as a functional Taylor series: The coefficients of this Taylor series are obtained by taking functional derivatives of Eq.(46) with respect to J, and then setting J = 0. For the first two orders we have: Each successive order relies only on solutions for correlation functions of lower order, and so the system of linear partial differential equations can be solved exactly, given a set of initial data.

Species Area Relationship
We now consider the time-independent probability P (N, L) that at late times there are N individuals in a given sample region extending from −L to +L. The generating function of this probability is: To underline the interpretation: P (N, L) is the (assumed time-independent solution for the) probability distribution that we will find N individuals in the region between −L and +L at late times. The second equality arises because this generating function can be obtained by setting J(x) = jRect(x, L) in the late-time, time-independent solution for Z[J, t], where Rect(x, L) is the rectangular function in 1d. I.e. we will set J(x) zero outside the sample region and equal to j inside. The expected number of species in the sampling region defined by L is proportional to 1 − P (0, L), and so this is the quantity we are aiming to solve for. If we can solve for this generating function then we have P (0, L) = ψ(−1, L). This is also known as the empty interval function, e.g. [4]. To find ψ(−1, L), we next define the modified moments which have insertions of e log(1+j) L n(x) compared with the usual moments: We note that the reason for using these functions is that: and hence if we can solve for f 1 (x, k, L) we will have the empty interval function, P (0, L).
We can obtain differential equations for these modified moments by taking successive functional derivatives of Eq.(46): etc. So setting J(x) = jRect(x, L) and time derivatives equal to zero in these equations we have a kind of moment hierarchy: So far there is no approximation. In our earlier paper we truncated and solved the first equation in this hierarchy: While giving qualitatively accurate description of the shape of the SAR, this earlier approximation becomes quantitatively inaccurate as a speciation rate becomes small. It also fails to give a good description of the Species Abundance Distribution.

Numerical Inversion of the SAD Generating Function
We implemented a method described in [5] for numerical contour integration using Cauchy's theorem. Our annotated and documented code will be freely available on the O'Dwyer lab GitHub repository.

Uncertainty in Estimating Speciation Rates
In our analysis of empirical data showin in Fig.3 of the main text, we draw from earlier work fitting speciation rate ν and dispersal length-scale σ using [6] using multiple 1ha plots distributed across Panama. One other important conclusion from this work was that the best fit parameter for ν was highly variable across different locations (e.g. in comparison with Panama, two South American locations were fitted with speciation rates that were orders of magnitude smaller). Here, we repeat the analysis of Figure 3, but using one of these lower fitted speciation rates (and the corresponding value of σ, alongside the local density of individuals in Ecuador). Clearly, this cannot be a complete analysis, as we are no longer providing the best fit to the large-scale Panamanian data. But it is interesting to ask whether these much lower speciation rates can significantly change either the local prediction for overall species richness, or whether the shape of the local SAD becomes more realistic. We show this analysis in Fig S1. The result is that we do find a larger expected species richness locally, but still not close to the observed number of species at BCI. At the same time, the shape of the species abundance distribution is even more skewed away from rare species. Comparison of Observed and Predicted Abundances Figure S1 Neutral predictions at BCI, using Yasuni fitted parameters. Neutral predictions are generated by fitting our spatial neutral model using large-scale data reported and analyzed in [6], but now using the lower reported speciation rate fitted using data surrounding the Yasuni CTFS plot in Ecuador (ν = 1.7.10 −11 , rather than ν = 5.10 8 ), alongside the other parameters σ and density ρ fitted at Yasuni (which were of the same order of magnitude as those in Panama). The results show that these large-scale fits still produce a local-scale prediction for species abundances that both underestimates local species richness, compared with observed data, and still skews abundances from rare to more abundant species.