Field theory for biogeography: a spatially explicit model for predicting patterns of biodiversity

Predicting the variation of biodiversity across the surface of the Earth is a fundamental issue in ecology, and in this article we focus on one of the most widely studied spatial biodiversity patterns: the species–area relationship (SAR). The SAR is a central tool in conservation, being used to predict species loss following global climate change, and is striking in its universality throughout different geographical regions and across the tree of life. In this article we draw upon the methods of quantum field theory and the foundation of neutral community ecology to derive the first spatially explicit neutral prediction for the SAR. We find that the SAR has three phases, with a power law increase at intermediate scales, consistent with decades of documented empirical patterns. Our model also provides a building block for incorporating non-neutral biological variation, with the potential to bridge the gap between neutral and niche-based approaches to community assembly. Ecology Letters (2010) 13: 87–95

This transformation means we can now take the limit as the separation between discrete spatial sites, ∆x and ∆y in the x and y directions respectively, goes to zero. To take this limit, we use dimensional analysis to assign the following scalings: θ → θ s ∆x∆y Q ij → Q(x i − x j , y i − y j )∆x∆y, so that the spatially-implicit neutral theory fundamental diversity parameter, θ, becomes a new fundamental diversity parameter per unit area. We have further replaced Q ij , the probability of dispersal from site i to site j, by a continuous function, known as a dispersal kernel, also with dimensions per unit area. The continuum limit of Eq. (1) is then the following functional differential equation: The functional derivatives denoted by δ/δH(x, y) can be thought of as a generalization to continuous space of the partial derivatives, ∂/∂h i , and are defined so that: The first term on the RHS of Eq.
(2) characterizes birth and dispersal, while the remaining terms represent mortality and speciation, respectively. The exponential terms in H(x, y) differentiate this from the canonical non-interacting quantum field theory (Zinn-Justin, 2002) and arise from the discrete nature of individuals.
So far we have allowed for all possible dispersal kernels, but we now model dispersal as a diffusion process. If the moments of the dispersal kernel Q(w, z) are finite, we may rewrite the dispersal term in Eq.
(2) as an expansion in spatial derivatives of δZ/δH, weighted by the moments of the dispersal kernel. To derive our central equation we keep only the first two terms in this expansion, corresponding to a diffusion approximation: where ∇ 2 = ∂/∂x 2 + ∂/∂y 2 , and the length-scale associated with the dispersal process is σ, equal to the square root of the variance of the underlying dispersal kernel. We expect the diffusion approximation to be valid beyond very small length-scales of order ∼ σ, so long as the moments of the dispersal kernel drop off sufficiently quickly, as discussed in (Chave & Leigh, 2002). Having made this diffusion approximation, we obtain the following functional differential equation for Z [H, t]: also given in the main text.

Equilibrium Partition Function and SAR
An equilibrium solution for the partition function is a time-independent functional, Z eq [H], which may be defined formally as the following expectation value: Z eq [H] then satisfies a local, time-independent version of Eq. (5), We can relate this equilibrium partition function to Φ(h, R), defined by Eq. 7 in the main text, by setting the function H(x, y) to be a constant h inside the circular sample region, and zero outside. Making this choice for H(x, y), Eq. (7) reduces to a differential equation: Rect(x, y, R) is a two-dimensional version of the rectangular function, so that it takes the value 1 within a circle of radius R and zero outside this radius, reflecting the extent of the sampled region. The other parameters in this equation are given by: This equation is useful because f (x, y, h, R) can be related to the generating function Φ(h, R) by integrating Eq. (9) over the sampled region: so that if we integrate the appropriate solution for f (x, y, h, R), we have a solution for ∂Φ/∂h as required.
The differential equation Eq. (8) can be solved by obtaining separate solutions inside and outside the sample region, and gluing them together by requiring continuity and continuity of first derivatives; the second derivative of f (x, y, h, R) is clearly discontinuous, owing to the form of Eq. (8). Given the symmetries of the problem we write f as a function of polar coordinates, r and φ and look for a solution independent of φ. Inside the sample region the solution takes the form: while outside the solution is where I 0 and K 0 are modified Bessel functions, and A and B are as yet undetermined coefficients. We fix the coefficients A and B by imposing that f in (R, h, R) = f out (R, h, R), and that the first derivatives of f in and f out are also equal at r = R. The relevant result is that Integrating this result from r = 0 to r = R, and from φ = 0 to φ = 2π, we obtain our solution for ∂Φ/∂h,

