Spatial patterns emerging from a stochastic process near criticality

There is mounting empirical evidence that many communities of living organisms display key features which closely resemble those of physical systems at criticality. We here introduce a minimal model framework for the dynamics of a community of individuals which undergoes local birth-death, immigration and local jumps on a regular lattice. We study its properties when the system is close to its critical point. Even if this model violates detailed balance, within a physically relevant regime dominated by fluctuations, it is possible to calculate analytically the probability density function of the number of individuals living in a given volume, which captures the close-to-critical behavior of the community across spatial scales. We find that the resulting distribution satisfies an equation where spatial effects are encoded in appropriate functions of space, which we calculate explicitly. The validity of the analytical formul{\ae} is confirmed by simulations in the expected regimes. We finally discuss how this model in the critical-like regime is in agreement with several biodiversity patterns observed in tropical rain forests.


I. INTRODUCTION
Several authors have showed that the parameters of models which describe biological systems are not located at random in their parameter space, but are preferably poised in the vicinity of a point or surface which separates regimes of qualitatively different behaviors [1]. In this sense, stationary states of living systems are not only far from equilibrium, but bring the hallmark of criticality. Although the connection between underpinning dynamics and measurable quantities is sometimes tenuous and hence conclusions about criticality loose, empirical examples span a wide range of biological organization, from gene expression in macrophage dynamics [2], to cell growth [3], relatively small networks of neurons [4], flocks of birds [5] and, possibly, tree populations in tropical forests [6].
In this article, we mainly focus on the spatial patterns emerging from a minimal model of population dynamics close to its critical point. This latter is identified as a singularity in the population size of the system, similarly to what happens in the theory of branching processes in the sub-critical regime [7], when the fluctuations play a crucial role. Therefore, in our model the critical point does not mark a transition between ordered and disordered phases sensu equilibrium statistical mechanics [8], although connections in a broader context may certainly exist. The emergent patterns are not calculated by using classical size-expansion methods, but introducing a parameter expansion which appropriately identifies criticality in the parameter space of the model.
The calculation of the probability distribution of large scale configurations emerging from the microscopic dynamics is challenging, even at stationarity [9,10]. When stochastic processes violate detailed balance, they have a generator which is not self-adjoint [11] and different states are coupled by probability currents at microscopic level [12]. These flows of probability among microstates break detailed balance, time symmetry and produce macroscopic non-equilibrium behavior. A common way to overcome these hurdles is to formulate some kind of effective Langevin equation which describes the dynamics of the mesoscopic variables of interest, loosing track, however, of the underlying microscopic dynamics [13][14][15]. Nonetheless, in this paper we study a model which, despite violating microscopic detailed balance [16,17], allows one to study analytically (stationary) out-of-equilibrium properties of spatial patterns. These latter emerge mainly because of the large intrinsic fluctuations of the local population sizes. Also, the model's mathematical amenability allows us to analyze in detail those spatial ecological patterns and to compare them with observation data for ecosystems with large species richness. The agreement between model predictions and empirical data highlights the usefulness of the approach and strengthens the connections between physics and theoretical ecology.
In this spatial metacommunity model, local communities (or, equivalently, sites or voxels) are located on a ddimensional regular graph (or lattice) where individuals are treated as well-mixed particles which undergo a birth and death process with local diffusion and constant colonization. We thoroughly investigate the spatial stochastic dynamics close to criticality and deduce the equation governing the evolution of the conditional distribution p(N |V ; t), the probability to find N individuals in a volume V at time t (in dimension d). Within this regime we map the equation of p(N |V ) of the out-of-equilibrium spatial model into an equation of a suitable stochastic process, which instead satisfies detailed balance. This model is described by functions of space, which we are able to calculate exactly. The exact stochastic simulations are always matched by our analytical formulae in the expected regimes.
The rest of the paper is organized as follows: in section arXiv:1907.08852v2 [physics.bio-ph] 6 Dec 2019 II we introduce the master equation of the model; in section III we calculate the mean and pairwise correlation; in section IV we study the generating function of the conditional probability density function of N (V ); in section V and VI we calculate the population variance and the conditional pdf of N (V ), respectively, along with a comparison between simulations and analytical formulae; in section VII we present an ecological application of the model; finally, section VIII includes some discussions and perspectives about the results.

II. MASTER EQUATION OF THE MODEL
This is a metapopulation model in which individuals live in local communities (or sites) located on a ddimensional lattice, L, whose linear side is a. If X i , i ∈ L, indicates an individual living in site i, the reactions defining the model's dynamics can be cast into the form where j indicates a site which is a nearest neighbor of the site i. In this model, individuals within local communities (or, equivalently, sites or voxels) are treated as diluted, well-mixed point-like particles which undergo a minimal stochastic demographic dynamics: each individual may die at a constant death rate r and give birth to an offspring at a constant rate b. The newborn individual remains in the same community with probability γ, whereas it may hop onto one of the 2d nearest neighbours with probability 1 − γ. Also, all communities are colonized by external individuals at a constant immigration rate b 0 , which prevents the system to end up into an absorbing state without individuals [18,19]. Notice that (for 0 ≤ γ < 1) spatial movement is always coupled to birth, so that only newborn individuals can move. This is because we have in mind an application to spatial ecology, where this model mimics the population dynamics in species-rich communities of trees, in which only seeds can move. However, it can be easily modified to include random walk behaviour or different dispersals -like those for bacteria or humans -in which the local birth and the hopping probability are in general decoupled. We have indeed verified that the generality of our final results does not depend on that coupling (see Appendix F).
In the language of chemical reaction kinetics, the first reaction represents an autocatalytic production; this and the hopping move are responsible for the break of detailed balance as shown in Appendix A. Therefore, stationary states of this process are non-equilibrium steady states, albeit the model is defined by linear birth and death rates. Indeed, it is worth emphasizing that each "patch" has not a maximum number of individuals which it can accommodate, but any population size is allowed, albeit large sizes have an exponentially small probability to occur. This is because the model has no intrinsic carrying capacity which leads the population to saturation. A carrying capacity has the advantage to bring in more realistic features, but it also makes the model mathematically more complicated because of non-linear terms.
Here we show that linear rates of a stochastic model per se can produce a relevant phenomenology within a stationary out-of-equilibrium model. Thus, since we wanted to focus on the regime near criticality, we have preferred to delve into a relatively simpler system, in which nonlinearities are neglected in a first approximation. More complicated dynamics are definitely important and will be studied in future works.
Finally, for the system to avert demographic explosion we have to assume that b 0 > 0 and 0 < b < r, but it turns out that the most interesting features emerge when b r, i.e., close to its critical point. Indeed, as we will show in Section VII, for comparable birth and death rates the model is able to describe several spatial patterns of tree species in tropical forests [20].
Let us now indicate with n i the number of individuals in site i. Assuming that within every site the spatial structure can be neglected and that we have perfect mixing, when the configuration of the system is {n} = {n i : i ∈ L}, the linear birth and death rates in site i, i.e., W + i ({n}) and W − i ({n}), read respectively Let P ({n}, t) be the probability to find the system in the configuration {n} at time t. Then the master equation for P ({n}, t) reads where the dots represent that all other occupation numbers remain as in {n} and it is intended that P (·) = 0 whenever any of the entrances is negative. The spatial generating function of the process is defined as FIG. 1: Illustration of the model. Individual trees are represented by dark green circles within local communities which are located on a regular graph (lattice). Each individual may die or give birth to an offspring with constant per capita rates. New individuals may either remain in the community where they were born with probability γ, or hop onto one of the 2d nearest neighbours with probability 1 − γ. Finally, all communities are colonized by external individuals at a constant immigration rate b 0 . The dynamics of the model is therefore defined by the jump rates defined in eq.(1).
where H i ≤ 0 for every i ∈ L. Multiplying both sides of eq.(2) by e k ∈L n k H k and summing over all configurations of the system, we find that ζ({H}, t) satisfies the equation This is the main equation of the model from which we will calculate the most important results. We were not able to find the full solution of this equation. However, one can gain a lot of information about the general properties of the process by looking into the probability distribution of the random variable N (V, t) = i∈V n i (t), where V is the set of sites in a d-dim volume. Before studying such a distribution, it is useful to calculate the mean number of individuals and the spatial correlation between any pair of sites.

III. MEAN AND PAIR CORRELATION
The equation for the mean number of individuals in the site k can be obtained by taking the partial derivative of both sides of eq.(4) with respect to H k and then setting {H} = 0: where µ := r − b and ∆ k is the discrete Laplace operator, which is defined as This finite difference equation can be solved in full generality and at stationarity, i.e. for t → ∞, we get simply n = b0 µ , regardless of any spatial location. The pairwise spatial correlation between sites k and l, i.e., n k n l , can also be obtained by taking the partial derivatives of both sides of eq.(4) with respect to H k and H l , and then setting {H} = 0 (see Appendix B for details): ∂ ∂t n k n l =D ∆ k n k n l + ∆ l n k n l + −2µ n k n l + 2b 0 n + (5) +δ k,l 2σ 2 n + b 0 + D∆ k n k where we have used the notation We also introduce which are dimensionless parameters and provide important information about how spatial diffusion intermingles with demographic dynamics. In order to solve eq.(5), we introduce a d-dim system of Cartesian coordinates where the coordinates of each site are given as a multiple of the lattice side a. Thus, a vector k indicates the corresponding position of a site. In this way, we can calculate the stationary solution of eq.(5) by writing n k n l as a Fourier series expansion. Exploiting translation invariance the final expression of the solution reads (see Appendix B) where p i is the i-th Cartesian component of the d-dim vector p and C is the hypercubic primitive unit cell of size 2π/a. Interestingly, pairs of sites de-correlate for γ = 1 or b = 0, when at stationarity we obtain n k n l − n 2 = cδ k,l , where c is a constant depending only on the demographic parameters and δ k,l is a Kronecker delta. However, c = n , pointing out that local fluctuations are non-Poissonian. Eq.(6) is amenable to a continuous spatial limit, obtained as a → 0, and provides a closed analytic form for the pair correlation. Indicating now with n(x) and n(y) the density of individuals on the sites located at x and y, respectively, in continuous space (and rescaling parameters accordingly), we find (see Appendix B) where |x − y| is the distance between the sites located at x and y, K ν is the modified Bessel function of the second kind of order ν [21], and we have defined whereD := Da 2 andb 0 := b 0 /a d for a → 0, and they are assumed to be finite in the limit. The first one is a standard scaling for spatial diffusivity, whereas the second scaling assumption comes from the requirement that spatial continuous constants are finite and non-trivial as a → 0 for x = y.
As the asymptotic behavior of K ν as z → ∞ is 2z ,λ is the correlation length of the system;ρ 2 has the dimensions of a d-dim volume and gives the local intensity of the correlations. Because eq.(7) is the continuum limit of eq.(6), this expression of the pair correlation function is a good approximation of the one in eq.(6) only when |x − y| a andλ a.

IV. GENERATING FUNCTION CLOSE TO THE CRITICAL POINT
In this section we introduce the parameters which allow us to identify a region close to the critical point of the model. This suggests an expansion which will lead to simplified equations which, nonetheless, carry a lot of information about the model. A simple way to make progress with the master equation in eq.(2) is the use of a formal Kramers-Moyal expansion [22]. It is well-known that there are limitations to this procedure and it has been criticized, because one cannot always pinpoint a small parameter for the correct expansion [22,23]. The system-size expansion solves these difficulties, but it has to be applied when the size of the system becomes large [23]. Here, however, it is not entirely evident what parameter should identify the size (population size or volume) of the system in the critical regime. Indeed, the model has no carrying capacity or maximum population size, and the total volume of the system could only provide us with the macroscopic equation, which has no interest in the present case.
In order to make analytical progress, we have introduced two dimensionless parameters, ε and η, which identify a non-trivial region when 0 < b < r, but b → r. This choice comes from the observation that communities of living organisms often appear to have very large demographic fluctuations. Several studies (see [6,24,25]) have pointed out that per capita birth and death rates are close to each other in seemingly different systems, thus suggesting that a sensible theoretical limit to study is when b approaches r. For example, a fit of eq.(7) to the empirical two-point correlation function of the tropical forest inventory of Pasoh natural reserve in Malaysia leads to empirical values of 2(r−b) r+b of the order of 10 −7 (see Fig. (4) and section VII for more details). We hence define ε := 2(r−b) r+b with the condition that b0 µ ε = O(1) as ε → 0 + ; in this way we obtain a constant ρ 2 = µ/b 0 ε, which in real systems is large because usually r/b 0 1. The parameter ε indicates how close the system is to the critical point, regardless of spatial diffusion. With the independent parameter η := D σ 2 , instead, we compare the importance of spatial diffusion with respect to demographic fluctuations. We will assume that η = O(ε) as ε → 0 + and hence η/ε = λ 2 , a new independent constant. This scaling assumption reflects that we want to explicitly analyze spatial effects in the critical limit. In fact, it can be easily proved that when η/ε → 0, the model is equivalent to the mean field model without spatial diffusion, whereas for η/ε → ∞ spatial diffusion dominates over the birth-death dynamics. When b and r are close to each other, we expect that population sizes can be well approximated with continuous random variables in each site. In order to understand when this is possible in relation to the parameters ε and η, we assumed that the generating function ζ({H}, t) is analytic at H i = 0 for any i and the most important contribution to the equation of ζ({H}, t) comes from a negative neighborhood of the origin with thickness O(ε). This is tantamount to introducing the change of variable H i = εS i into eq.(4) and to expanding in powers of ε, assuming S i = O(1) and S i ≤ 0. Up to linear order in η and ε we obtain where we have introduced the dimensionless time T := µt and now ζ, with a slight abuse of notation, indicates the generating function corresponding to continuous (and dimensionless) random variables. Therefore, the evolution equation for the generating function of the population sizes becomes where H i is the variable conjugated to the continuous random variable n i . The parameters also correspond to this continuum limit and now the population sizes n i have an exponential cut-off with a (large) characteristic scale given by σ 2 /µ. Eq.(10) leads to the following Fokker-Planck equation: where now {n} are continuous random variables and P ({n}, T ) is the corresponding probability density function. It is interesting to note that this is not exactly equivalent to the naïve Kramers-Moyal expansion of eq.(2), which would entail additional terms in the diffusive part. Nevertheless it is a diffusive approximation of the process in the regime identified by the parameters η and ε.