Taylor Series Approximation for ∂Φ/∂h
We need to integrate Eq.(15) with respect to h, from h = 0 to h = −∞ to obtain the SAR, as explained in the main text. This integral may be completed numerically, but in order to arrive at a closed form analytical result we Taylor-expanded the Bessel functions in Eq.(15) and integrated analytically. To compare with empirical data we could simply complete the numerical integral for appropriate parameter values of α and σ, but the Taylor-series result is an excellent approximation to the numerical integration, as we show in Figure 1 for α = 1, α = 10 −6 and α = 10 −10 :

Correlation Functions and Beta Diversity
Functional derivatives of the partition function are by definition equal to spatial correlation functions. We can therefore expand the partition function as a power series in H(x, y), and extract correlation functions from the coefficients. We focus here on the first two terms in the expansion, for the one-point and two-point correlation functions: Z[H, t] = dxdy n(x, y, t) H(x, y) + n(x 1 , y 1 , t)n(x 2 , y 2 , t) H(x 1 , y 1 )H(x 2 , y 2 ) + . . .
Using our fundamental master equation, Eq. (5), we can derive differential equations for these correlation functions, known as Schwinger-Dyson equations in quantum field theory. To obtain these equations we insert the expansion (16) into Eq. (5), take an appropriate number of functional derivatives with respect to H(x, y), and then set H(x, y) = 0 everywhere. For the one and two point correlation functions these equations are respectively: ∂ ∂t n(x 1 , y 1 , t)n(x 2 , y 2 , t) = ∇ 2 1 + ∇ 2 2 + 2b − 2d n(x 1 , y 1 , t)n(x 2 , y 2 , t) It is similarly straightforward to derive differential equations for higher order correlation functions, by expanding the partition function to higher order in H.
Because our spatial neutral theory is non-interacting, we found only single functional derivatives in Eq. (5), and as a consequence we can solve exactly for these correlation functions, order by order. The equilibrium solution for the one-point function, which is the expectation value for the abundance of a single species, is a constant density, given by: We should expect that this expectation value is constant across space, because there are no preferred spatial locations in our spatially-homogeneous system. This result also means that the expectation value for total density of individuals, summed across all extant species, is given by where the dimensionless parameter α = d b − 1.
The second coefficient in the expansion of Z[H, t] is the two-point correlation function, n(x 1 , y 1 , t)n(x 2 , y 2 , t) , which characterizes the spatial clustering of individuals belonging to a single species. The equilibrium solution of Eq. (18) with appropriate boundary conditions is: where r = (x 1 − x 2 ) 2 + (y 1 − y 2 ) 2 is the separation between the two points in space and K 0 is a modified Bessel function of the second kind. The two-point function (21) characterizes the spatial correlation of individuals belonging to a single species, but we can use our model to predict a community level measure of beta-diversity, the function F (r) introduced by Chave & Leigh (Chave & Leigh, 2002). F (r) is defined as the probability that two individuals picked at random from a community, but separated by distance r, belong to the same species.
To derive F (r), let us first focus on one particular species, and consider the probability that two randomly-sampled individuals belong to it. This is the expectation value of the fraction of individuals belonging to this species at one spatial location, multiplied by the fraction of individuals belonging to the same species at a second spatial location. If the variation in the total density of individuals across space is negligible, this expectation value is given by: n(x 1 , y 1 ) J n(x 2 , y 2 ) J .
To obtain F (r), we simply sum over all species, and use our assumption of neutrality. The result is which has the same functional form derived in (Chave & Leigh, 2002).