A. Conditional generating function
Eq.(11) cannot be solved in full generality, but we can better understand the underlying dynamics by studying the distribution of the population sizes in arbitrary volumes of space. Let us indicate with V the set of sites belonging to a d-dim volume |V| = V , and let us introduce the random variable N (V, t) = i∈V n i (t), i.e. the total number of individuals in V at time t. By indicating with P (N |V, t) the corresponding probability density function of N , we define the conditional generating function where h ≤ 0. We obtain the corresponding equation for Z by specifying H i , i.e.
and substituting this into eq.(10). Thus, This equation depends on f (i, h, V, T ) which in general is unknown. An equation for f can be derived by differentiating both sides of eq.(10) with respect to H k . Assuming as before that k ∈ V, we finally obtain the following equation In eq.(15) the function g is still unknown, but it is possible to show that ∆ i g(i, k, h, T ) = ∆ k g(i, k, h, T ) to leading order as a → 0 (see Appendix C). This allows us to obtain a closed equation for f . Indeed, we have after using the identity (16).
The spatial continuous limit of eq.(13) as a → 0 reads where in n we have substituted b 0 withb 0 . This equation is of crucial importance in what follows. Similarly, eq. (15) becomes where with the continuous coordinates (n i → n(x) and When there are no spatial effects, i.e.,λ = 0, eq.(17) reads which has the following solution at stationarity where h ≤ 0 and µ n V /σ 2 = V /ρ 2 . This is the generating function of a gamma distribution [23] with mean n V and variance n 2 Vρ 2 . Also, it is not difficult to verify that whenλ = 0, the solution of eq. (18) is given by In the next section we calculate the stationary variance of the population sizes in a volume V , i.e., of the random variable N (V ).

V. POPULATION VARIANCE
In this section we calculate the stationary variance of the population size, i.e., N (V ) = V dx n(x), where V is a finite volume. This quantity is of pivotal importance in the approximation that we are going to develop in the following sections and it is key when introducing spatial information in the equations. Thus we outline here the main calculations, leaving to Appendix D further details.
We start from the spatial continuous limit of eq.(5) as a → 0, which at stationarity reads  7).
where δ(x − y) is a Dirac delta and we have included only the leading terms as ε → 0. By integrating both sides of eq.(22) with respect to y ∈ V, we obtain an equation for n(x)N (V ) . We plan to find an explicit expression for this quantity, because it plays a key role in the following when we will calculate Z(h|V ). Henceforth, we will take V to be a d-dim ball of radius R and we will assume that the origin of the Cartesian coordinates is at its center.
Thus we indicate with |x| the distance from the origin of the site located at x in this coordinate system. With this notation and using the symmetry of n(x)n(y) with respect to x and y, we obtain at stationaritȳ This linear ODE has to be solved separately for |x| > R and |x| < R. The continuity of n(x)N (R) and its first derivative at the boundary |x| = R provide the solution.
The final results is (see Appendix D) where the function Ψ takes the following form for |x| ≤ R where I ν (z) and K ν (z) are the modified Bessel functions of the first and second kind, respectively [21]. Integrating both sides of eq.(24) with respect to x ∈ V, we obtain the equation for the second moment of N (R), i.e.
where ψ(R/λ) takes the following explicit form in dimension d (see Appendix D for its behavior): These two functions, namely n(x)N (R) in eq. (24) and N (R) 2 in eq.(26) will be used in the next sections to calculate a first order approximation of the spatially explicit probability density function of N (V ) in the vicinity of the critical point. Finally, these solutions allow us to write down the analytic form of the variance of N (R) in dimension d, i.e.

A. The spatial Taylor's law
Taylor's law was first observed in ecological communities [26,27], where natural populations show some degree of spatial aggregation. This was phenomenologically captured by assuming a scaling relationship between the variance and mean of population sizes in different areas. More generally, and recently, Taylor's law denotes any power relation between the variance and the mean of random variables in complex systems [28,29]. The law postulates a relation of the following form where C is a positive constant and α typically assumes values between one and two [27,29]. The spatial model which we have introduced can predict the behavior of this relation across scales without making specific assumptions. If we focus on the two dimensional case and fixλ, we obtain Var[N ( This latter situation corresponds to the mean field case in which C 2 = σ 2 /µ. So for areas of radius much smaller than the correlation length the model is characterized by α = 2 with logarithmic corrections, while in the case of radii much larger thanλ we obtain α = 1. This is in agreement with previous studies [26][27][28]. The model predicts α = 1 at large scales regardless of the dimension of the system, whereas for small areas α strongly depends on d, e.g. α = 2 (without logarithmic corrections) in the one dimensional case, and α = 5/3 for d = 3. Therefore, from small to large spatial scales the model predicts a cross-over of exponents which is difficult to explain without a spatially-explicit framework. Recently, considerable attention has been devoted to the origin of curvatures in scaling relationships, such as those relating body size and the metabolic rates of living organisms (e.g species of mammals [30] or freshwater phytoplankton [31]). Unlike what has been found in these latter works, in our model the cross-over in the scaling exponent α is a result of the interplay between space and dispersal at different spatial scales. When considering small areas, local communities appear strongly correlated with each other, while at larger areas these communities are completely independent as they are located at distances much larger than the correlation length. When they are decoupled, the relation is simply It has been proved that Taylor's law can emerge in a much more general class of stochastic processes, for example when the dynamical rates of the model are affected by environmental variability [28,29]. Here we do not consider this effect, but it would be certainly interesting to investigate environmental stochasticity within a spatially explicit framework.

VI. SOLUTION OF THE CONDITIONAL PDF
The goal of this section is to derive the main result of the paper. With the spatial population variance obtained in the previous section along with an appropriate approximation, we are now able to close and solve eq.(17) for the conditional generating function Z(h|V ). We will then invert this latter for the conditional probability density function P (N |V ), from which several characteristics of the spatial patterns can be deduced.
In order to calculate the explicit form of Z(h|V, T ) from eq.(17), we first need to calculate the solution of eq. (18), which in turn has to satisfy the identity in eq. (14). We will first make use of this latter in the form It is not difficult to verify that f can be expressed as where the functions A i are such that for any i. One could calculate the explicit form of A i , which depends on the spatial position, by substituting the expression in eq.(30) into eq. (18). However, the meaning of these functions provides a more efficient way for the calculation. Because of the definition of f , we readily obtain f (x, h = 0) = n and for m = 1, 2, . . ., from which we can make explicit the expression of A i by using eq. (30). For instance, it is not difficult to show (see Appendix E) that at stationarity where Ψ and ψ have been defined in eqs. (25) and (27), respectively. Similar relations, though more complicated, hold for i = 2, 3, . . .. Actually, n(x)N i can be calculated by integrating i times over the volume V the (i + 1)-th spatial correlation function. Nonetheless, note that this result depends on the symmetry of the d-dim volume V, and on its connectedness. This is important for exploiting the spherical symmetry of the system when introducing polar coordinates, and in selecting its center as the origin (see also the previous section and Appendix D) .
Since we are interested in relatively large population sizes when the correlation length is either large or small (λ → 0, ∞ but finite), we retain only the first two terms in the bracket of eq. (30). Thus, f at stationarity turns into where A 1 is the one we have obtained before. It is remarkable that, when substituting eq.(33) into eq. (17), at stationarity one obtains (see Appendix E) where Σ(R) = σ 2 ψ(R/λ)/µ is the spatial Fano factor defined in the previous section. The form of eq.(34) and that of eq. (19) at stationarity are the same, provided we replace σ 2 /µ with Σ(R). Therefore the solution of eq. (34) is which, when inverted for the probability density function, gives a gamma distribution of the form . (36) Thus, while without spatial effects the characteristic scale of the population size is −1 σ 2 /µ 1, in this regime space introduces a space-dependent scale for fluctuations which is quantified by Σ(R) = −1 ψ(R/λ), where ψ is the function defined in eq. (27).
As a further insight, if one defines the new process for the random variable N (R, T ) in the volume V of fixed radius R by the stochastic differential equatioṅ where ξ(t) is a zero-mean Gaussian white noise and ξ(t)ξ(t ) = 2δ(t−t ), then the stationary pdf of N (R) is exactly eq. (36). Notice that space is taken into account only implicitly through the functions V (R) and ψ(R/λ).
Eq.(37) can also be obtained as an -limit of a spatiallyimplicit master equation along the lines we have showed in Section IV. This process -unlike the spatial onesatisfies the detailed balance condition at stationarity as the flux at N = 0 is set to zero. This result suggests that there are some families of spatially-explicit processes which, when restricted to a finite volume, can be well approximated by spatially-implicit processes. Whilst the former brakes detailed balance, the latter turns out to be simpler and satisfies the detailed balance condition. In this model the region of this approximation is close to the critical point of the process.
It remains to understand when the f in the form of eq.(33) solves eq.(18), i.e. the original equation for f . In Appendix E we show that, for a fixed radius R, at stationarity f is a solution when h is small andλ is either very large or very small (but finite) compared to R, regardless of the spatial dimension d. The limitsλ → 0 and λ → +∞ of eq.(36) lead to the mean-field expressions, respectively eq. (19) and P (N |R) = δ(N − n V (R)) at leading order. Thus, eq.(36) captures the leading behavior of the distribution of the random variable N (R) in the large population regime and in the vicinity of the critical point. The simulations indeed confirm this with very good accuracy as shown in Figs. (2) and (3).
Because we assumed that A 1 (x) is at stationarity, eq.(33) does not give all the correct terms for a time evolution of the process N (R, T ). However, relatively close to stationarity, even the temporal dynamics is accurately described by eq.(37). We have checked this idea and compared the exact simulations of the process -as provided by the Doob-Gillespie algorithm (see also the following section) -with the corresponding analytic temporal evolution as obtained from solving eq. (37). Assuming that initially in the volume V there are N 0 individuals, we find the following time-dependent solution (see [32] for the details of the derivation): where we used reflecting boundary conditions at N = 0 at any t > 0. This pdf indeed tends to the stationary solution in eq.(36) as t → ∞. The agreement between simulations and eq.(38) is showed in Fig.(7) and (8).

Simulations with an efficient algorithm
These families of spatial stochastic models are difficult to simulate with parameters in arbitrary regimes and even so, usually only on relatively small lattices in low dimensions. Indeed, it is difficult to assess whether or not simulations have reached stationarity (because of the effect of large fluctuations) and how many replicates are necessary to get reliable predictions for the spatial mo-ments. In this light, it is even more important to know the analytical behavior of some quantities, which could not have been guessed from the simulations only. We have therefore compared our analytic predictions to simulations as obtained from a range of different parameter sets and from two different simulation schemes.
The first one is the Doob-Gillespie's algorithm [33] for producing exact trajectories of Markovian processes. We have used periodic boundary conditions in 1-d lattices of various sizes. Parameters were chosen so that the correlation length of the system was much smaller than the total size of the lattice. We have analyzed different sets of parameters and for each one we have run 50,000 independent realizations: error bars were calculated by grouping the results into 50 sets made of 1000 realizations each. In principle this simulation scheme allows us to obtain the exact trajectories of the system at any time, from an initial configuration up to stationarity. However, it is computationally very expensive, and therefore we have been forced to choose relatively small lattice sizes to investigate significant changes from the initial configuration. The results of this are shown in Figs. (7) and (8).
In order to analyze a wider set of parameters and larger lattices, we have simulated the process by using a new algorithm which generates the stationary random field obtained from the stochastic partial differential equation defined on the lattice. This was introduced in [34], and modifies a previous scheme that was used for simulating models of directed percolation [35]. Here we briefly summarize the main steps of the pseudo-code. From the definition of the discrete Laplace operator, we split the term λ 2 ∆ i n i in eq.(11) into one part depending only on n i (2 d λ 2 n i ) and another one depending on the densities in the nearest neighboring sites (λ 2 j:|i−j|=1 n j ). Conditional on the values of n j for j = i, the second term is constant and thus P (n i |n j ) can be obtained as a gamma distribution at stationarity [34]. Starting from a random initial configuration, at each iteration m, we randomly select a site i and update the value of n m i in the next step by sampling from n m+1 i ∼ P (n i |n m j ), conditional on the values of n m j (i.e. n m+1 j = n m j ). These steps are repeated until lattice configurations become independent of initial conditions. While standard stochastic integration schemes fail to preserve the positivity of {n} at any step, this algorithm produces non-negative populations by construction. We have verified that the simulated distribution thus obtained for N (R) matches the exact simulations of the Doob-Gillespie's algorithm in d = 1 and for different lengths, when the comparison was feasible (see Fig.(9)). By means of this algorithm we were able to study much larger lattice sizes at stationarity, and compare simulations against the predictions of the analytic solutions. The agreement was excellent in all the expected regimes (see Figs. (2) and (3)). the prediction based on the model. Histograms and black points represent the empirical data, the green solid lines are the predictions obtained from the solution of the model, i.e. from eq.(36), and the blue solid lines are those from the best-fit of the Fisher log-series. We highlight that the green lines are not best fits to the empirical data, but genuine predictions that are formulated from the empirical measures of the pair correlation function, as described in the main text (see eq. (7)). The radii of the areas are 15, 70 and 200 meters as reported on the corresponding panels. For the statistical analysis see Appendix G. In panel (d) we compare the empirical pair correlation function (black dots) with the best fit of eq.(7). The correlation length that is thus calculated isλ ≈ 2.5 × 10 3 meters, with ρ ≈ 8.9 × 10 3 and n ≈ 6.1 × 10 −4 trees per square meter for each species. The total number of species in the whole 50Ha forest stand is ≈ 900.

VII. AN ECOLOGICAL APPLICATION
A simple, but far from trivial, application of the mathematical model we have previously described is the modelling of spatial patterns in ecosystems with a large number of species. Examples include the species richness and abundance distribution of coral reefs, bees across landscapes or vascular plant species in tropical forest inventories. Starting with the crude approximation that species are independent (at least at some relatively large scale), this model can be used to predict the abundance distribution of species from measures of the two-point correlation function (PCF) and mean abundance per species. Because these latter descriptors are relatively easy to calculate, it turns out that we can obtain an estimate of how many rare (or abundant) species live in a region without surveying the entire system. Therefore, as well as being theoretically interesting, this approach has also an important practical advantage, because it ultimately allows one to infer the total number of species within a very large spatial region by utilizing only scattered and smallscale samples of the region itself. This is a long-lasting problem which has received a lot of attention recently [6,15,20].
Our goal here is not to explain the upscaling method, but only show to what degree the spatial model is in agreement with ecological empirical data. The two-point correlation function (PCF) specifies how similar individuals are distributed as a function of the geographic distance. Since species are assumed to be independent, the spatial distribution of individuals belonging to the same species can be considered as an independent realization of the stochastic process. As a consequence, the term n(x)n(y) / n 2 in eq.(7) can be easily calculated as the product of local abundances of individuals in sites that lay at the same distance, averaged across all species. Usually, it is more likely that close-by individuals belong to the same species than individuals that live farther apart. This translates into a PCF that is always positive, but decays with distance [36,37]. A decline in similarity with increasing geographic distance indicates that the individuals of a community are spatially aggregated. Therefore, a simple random placement of individuals in space is not a good approximation of the configurations of the community as thought in the past [38]. On the contrary, the stationarity PCF of this model decays with distance, is always positive (see eq. (7)) and has two free parameters (λ andρ). When we calculate these latter from a best fit of the empirical data, we obtain a good agreement as shown in Fig.(4d).
The mean abundance per species is readily available, because the total number of species and individuals are known in this forest plot. Since at stationarity the model is fully specified by these three parameters ( n ,λ and ρ), we are henceforth able to predict all the stationary patterns that we like to compare with those of the empirical ecosystem. One of them is the probability that a species has a given number of individuals within a specific region. In the ecological literature it is often referred to as species abundance distribution (SAD) [39]. In our model it is given by eq.(36) and the histograms of this pattern are reported in Fig.(4a-c) for different radii. The SAD represents one of the most commonly used static measures for summarizing information on ecosystem's diversity. Interestingly, the shape of the SAD in tropical forests has been observed to maintain similar features, regardless of the geographical location or the details of species interactions. Indeed, it often displays a unimodal shape at larger scales and a peak at small abundances at relatively small spatial scales. Numerous papers have focused on evaluating the processes that generate and maintain such observed characteristics [20,40], but only few of them considered spatial effects in an explicit framework [19,41].
The panels in Fig.4 show a comparison between the empirical data from the forest inventory of Pasoh Natural reserve in Malaysia (year 2005) and the predicted abundance distribution of species obtained from the model we have described in the previous sections. It is remarkable that the predicted curves in the first three panels are genuine inferences obtained from the mean abundance per species (i.e., n ) and the PCF (i.e.,λ andρ), and not best fitted curves to empirical SAD. As a comparison, we have included the best fit to data of a probability distribution commonly used as null model in the literature, the Fisher log-series (see Appendix G) [20]. The two methods have comparable accuracies at smaller scales, while at larger scales our method outperforms the Fisher logseries and captures the empirical SAD with much higher accuracy.
Such an agreement confirms that abundances of species are indeed characterized by very large fluctuations and that local populations appear correlated over very large spatial scales. This entails that complex ecosystems may comprise a large number of rare species, whereas only a few have large abundances (hyperdominant species). Because this model does not explicitly include interactions nor environmental forcing, we cannot single out which ecological processes bring about this finding. Nevertheless, the model is able to explain this separation of population size scales by poising a system close to criticality.
In terms of new physical insight, the theoretical framework that we have presented here allows us to connect exactly the two-point correlation function to the Species Abundance Distribution (SAD) and the Species Area Relationship (SAR) of a system. This explains quantitatively why and how SAR (known as α-diversity in the ecological literature) and species spatial turnover (known as β-diversity in the ecological literature) are related. For instance, it is usually assumed that the SAR is a power law function of the area [42], i.e. S(A) = cA z , however our results show that this can only be an approximation which works on a given range of spatial scales.
The simplicity and generality of the model make it suitable for the description of patterns in other biological systems. Indeed, numerous recent studies have showed that spatial bio-geographical patterns emerge in marine ecosystems [43,44], microorganisms [45], including bacteria [46,47], archaea, viruses, fungi [48] and eukaryotes [49,50]. Future work will include the application of the current framework to those biological communities.

VIII. CONCLUSIONS
In this paper we have studied a spatial stochastic model which can be fruitfully used to describe the main large scale characteristics of species-rich ecosystems. We have shown how to calculate analytically some of the most important spatial patterns when the system is close to criticality, which is the regime where the most important features emerge.
The model encapsulates birth, death, immigration and local hopping of individuals. It describes the dynamics of point-like and well-mixed individuals living in a metacommunity defined on a d-dimensional regular graph. The model is also minimal, meaning that, without one of its components (i.e., birth, death, nearest-neighbor hopping and external immigration), it yields either trivial or well-known results. Despite its simplicity, however, it violates detailed balance (see Appendix A) and generates a remarkable phenomenology of patterns, when the birth and death rates are comparable (criticality), leading to its properties being governed by large fluctuations, whose effect of course escape any classical mean-field analysis. These patterns and fluctuations entail strong correlations on large spatial and temporal scales.
The linearity of the rates is not sufficient to derive a full solution (in a weak sense) of the model. However, in applications one is usually interested in the analytical properties of processes that are much less general than the spatial random field. Thus we restricted our analysis to the conditional distribution, p(N |V ), that N individuals are found in a volume V . This quantity is sufficiently general to describe a wealth of patterns in several systems. We have found that in the close-to-critical regime, p(N |V ) satisfies an aptly derived equation which has the form of the corresponding mean-field equation of the process (i.e., without space). This equation includes functions of V , which we have exactly calculated. Such spatial redefinition of the parameters introduces strong deviations from the corresponding mean-field solutions, as confirmed by the exact stochastic simulations. Also, this shows that the process that governs the random variable N satisfies detailed balance in a first approximation and close to the critical point, thus being considerably simpler than the distribution of the random field. This result suggests that not only it is possible to make considerable analytical progress in this model, but there may be other, more general, models close to criticality which can be studied with a similar approach.
Indeed, our results suggest the tantalizing hypothesis that the conditional distribution, p(N |V ), provided by some families of spatial stochastic processes, is described by much simpler processes within a specific region of the parameter space. In our model the region of this approximation is close to the critical point of the process and the distribution of the simpler process holds the same functional shape across all spatial scales. Of course, the hypothesis requires much more scrutiny, especially when spatial models include nonlinear terms which can jeop-ardize the methods developed here. On the other side, our approach allows for the analysis of several generalizations, including non-local dispersal kernels and different sources of noise (e.g., environmental noise).
On a more applied side, our framework explains how the most important patterns in macro-ecology (e.g., SAR, SAD and PCF) are intrinsically connected with each other (see eqs. (7) and (36) and their relation through Σ(R) in eq.(28)), and also accounts for a crossover of power law behaviors in the spatial Taylor's law (see eq.(28)). We have compared the predictions against those of the Fisher log-series, commonly used as null model in the ecological literature. Despite the free parameter of the Fisher log-series distribution was calculated directly from the empirical data of the species abundance distribution (best fit) at each spatial scale, our model performed much better (see Fig.4), even though we employed the parameters obtained from the PCF curve, which does not provide information about the distribution of abundances of species. Discrepancies between empirical data and theoretical predictions will potentially inform us on the importance of alternative physical effects which will be included in more realistic models.
We have finally shown that, when applying the model to large ecosystems, species are predicted to display a broad range of abundances as a consequence of the critical regime. Therefore, many of them are rare and only a few are very abundant. These large demographic fluctuations are correlated across large spatial scales as well as over long times, in agreement with several empirical datasets collected in well-known tropical forest inventories [6,34,51].
Of course, these findings do not prove that these biological systems are close to criticality, but suggest that it is worth pursuing that route further. In this article we have not investigated the reasons why those systems look nearly critical, nor how they can operate within such tiny regions of their parameter space. For this, one needs to look into how inter-and intra-interactions dynamically lead living systems towards the correct region, in which states are biologically meaningful -however, see [52][53][54][55]. Nonetheless, our results could help understand what key factors drive such dynamics, and possibly shed light on the importance and effects of non-linearities among interacting agents.
We thank Amos Maritan for insightful discussions.

Appendix A BROKEN DETAILED BALANCE
We here recall the main general properties of detailed balance. Let us denote with c the configuration of a generic stochastic process and indicate with p(c, t|c 0 , t 0 ) the probability that the configuration c is seen at time t, given that the configuration at time t 0 was c 0 (abbreviated p(c, t)). We introduce W(c |c), which is the (time-independent) rate to transit from state c to c . If we consider Markovian dynamics, the evolution of p(c, t) is given by the following master equation (ME) [22,23] ∂p(c, t) ∂t = c W(c|c )p(c , t)−W(c |c)p(c, t) . (39) If it happens that W(c|c )P (c ) − W(c |c)P (c) = 0 for all configurations (this is the detailed balance (DB) condition), then the probability distribution P (c) is also a stationary solution of the ME, as we see from eq.(39).
On the other hand, it is clear that not all stationary distributions satisfy the detailed balance condition.
FIG. 5: Broken Detailed Balance. The picture represents a closed path in the space of configurations of the process which has different probabilities depending on the direction of the path (see also [17]). We take two neighboring sites, which initially have n and m individuals. The rates of jumping into the next configuration of the path are reported in the image and the arrows indicate the direction. In this simple case D = b(1 − γ). In the clockwise direction (red arrows) the total rate is [bγn + Dm + b 0 ][bγm + D(n + 1) + b 0 ][r(n + 1)][r(m + 1)].
In the anti-clockwise direction (blue arrows) the total rate is These two total rates must be equal for the detailed balance to hold. However, for any arbitrary configuration (m = n) this is true only when D(D − bγ)r = 0, i.e. for b = 0, r = 0, γ = 1, 1/2.
It is possible to show [12] that a condition for the validity of DB is that the probability of following a closed path in the space of configurations does not depend on the orientation of the path. More precisely, DB is satisfied if and only if for any choice of a closed path {c 1 , ..., c m }, with m an arbitrary number, the following holds W(c 1 |c 2 )W(c 2 |c 3 ) · · · W(c m |c 1 ) = (40) =W(c 1 |c m )W(c m |c m−1 ) · · · W(c 2 |c 1 ) .
This equation corresponds to microscopic reversibility and, when it is violated, the system can be found in nonequilibrium steady states.
Violation of DB is key to living systems and is the subject of intense recent interest [56,57]. In fig. (5) we show with a simple example that for this model the condition in eq.(40) is not satisfied, hence the detailed balance condition does not hold: Fig. (5) shows a path in the space of configurations for which the total rate does depend on the orientation of the closed path. Notice that such counterexample does not hold when spatial dispersal is switched off (i.e. D = 0 or γ = 1), or when autocatalytic production and spatial dispersal rates are equal (γ = 1/2). From eq.(2), multiplying both sides through by e s∈L nsHs and summing over all states, we find eq.(4) of the main text. If we differentiate by H k and impose {H} = 0, we find the equation for the mean number of individuals n k , reported in the main text. Taking another derivative with respect to H l and setting {H} = 0, we obtain the equation for the spatial two-point correlation function (PCF) among the sites k and l ∂ ∂t n k n l = D ∆ l n k n l + ∆ k n k n l − 2µ n k n l + +2b 0 n + δ k,l 2σ 2 n + b 0 + D∆ k n k , where σ 2 := b+r 2 , D := b(1−γ) 2d , δ k,l is a Kronecker delta, and ∆ is the discrete Laplace operator as defined in the main text. Considering stationary patterns, because of homogeneity we have ∆ k n k = 0 and, introducing G k,l = n k n l − n 2 , we obtain D(∆ k G k,l +∆ l G k,l ) − 2µG k,l + + 2σ 2 n + b 0 δ k,l = 0 .
We will now consider the Fourier series expansion of G k,l , which we will write as where k, l are the Cartesian coordinates of the locations of sites on the lattice (in a units), p is a vector with d components which belongs to C, the hypercubic ddimensional primitive unit cell of size 2π/a. Upon substituting G k,l in the stationary equation, we get an expression forĜ(p), which iŝ where p i is the i-th component of p. Therefore .
We can obtain a good deal of simplification by taking a continuous spatial limit (a → 0). Renaming p → p, k → x and l → y (now continuous variables in R d ), and appropriately rescaling the constants as explained in the main text, we arrive at the expression for the pairwise spatial correlation in dim d, i.e.
where we have replaced b 0 withb 0 in n and used the definitions ofρ,λ given in the main text. Actually, the ddim integral in the previous expression can be calculated dp e ip·(x−y) 1 +λ 2 p 2 = ∞ 0 ds dp e ip·(x−y)−(1+λ 2 p 2 )s where K a (z) is the modified Bessel function of the second kind and in the last step we used 9.6.24 from [58]. Eventually, the PCF is Because K a (z) decays exponentially for large z regardless of a,λ plays the role of a correlation length of the spatial system. Instead,ρ 2 has dimensions of a d-dim volume and provides a characteristic volume of local demographic fluctuations. Z and to the function g(i, k, h, t) = n i n k e hN (V ) . In this section we want to show that, at leading order in the limit a → 0, the value of ∆ i g(i, k, h, t) approaches ∆ k g(i, k, h, t). This makes the equation for f (k, h, t) much simpler, basically decoupling it with the other quantities. For simplicity we will only consider the one dimensional case, but this claim holds true in higher dimensions as well. Thus V will be an interval of length 2R and the origin of the coordinate system will be located at its center. By changing slightly the notation, we will now indicate N (V ) = N (−R, R), thus making explicit that V extends from site −R to R. Since the system is spatially homogeneous, we can write n i−a n k e hN (−R,R) = = n i n k+a e hN (−R+a,R+a) As a → 0 and at leading order, we can neglect a in the argument of N (−R, R), thus obtaining g(i − a, k, h, t) = n i−a n k e hN (−R,R) = = n i n k+a e hN (−R,R) = =g(i, k + a, h, t) Similarly n i+a n k e hN (−R,R) = n i n k−a e hN (−R,R) . Thus, as a → 0 at leading order we can write This approximation makes it possible to obtain eq.(18) in the continuous limit.

Appendix D MULTIDIMENSIONAL VARIANCE
In the main text we have outlined how to calculate the second moment of the random variable N (V ). Here we provide some more details. We will take V to be a d-dimensional ball with the origin of the Cartesian coordinates in its center. By integrating n(x)n(y) over y ∈ V in eq. (22) and using the symmetry of n(x)n(y) in x and y, we obtain an equation for n(x)N (V ) , which at stationarity reads where Θ(z) is the Heaviside step function and V is the with radius R. Eq.(42) must be solved for R < |x| and R > |x| separately. Boundary conditions for n(x)N (V ) and its first derivative will fix the values of the integrating constants. For |x| < R we obtain The constants A and B will be fixed using the aforementioned continuity conditions. Upon explicit calculation, for |x| ≤ R we obtain where Ψ(a, b) was defined in eq.(25) of the main text and I ν (z) and K ν (z) are the modified Bessel functions of the first and second kind, respectively [21]. Integrating with respect to x ∈ V and using the properties of I ν (z), we can readily obtain the explicit form of the second moment of N (V ), i.e.

Appendix E EVALUATION OF THE REGIMES OF ACCURACY OF THE METHOD
In this section we show that the truncation of f (x, V, h) as in eq.(33) yields an accurate approximation for the conditional probability distribution of the model. By retaining the first two terms in the square brackets of eq.(30), at stationarity we are left with the following where Z(h) is the conditional generating function at stationarity. By taking the derivative with respect to h of both sides of eq.(45) and setting h = 0, we obtain which gives Because in Appendix D we have already calculated n(x)N (V ) and N (V ) 2 (see eqs. (43) and (44)), A 1 (x, V ) is known explicitly. Substituting f (x, V, h) in eq.(45) with A 1 (x, V ) obtained before into eq.(17) at stationarity, we get Now we can readily simplify the termλ 2 ∇ 2 x n(x)N (V ) by making use of eq. (42). Integrating this latter with respect to x in V, since N (R) = n V and Σ(R) = ( N (R) 2 − N (R) 2 )/ N (R) (see Appendix D), we are therefore left with the following equation which therefore provides the equation for the generating function of N (R) up to terms O(h).
Similarly to what we have done so far, we can get further insight into the evolution of f (x, V, h). We substitute f (x, V, h) from eq.(45) into eq.(18) and use eq.(42) to obtainλ 2 ∇ 2 x n x N (V ) . Eventually, at stationarity this yields The first addend of eq.(47) is zero because of eq.(46), and the remaining terms are negligible when hA 1 is small.
Fixing the values of |x|, R and σ 2 /µ, we can calculate explicitly the regimes ofλ where A 1 (x, R) approaches zero. From eqs. (25) and (27) in the main text it is easy to rewrite A 1 (x, R) as The asymptotic expansion of the modified Bessel functions is [21] I ν (z) = e z √ 2πz when z → ∞. From this it is not difficult to verify that asλ → 0 we have The case ofλ → ∞ is more elaborate. Let's call z = R/λ. After some lengthy but otherwise straightforward calculations we can verify that as z → 0 we can write and at leading order ψ(z) is proportional to z 2 for d ≥ 3.
For the case of Ψ let us write z x = |x|/λ. We can verify that and at leading order Ψ is proportional to z 2 for d ≥ 3. Thus again A 1 (x, R) → 0. These findings follow from that f goes into mean-field regimes asλ → 0, ∞ (but finite) and hence all spatial terms go to zero.

Appendix F MODEL WITH INDEPENDENT DISPERSAL
In this section we show how the analysis that is undertaken in the main text can be extended also to a model where spatial dispersal is independent of birth.
As in the main text, the new model consists of a spatial meta-community where local communities are located on a d-dimensional regular graph (or lattice), and within each community individuals are treated as diluted, wellmixed and point-like particles.
The model is defined by the following birth-death dynamics: each individual dies at a constant death rate r and gives birth at a constant rate b, and communities are also colonized from the outside at a constant immigration rate b 0 However, individuals can now jump from one site to any of its 2d nearest neighboring sites at any time (not just after birth), and this happens with rate D.
Indicating with X i , i ∈ L, the individual living in site i, the reactions defining the model's dynamics are the following where j indicates a nearest neighbour of site i. Comparison with the reactions reported in section II shows that, indeed, now spatial movement is decoupled from birth events.
Let us now indicate with n i the number of individuals in site i and with P ({n}, t) the probability to find the system in the configuration {n} at time t. The master equation for P ({n}, t) then reads where the dots represent that all other occupation numbers remain as in {n} and it is intended that P (·) = 0 whenever any of the entrances is negative.
Similarly to the procedure followed in the main text, as a first step we introduce the spatial generating function of the model, defined in eq.(3). Multiplying through eq.(48) by e k∈L n k H k and averaging, we obtain the equation for ζ({H}, t), which reads The next steps follow exactly the same procedure undertaken in the main text. Thus, we define the parameter ε = 2(r−b) r+b and assume the following parameter scaling: b0 µ ε = O(1) as ε → 0 + . We introduce the parameter η = D σ 2 and fix the scaling η = O(ε) as ε → 0 + . We further assume that the generating function ζ({H}, t) is analytic at H i = 0 for any i and that the most important contribution to the equation of ζ({H}, t) comes from a negative real neighborhood of the origin with thickness O(ε). Taking the change of variables H i = εS i , we expand eq.(49) in powers of ε, assuming S i = O(1) and S i ≤ 0. With the definition of the following constants which are formally equivalent to those defined in section III, and retaining only the leading order in ε at eq.(49), we obtain with ∆ i the discrete Laplace operator. Dividing through by ε and rescaling time as T := µt we finally obtain which is exactly equal to eq.(9). Since all the results of this paper stem from this equation (or, equivalently, eq.(10)), the conclusions drawn in the main text for the model with 'seed dispersal' also hold true for the diffusion model presented in this section, provided the systems are close to the critical point (i.e., as ε → 0).

Appendix G A NULL MODEL: THE FISHER LOG-SERIES
The Log-series distribution was first proposed in 1943 by the statistician Ronald Fisher to describe the empirical abundance distribution of British moths and Malaysian butterflies [59]. From a set of assumptions involving the independence of species and the absence of spatial dispersal, Fisher derived the following formula, which provides the probability, P n , that a species has n individuals within an ecosystem: x n n (52) in which n > 0 and 0 < x < 1 is a free parameter that is usually determined via a best-fit to data. The Fisher Log-series is still largely used as a null-model in theoretical ecology [18,42], because of the low number of free parameters (only one) and of its straightforward mathematical derivation. Nonetheless, many limitations have emerged [19,20] in more recent years. Among these, the lack of a spatial structure has strongly limited its accuracy in describing spatial ecological patterns. In Fig.4 we have compared the predictions of our model to the best-fits of eq.(52) (rescaled by the total number of species) for three areas of different sizes. The histograms are log-scaled in the x-axis, meaning that the i-th bin counts the number of species that have abundances be-tween e i−1 and e i . In order to compare these empirical data to the analytical predictions, we first need to compute the probability that a species has abundances between e i−1 and e i , which is straightforward from eq.(36) and from eq. (52).
Our model at stationarity has a total of three free parameters (namely, n ,λ andρ). n has been obtained directly from the empirical data, whereasλ andρ have been calculated from the best fit of the theoretical PCF, i.e. eq.(7), to the empirical PCF, which is unrelated to the species abundance distribution. Once obtained the three parameters, we have predicted the species abundance distribution by using eq.(36), which is the curve showed in Fig.(4) along with the histograms of the empirical data.
Looking at the plots in Fig.(4), we see that our model prediction is very accurate in all three cases, whereas the Fisher log-series fails at fitting the distribution in panels (b) and (c). This is confirmed by standard chi-square analysis: 0.91 p-value in the first panel, and p-values of 0.7 in the other two cases; while the Fisher Log-series yields a p-value of 0.46 and 10 −6 in panels (a) and (b), and of ≈ 10 −21 in panel (c). More importantly, our model can predict the SAD at all spatial scales without any further assumptions on the system, whereas a method relying on the fit of eq.(52) introduces new parameters at each spatial scale (which in general are incompatible with each other) and is inapplicable at scales where we do not have empirical data. where each site initially contained exactly 100 individuals, with a total number of 200 lattice sites; parameters are the same as in Fig.(7). The three panels show comparisons of simulated data from the Doob-Gillespie algorithm and the analytic formula presented in eq. (38), and at different segment lengths (20,40 and 60 sites, respectively), and within each panel the three plots refer to different times T (T = 0.05 for red histograms, T = 0.5 for green histograms and T = 2 for blue histograms, units of µ −1 ). Histograms represent data from simulations, solid lines are the analytic predictions and error bars have length twice as much the standard deviation. For small T predictions and simulated data differ, as expected. However, as we approach stationarity the prediction improves significantly, and already at T = 0.5 the analytical and simulated distribution match very well.