Spatial patterns emerging from a stochastic process near criticality

TL;DR

This paper introduces a minimal stochastic model for community dynamics near criticality, analytically deriving spatial distribution properties that align with empirical biodiversity patterns in tropical rain forests.

Contribution

It provides an analytical framework for understanding spatial patterns in communities near critical points, incorporating fluctuations and violating detailed balance.

Findings

Analytical probability density functions for local populations derived.

Model predictions match observed biodiversity patterns.

Simulations confirm analytical results in relevant regimes.

Abstract

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.

Equations

$$X_i\stackrel{b\gamma }{\to }$$

$$X_i\stackrel{\frac{b(1-\gamma )}{2d}}{\to }$$

$$X_i\stackrel{r}{\to }$$

$$\varnothing \stackrel{b_0}{\to }$$

$${\mathcal{W}}_i^{+}(\{n\})=$$

$${\mathcal{W}}_i^{-}(\{n\})=$$

$$\frac{\partial }{\partial t}P(\{n\},t)=$$

$${\mathcal{W}}_i^{+}(\{...,n_i-1,...\})P(\{...,n_i-1,...\},t)+$$

$$-$$

$$+$$

$$-$$

$$\zeta (\{H\},t)=$$

$$=$$

$$\frac{\partial }{\partial t}\zeta (\{H\},t)=$$

$$\displaystyle+b\gamma\frac{\partial\zeta}{\partial H_{i}}+b_{0}\ \zeta\Big{]}+r(e^{-H_{i}}-1)\frac{\partial\zeta}{\partial H_{i}}\Big{\}}\;.$$

$$\frac{\partial ⟨n_k⟩}{\partial t}=\frac{b(1-\gamma )}{2d}{\Delta }_k⟨n_k⟩-\mu ⟨n_k⟩+b_0$$

$$\Delta_{k}f(k)=\sum_{j:\mid k-j\mid=1}\Big{(}f(j)-f(k)\Big{)}\quad.$$

$$\frac{\partial }{\partial t}⟨n_k{n}_l⟩=$$

$$-$$

$$+$$

$$D:=\frac{b(1-\gamma )}{2d}\text{and}{\sigma }^2:=\frac{r+b}{2}.$$

$$\lambda :=\sqrt{\frac{D}{\mu }}\text{and}\rho :=\sqrt{\frac{{\sigma }^2}{b_0}},$$

$$⟨n_{\text{k}}{n}_{\text{l}}⟩$$

$$\displaystyle\times\Big{(}\frac{a}{2\pi}\Big{)}^{d}\int_{\mathcal{C}}\,\mathrm{d}\textbf{p}\frac{e^{i\textbf{p}\cdot(\textbf{k}-\textbf{l})}}{1+2\lambda^{2}\sum_{i=1}^{d}(1-\cos(p_{i}a))}$$

$$\frac{⟨n(x)n(y)⟩}{⟨n{⟩}^2}=1+$$

$$\times$$

$$\bar{\lambda }:=\sqrt{\frac{\bar{D}}{\mu }}\bar{\rho }:=\sqrt{\frac{{\sigma }^2}{{\bar{b}}_0}},$$

$$\displaystyle\frac{\partial}{\partial t}\zeta(\{S\},t)=\sum_{i\in\mathbb{L}}\sigma^{2}S_{i}\Big{\{}$$

$$\displaystyle+\varepsilon S_{i}\frac{\partial\zeta}{\partial S_{i}}\Big{\}}$$

$$\displaystyle\frac{\partial}{\partial T}\zeta(\{S\},t)=\sum_{i\in\mathbb{L}}S_{i}\Big{\{}$$

$$\displaystyle+S_{i}\frac{\partial\zeta}{\partial S_{i}}\Big{\}}\quad,$$

$$\frac{\partial\zeta}{\partial T}=\sum_{i\in\mathbb{L}}H_{i}\Big{\{}\lambda^{2}\Delta_{i}\frac{\partial\zeta}{\partial H_{i}}-\frac{\partial\zeta}{\partial H_{i}}+\langle n\rangle\zeta+\frac{\sigma^{2}}{\mu}H_{i}\frac{\partial\zeta}{\partial H_{i}}\Big{\}}\quad,$$

$$\frac{\partial }{\partial T}P(\{n\},T)=$$

$$\displaystyle\sum_{i\in\mathbb{L}}\Big{\{}-\frac{\partial}{\partial n_{i}}$$

$$+$$

$$Z(h|V,t)=⟨e^{hN(V,t)}⟩={\int }_0^{\infty }\mathrm{d}N{e}^{hN}P(N|V,t)$$

$$H_i=\{\begin{array}{l}h\text{ if }i\in \mathcal{V}\\ 0\text{ if }i\notin \mathcal{V}\end{array}$$

$$\displaystyle\frac{\partial}{\partial T}Z(h|V,T)=h\Big{[}\lambda^{2}$$

$$+$$

$$\sum _{i\in \mathcal{V}}⟨n_i{e}^{hN}⟩=\frac{\partial Z}{\partial h}(h|V,T).$$

Full text

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 mora2011biological . 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 nykter2008gene , to cell growth furusawa2012adaptation , relatively small networks of neurons schneidman2006weak , flocks of birds cavagna2010scale and, possibly, tree populations in tropical forests tovo2017upscaling .

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 athreya2004branching , 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 cardy1996scaling , 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 jahnke2007solving ; mellor2016characterization . When stochastic processes violate detailed balance, they have a generator which is not self-adjoint pavliotis2014stochastic and different states are coupled by probability currents at microscopic level zia2007probability . 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 henkel2008non ; garcia2012noise ; shem2017solution . Nonetheless, in this paper we study a model which, despite violating microscopic detailed balance krapivsky2010kinetic ; grilli2012absence , 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 $d$-dimensional 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 formulæ in the expected regimes.

The rest of the paper is organized as follows: in section 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 formulæ; 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 $d$-dimensional lattice, $\mathbb{L}$, whose linear side is $a$. If $X_{i}$, $i\in\mathbb{L}$, indicates an individual living in site $i$, the reactions defining the model’s dynamics can be cast into the form

[TABLE]

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 $\gamma$, whereas it may hop onto one of the 2$d$ nearest neighbours with probability $1-\gamma$. 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 o2010field ; o2018cross .

Notice that (for $0\leq\gamma<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 non-linearities 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\simeq 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 azaele2016statistical .

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\in\mathbb{L}\}$, the linear birth and death rates in site $i$, i.e., $\mathcal{W}^{+}_{i}(\{n\})$ and $\mathcal{W}^{-}_{i}(\{n\})$, read respectively

[TABLE]

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

[TABLE]

where the dots represent that all other occupation numbers remain as in $\{n\}$ and it is intended that $P(\cdot)=0$ whenever any of the entrances is negative. The spatial generating function of the process is defined as

[TABLE]

where $H_{i}\leq 0$ for every $i\in\mathbb{L}$. Multiplying both sides of eq.(2) by $e^{\sum_{k^{\prime}\in\mathbb{L}}n_{k^{\prime}}H_{k^{\prime}}}$ and summing over all configurations of the system, we find that $\zeta(\{H\},t)$ satisfies the equation

[TABLE]

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)=\sum_{i\in\mathcal{V}}n_{i}(t)$, where $\mathcal{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$:

[TABLE]

where $\mu:=r-b$ and $\Delta_{k}$ is the discrete Laplace operator, which is defined as

[TABLE]

This finite difference equation can be solved in full generality and at stationarity, i.e. for $t\to\infty$, we get simply $\langle n\rangle=\frac{b_{0}}{\mu}$, regardless of any spatial location.

The pairwise spatial correlation between sites $k$ and $l$, i.e., $\langle n_{k}n_{l}\rangle$, 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):

[TABLE]

where we have used the notation

[TABLE]

We also introduce

[TABLE]

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 $\langle n_{\textbf{k}}n_{\textbf{l}}\rangle$ as a Fourier series expansion. Exploiting translation invariance the final expression of the solution reads (see Appendix B)

[TABLE]

where $p_{i}$ is the $i$-th Cartesian component of the $d$-dim vector p and $\mathcal{C}$ is the hypercubic primitive unit cell of size $2\pi/a$. Interestingly, pairs of sites de-correlate for $\gamma=1$ or $b=0$, when at stationarity we obtain $\langle n_{\textbf{k}}n_{\textbf{l}}\rangle-\langle n\rangle^{2}=c\delta_{\textbf{k},\textbf{l}}$, where $c$ is a constant depending only on the demographic parameters and $\delta_{\textbf{k},\textbf{l}}$ is a Kronecker delta. However, $c\neq\langle n\rangle$, pointing out that local fluctuations are non-Poissonian.

Eq.(6) is amenable to a continuous spatial limit, obtained as $a\to 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)

[TABLE]

where $|x-y|$ is the distance between the sites located at $x$ and $y$, $K_{\nu}$ is the modified Bessel function of the second kind of order $\nu$ lebedev2012special , and we have defined

[TABLE]

where $\bar{D}:=Da^{2}$ and $\bar{b}_{0}:=b_{0}/a^{d}$ for $a\to 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\to 0$ for $x\neq y$.

As the asymptotic behavior of $K_{\nu}$ as $z\to\infty$ is $K_{\nu}(z)\sim e^{-z}\sqrt{\frac{\pi}{2z}}$, $\bar{\lambda}$ is the correlation length of the system; $\bar{\rho}^{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|\gg a$ and $\bar{\lambda}\gg 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 gardiner2004handbook . 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 gardiner2004handbook ; van1992stochastic . The system-size expansion solves these difficulties, but it has to be applied when the size of the system becomes large van1992stochastic . 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, $\varepsilon$ and $\eta$, which identify a non-trivial region when $0<b<r$, but $b\to r$. This choice comes from the observation that communities of living organisms often appear to have very large demographic fluctuations. Several studies (see Volkov2005 ; Volkov2007 ; tovo2017upscaling ) 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 $\frac{2(r-b)}{r+b}$ of the order of $10^{-7}$ (see Fig. (4) and section VII for more details). We hence define $\varepsilon:=\frac{2(r-b)}{r+b}$ with the condition that $\frac{b_{0}}{\mu}\varepsilon=\mathcal{O}(1)$ as $\varepsilon\to 0^{+}$; in this way we obtain a constant $\rho^{2}=\mu/b_{0}\varepsilon$, which in real systems is large because usually $r/b_{0}\gg 1$. The parameter $\varepsilon$ indicates how close the system is to the critical point, regardless of spatial diffusion. With the independent parameter $\eta:=\frac{D}{\sigma^{2}}$, instead, we compare the importance of spatial diffusion with respect to demographic fluctuations. We will assume that $\eta=\mathcal{O}(\varepsilon)$ as $\varepsilon\to 0^{+}$ and hence $\eta/\varepsilon=\lambda^{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 $\eta/\varepsilon\to 0$, the model is equivalent to the mean field model without spatial diffusion, whereas for $\eta/\varepsilon\to\infty$ 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 $\varepsilon$ and $\eta$, we assumed that the generating function $\zeta(\{H\},t)$ is analytic at $H_{i}=0$ for any $i$ and the most important contribution to the equation of $\zeta(\{H\},t)$ comes from a negative neighborhood of the origin with thickness $\mathcal{O}(\varepsilon)$. This is tantamount to introducing the change of variable $H_{i}=\varepsilon S_{i}$ into eq.(4) and to expanding in powers of $\varepsilon$, assuming $S_{i}=\mathcal{O}(1)$ and $S_{i}\leq 0$. Up to linear order in $\eta$ and $\varepsilon$ we obtain

[TABLE]

or

[TABLE]

where we have introduced the dimensionless time $T:=\mu t$ and now $\zeta$, 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

[TABLE]

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 $\sigma^{2}/\mu$. Eq.(10) leads to the following Fokker-Planck equation:

[TABLE]

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 $\eta$ and $\varepsilon$.

IV.1 Conditional generating function

Eq.(IV) 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 $\mathcal{V}$ the set of sites belonging to a $d$-dim volume $|\mathcal{V}|=V$, and let us introduce the random variable $N(V,t)=\sum_{i\in\mathcal{V}}n_{i}(t)$, i.e. the total number of individuals in $\mathcal{V}$ at time $t$. By indicating with $P(N|V,t)$ the corresponding probability density function of $N$, we define the conditional generating function

[TABLE]

where $h\leq 0$. We obtain the corresponding equation for $Z$ by specifying $H_{i}$, i.e.

[TABLE]

and substituting this into eq.(10). Thus,

[TABLE]

where $f(i,h,V,T)=\langle n_{i}(T)e^{hN(V,T)}\rangle$ and we have used the identity

[TABLE]

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\in\mathcal{V}$, we finally obtain the following equation for $f(k,h,V,T)$ (abbreviated $f(k,h)$)

[TABLE]

where $g(i,k,h,V,T):=\langle n_{i}(T)n_{k}(T)e^{hN(V,T)}\rangle$ and where we have used the identity

[TABLE]

In eq.(15) the function $g$ is still unknown, but it is possible to show that $\Delta_{i}g(i,k,h,T)=\Delta_{k}g(i,k,h,T)$ to leading order as $a\to 0$ (see Appendix C). This allows us to obtain a closed equation for $f$. Indeed, we have

[TABLE]

after using the identity (16).

The spatial continuous limit of eq.(13) as $a\to 0$ reads

[TABLE]

where in $\langle n\rangle$ we have substituted $b_{0}$ with $\bar{b}_{0}$. This equation is of crucial importance in what follows. Similarly, eq.(15) becomes

[TABLE]

where with the continuous coordinates ($n_{i}\to n(x)$ and $n_{k}\to n(y)$), we have also $f(x,h,V,T)=\langle n(x)e^{hN}\rangle$ and $N(V,t)=\int_{\mathcal{V}}\,\mathrm{d}x\ n(x,t)$. When there are no spatial effects, i.e., $\bar{\lambda}=0$, eq.(17) reads

[TABLE]

which has the following solution at stationarity

[TABLE]

where $h\leq 0$ and $\mu\langle n\rangle V/\sigma^{2}=V/\bar{\rho}^{2}$. This is the generating function of a gamma distribution van1992stochastic with mean $\langle n\rangle V$ and variance $\langle n\rangle^{2}V\bar{\rho}^{2}$. Also, it is not difficult to verify that when $\bar{\lambda}=0$, the solution of eq.(18) is given by

[TABLE]

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)=\int_{\mathcal{V}}\,\mathrm{d}x\,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\to 0$, which at stationarity reads

[TABLE]

where $\delta(x-y)$ is a Dirac delta and we have included only the leading terms as $\varepsilon\to 0$. By integrating both sides of eq.(22) with respect to $y\in\mathcal{V}$, we obtain an equation for $\langle n(x)N(V)\rangle$. 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 $\mathcal{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 $\langle n(x)n(y)\rangle$ with respect to $x$ and $y$, we obtain at stationarity

[TABLE]

This linear ODE has to be solved separately for $|x|>R$ and $|x|<R$. The continuity of $\langle n(x)N(R)\rangle$ and its first derivative at the boundary $|x|=R$ provide the solution. The final results is (see Appendix D)

[TABLE]

where the function $\Psi$ takes the following form for $|x|\leq R$

[TABLE]

where $I_{\nu}(z)$ and $K_{\nu}(z)$ are the modified Bessel functions of the first and second kind, respectively lebedev2012special . Integrating both sides of eq.(24) with respect to $x\in\mathcal{V}$, we obtain the equation for the second moment of $N(R)$, i.e.

[TABLE]

where $\psi(R/\bar{\lambda})$ takes the following explicit form in dimension $d$ (see Appendix D for its behavior):

[TABLE]

These two functions, namely $\langle n(x)N(R)\rangle$ in eq.(24) and $\langle N(R)\rangle^{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.

[TABLE]

where $\langle N(R)\rangle=\langle n\rangle V=\bar{b}_{0}V/\mu$ and $\Sigma(R):=\sigma^{2}\psi(R/\bar{\lambda})/\mu$. The function $\Sigma(R)$ is the spatial Fano factor and quantifies the deviations of the fluctuations from a Poisson process. Since $\sigma^{2}/\mu=\mathcal{O}((r-b)^{-1})$, when the system is close to the critical point – for fixed $R$ and $\bar{\lambda}$ – the system has large fluctuations on all scales larger than the correlation length $\bar{\lambda}$. Also, in the regime $R/\bar{\lambda}\to+\infty$ we obtain $\mu\Sigma(R)/\sigma^{2}=1+\mathcal{O}[(R/\bar{\lambda})^{-1}]$, thus recovering the mean-field fluctuations as predicted by eq.(20).

V.1 The spatial Taylor’s law

Taylor’s law was first observed in ecological communities taylor1961aggregation ; taylor1977aggregation , 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 eisler2008fluctuation ; james2018zipf . The law postulates a relation of the following form

[TABLE]

where $C$ is a positive constant and $\alpha$ typically assumes values between one and two taylor1977aggregation ; james2018zipf . 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 $\bar{\lambda}$, we obtain $\text{Var}[N(R)]=C_{1}\langle N(R)\rangle^{2}\log(\langle N(R)\rangle)$ for $R\ll\bar{\lambda}$, while $\text{Var}[N(R)]=C_{2}\langle N(R)\rangle$ for $R\gg\bar{\lambda}$. This latter situation corresponds to the mean field case in which $C_{2}=\sigma^{2}/\mu$. So for areas of radius much smaller than the correlation length the model is characterized by $\alpha=2$ with logarithmic corrections, while in the case of radii much larger than $\bar{\lambda}$ we obtain $\alpha=1$. This is in agreement with previous studies taylor1961aggregation ; taylor1977aggregation ; eisler2008fluctuation . The model predicts $\alpha=1$ at large scales regardless of the dimension of the system, whereas for small areas $\alpha$ strongly depends on $d$, e.g. $\alpha=2$ (without logarithmic corrections) in the one dimensional case, and $\alpha=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 r01748 or freshwater phytoplankton Zaoli17323 ). Unlike what has been found in these latter works, in our model the cross-over in the scaling exponent $\alpha$ 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 $\text{Var}[N(R)]=b\sigma^{2}V(R)/\mu^{2}$, where $V(R)$ is the volume of a $d$-dimensional sphere.

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 eisler2008fluctuation ; james2018zipf . 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

[TABLE]

It is not difficult to verify that $f$ can be expressed as

[TABLE]

where the functions $A_{i}$ are such that

[TABLE]

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)=\langle n\rangle$ and

[TABLE]

for $m=1,2,\ldots$, 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

[TABLE]

where $\Psi$ and $\psi$ have been defined in eqs.(25) and (27), respectively. Similar relations, though more complicated, hold for $i=2,3,\ldots$. Actually, $\langle n(x)N^{i}\rangle$ can be calculated by integrating $i$ times over the volume $\mathcal{V}$ the $(i+1)$-th spatial correlation function. Nonetheless, note that this result depends on the symmetry of the $d$-dim volume $\mathcal{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 ($\bar{\lambda}\to 0,\infty$ but finite), we retain only the first two terms in the bracket of eq.(30). Thus, $f$ at stationarity turns into

[TABLE]

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)

[TABLE]

where $\Sigma(R)=\sigma^{2}\psi(R/\bar{\lambda})/\mu$ 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 $\sigma^{2}/\mu$ with $\Sigma(R)$. Therefore the solution of eq.(34) is

[TABLE]

which, when inverted for the probability density function, gives a gamma distribution of the form

[TABLE]

Thus, while without spatial effects the characteristic scale of the population size is $\epsilon^{-1}\simeq\sigma^{2}/\mu\gg 1$, in this regime space introduces a space-dependent scale for fluctuations which is quantified by $\Sigma(R)=\epsilon^{-1}\psi(R/\bar{\lambda})$, where $\psi$ 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 equation

[TABLE]

where $\xi(t)$ is a zero-mean Gaussian white noise and $\langle\xi(t)\xi(t^{\prime})\rangle=2\delta(t-t^{\prime})$, 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 $\psi(R/\bar{\lambda})$. Eq.(37) can also be obtained as an $\epsilon$-limit of a spatially-implicit master equation along the lines we have showed in Section IV. This process – unlike the spatial one – satisfies 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 $\bar{\lambda}$ is either very large or very small (but finite) compared to $R$, regardless of the spatial dimension $d$. The limits $\bar{\lambda}\to 0$ and $\bar{\lambda}\to+\infty$ of eq.(36) lead to the mean-field expressions, respectively eq.(19) and $P(N|R)=\delta(N-\langle n\rangle 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 azaele2006dynamical for the details of the derivation):

[TABLE]

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\to\infty$. 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 moments. 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 Gillespie1977exact 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 peruzzo2017phenomenological , and modifies a previous scheme that was used for simulating models of directed percolation Dornic2005integration . Here we briefly summarize the main steps of the pseudo-code. From the definition of the discrete Laplace operator, we split the term $\lambda^{2}\Delta_{i}n_{i}$ in eq.(IV) into one part depending only on $n_{i}$ ($2^{d}\ \lambda^{2}n_{i}$) and another one depending on the densities in the nearest neighboring sites ($\lambda^{2}\sum_{j:|i-j|=1}n_{j}$). Conditional on the values of $n_{j}$ for $j\neq i$, the second term is constant and thus $P(n_{i}|n_{j})$ can be obtained as a gamma distribution at stationarity peruzzo2017phenomenological . 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}\sim 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)).

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 small-scale samples of the region itself. This is a long-lasting problem which has received a lot of attention recently tovo2017upscaling ; shem2017solution ; azaele2016statistical .

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 $\langle n(x)n(y)\rangle/\langle n\rangle^{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 condit2002beta ; chave2002spatially . 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 coleman1981random . On the contrary, the stationarity PCF of this model decays with distance, is always positive (see eq.(7)) and has two free parameters ($\bar{\lambda}$ and $\bar{\rho}$). 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 ($\langle n\rangle$, $\bar{\lambda}$ and $\bar{\rho}$), 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) mcgill2007species . 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 black2012stochastic ; azaele2016statistical , but only few of them considered spatial effects in an explicit framework pigolotti2018stochastic ; o2018cross .

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., $\langle n\rangle$) and the PCF (i.e., $\bar{\lambda}$ and $\bar{\rho}$), 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) azaele2016statistical . The two methods have comparable accuracies at smaller scales, while at larger scales our method outperforms the Fisher log-series 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 $\alpha$-diversity in the ecological literature) and species spatial turnover (known as $\beta$-diversity in the ecological literature) are related. For instance, it is usually assumed that the SAR is a power law function of the area hubbell2001unified , 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 rinaldo2002cross ; ser2018ubiquitous , microorganisms woodcock2007neutral , including bacteria giometto2014emerging ; giometto2015sample , archæa, viruses, fungi hanson2012beyond and eukaryotes altermatt2015big ; lynch2015ecology . 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 jeopardize 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 $\Sigma(R)$ in eq.(28)), and also accounts for a cross-over 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 peruzzo2017phenomenological ; tovo2017upscaling ; azaele2015towards .

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 mastromatteo2011criticality ; goudarzi2012emergent ; stieg2012emergent ; hidalgo2014information . 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.

IX Acknowledgments

FP thanks the NERC SPHERES DTP (NE/L002574/1) for funding his studentship. 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(\textbf{c},t|\textbf{c}_{0},t_{0})$ the probability that the configuration c is seen at time $t$, given that the configuration at time $t_{0}$ was $\textbf{c}_{0}$ (abbreviated $p(\textbf{c},t)$). We introduce $\mathcal{W}(\textbf{c}^{\prime}|\textbf{c})$, which is the (time-independent) rate to transit from state c to $\textbf{c}^{\prime}$. If we consider Markovian dynamics, the evolution of $p(\textbf{c},t)$ is given by the following master equation (ME) van1992stochastic ; gardiner2004handbook

[TABLE]

If it happens that $\mathcal{W}(\textbf{c}|\textbf{c}^{\prime})P(\textbf{c}^{\prime})-\mathcal{W}(\textbf{c}^{\prime}|\textbf{c})P(\textbf{c})=0$ for all configurations (this is the detailed balance (DB) condition), then the probability distribution $P(\textbf{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.

It is possible to show zia2007probability 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 $\{\textbf{c}_{1},...,\textbf{c}_{m}\}$, with $m$ an arbitrary number, the following holds

[TABLE]

This equation corresponds to microscopic reversibility and, when it is violated, the system can be found in non-equilibrium steady states.

Violation of DB is key to living systems and is the subject of intense recent interest gnesotto2018broken ; battle2016broken . 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 $\gamma=1$), or when autocatalytic production and spatial dispersal rates are equal ($\gamma=1/2$).

Appendix B Spatial correlation of the birth-death Markov process

In the main text we have defined the spatial generating function of the model

[TABLE]

From eq.(2), multiplying both sides through by $e^{\sum_{s\in\mathbb{L}}n_{s}H_{s}}$ 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 $\langle n_{k}\rangle$, 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$

[TABLE]

where $\sigma^{2}:=\frac{b+r}{2}$, $D:=\frac{b(1-\gamma)}{2d}$, $\delta_{k,l}$ is a Kronecker delta, and $\Delta$ is the discrete Laplace operator as defined in the main text.

Considering stationary patterns, because of homogeneity we have $\Delta_{k}\langle n_{k}\rangle=0$ and, introducing $G_{k,l}=\langle n_{k}n_{l}\rangle-\langle n\rangle^{2}$, we obtain

[TABLE]

We will now consider the Fourier series expansion of $G_{\textbf{k},\textbf{l}}$, which we will write as

[TABLE]

where $\textbf{k},\textbf{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 $\mathcal{C}$, the hypercubic $d$-dimensional primitive unit cell of size $2\pi/a$. Upon substituting $G_{\textbf{k},\textbf{l}}$ in the stationary equation, we get an expression for $\hat{G}(\textbf{p})$, which is

[TABLE]

where $p_{i}$ is the $i$-th component of p. Therefore

[TABLE]

We can obtain a good deal of simplification by taking a continuous spatial limit ($a\to 0$). Renaming $\textbf{p}\to p$, $\textbf{k}\to x$ and $\textbf{l}\to y$ (now continuous variables in $\mathbb{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.

[TABLE]

where we have replaced $b_{0}$ with $\bar{b}_{0}$ in $\langle n\rangle$ and used the definitions of $\bar{\rho},\bar{\lambda}$ given in the main text. Actually, the $d$-dim integral in the previous expression can be calculated

[TABLE]

where $K_{a}(z)$ is the modified Bessel function of the second kind and in the last step we used 9.6.24 from abramowitz1965handbook . Eventually, the PCF is

[TABLE]

Because $K_{a}(z)$ decays exponentially for large $z$ regardless of $a$, $\bar{\lambda}$ plays the role of a correlation length of the spatial system. Instead, $\bar{\rho}^{2}$ has dimensions of a $d$-dim volume and provides a characteristic volume of local demographic fluctuations.

Appendix C Approximating $g(i,j,h,t)$ in the continuous spatial limit

Equation (15) shows that the evolution of $f(k,h,t)$ is not closed, being coupled to the generating function $Z$ and to the function $g(i,k,h,t)=\langle n_{i}n_{k}e^{hN(V)}\rangle$. In this section we want to show that, at leading order in the limit $a\to 0$, the value of $\Delta_{i}g(i,k,h,t)$ approaches $\Delta_{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 $\mathcal{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

[TABLE]

As $a\to 0$ and at leading order, we can neglect $a$ in the argument of $N(-R,R)$, thus obtaining

[TABLE]

Similarly $\langle n_{i+a}n_{k}e^{hN(-R,R)}\rangle=\langle n_{i}n_{k-a}e^{hN(-R,R)}\rangle$. Thus, as $a\to 0$ at leading order we can write

[TABLE]

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 $\mathcal{V}$ to be a $d$-dimensional ball with the origin of the Cartesian coordinates in its center. By integrating $\langle n(x)n(y)\rangle$ over $y\in\mathcal{V}$ in eq.(22) and using the symmetry of $\langle n(x)n(y)\rangle$ in $x$ and $y$, we obtain an equation for $\langle n(x)N(V)\rangle$, which at stationarity reads

[TABLE]

where $\Theta(z)$ is the Heaviside step function and $V$ is the $d$-dim volume

[TABLE]

with radius $R$. Eq.(D) must be solved for $R<|x|$ and $R>|x|$ separately. Boundary conditions for $\langle n(x)N(V)\rangle$ and its first derivative will fix the values of the integrating constants. For $|x|<R$ we obtain

[TABLE]

while for $|x|>R$ instead

[TABLE]

The constants $A$ and $B$ will be fixed using the aforementioned continuity conditions. Upon explicit calculation, for $|x|\leq R$ we obtain

[TABLE]

where $\Psi(a,b)$ was defined in eq.(25) of the main text and $I_{\nu}(z)$ and $K_{\nu}(z)$ are the modified Bessel functions of the first and second kind, respectively lebedev2012special . Integrating with respect to $x\in\mathcal{V}$ and using the properties of $I_{\nu}(z)$, we can readily obtain the explicit form of the second moment of $N(V)$, i.e.

[TABLE]

where we can read off the explicit expression for $\psi(R/\lambda)$, also reported in eq.(27). Finally, since $\text{Var}(N(V))=\langle N(V)^{2}\rangle-\langle N(V)\rangle^{2}$ we can explicitly write down the variance of $N(V)$, which takes the form

[TABLE]

where $\Sigma(R/\bar{\lambda})=\sigma^{2}\psi(R/\bar{\lambda})/\mu$. $\Sigma(R/\bar{\lambda})$ is the spatial Fano factor and quantifies the relative importance of fluctuations in the system. Finally, Fig.(6) shows the behavior of $\psi(x)$.

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

[TABLE]

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

[TABLE]

which gives

[TABLE]

Because in Appendix D we have already calculated $\langle n(x)N(V)\rangle$ and $\langle N(V)^{2}\rangle$ (see eqs.(D) 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

[TABLE]

Now we can readily simplify the term $\bar{\lambda}^{2}\nabla^{2}_{x}\langle n(x)N(V)\rangle$ by making use of eq.(D). Integrating this latter with respect to $x$ in $\mathcal{V}$, since $\langle N(R)\rangle=\langle n\rangle V$ and $\Sigma(R)=(\langle N(R)^{2}\rangle-\langle N(R)\rangle^{2})/\langle N(R)\rangle$ (see Appendix D), we are therefore left with the following equation

[TABLE]

which therefore provides the equation for the generating function of $N(R)$ up to terms $\mathcal{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.(D) to obtain $\bar{\lambda}^{2}\nabla^{2}_{x}\langle n_{x}N(V)\rangle$. Eventually, at stationarity this yields

[TABLE]

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 $\sigma^{2}/\mu$, we can calculate explicitly the regimes of $\bar{\lambda}$ 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

[TABLE]

The asymptotic expansion of the modified Bessel functions is lebedev2012special

[TABLE]

when $z\to\infty$. From this it is not difficult to verify that as $\bar{\lambda}\to 0$ we have

[TABLE]

and so for $|x|<R$ indeed $A_{1}(x,R)\to 0$.

The case of $\bar{\lambda}\to\infty$ is more elaborate. Let’s call $z=R/\bar{\lambda}$. After some lengthy but otherwise straightforward calculations we can verify that as $z\to 0$ we can write

[TABLE]

and at leading order $\psi(z)$ is proportional to $z^{2}$ for $d\geq 3$.

For the case of $\Psi$ let us write $z_{x}=|x|/\bar{\lambda}$. We can verify that

[TABLE]

and at leading order $\Psi$ is proportional to $z^{2}$ for $d\geq 3$. Thus again $A_{1}(x,R)\to 0$. These findings follow from that $f$ goes into mean-field regimes as $\bar{\lambda}\to 0,\infty$ (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, well-mixed 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 2$d$ nearest neighboring sites at any time (not just after birth), and this happens with rate $D$.

Indicating with $X_{i}$, $i\in\mathbb{L}$, the individual living in site $i$, the reactions defining the model’s dynamics are the following

[TABLE]

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

[TABLE]

where the dots represent that all other occupation numbers remain as in $\{n\}$ and it is intended that $P(\cdot)=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^{\sum_{k\in\mathbb{L}}n_{k}H_{k}}$ and averaging, we obtain the equation for $\zeta(\{H\},t)$, which reads

[TABLE]

The next steps follow exactly the same procedure undertaken in the main text.

Thus, we define the parameter $\varepsilon=\frac{2(r-b)}{r+b}$ and assume the following parameter scaling: $\frac{b_{0}}{\mu}\varepsilon=\mathcal{O}(1)$ as $\varepsilon\to 0^{+}$. We introduce the parameter $\eta=\frac{D}{\sigma^{2}}$ and fix the scaling $\eta=\mathcal{O}(\varepsilon)$ as $\varepsilon\to 0^{+}$. We further assume that the generating function $\zeta(\{H\},t)$ is analytic at $H_{i}=0$ for any $i$ and that the most important contribution to the equation of $\zeta(\{H\},t)$ comes from a negative real neighborhood of the origin with thickness $\mathcal{O}(\varepsilon)$. Taking the change of variables $H_{i}=\varepsilon S_{i}$, we expand eq.(49) in powers of $\varepsilon$, assuming $S_{i}=\mathcal{O}(1)$ and $S_{i}\leq 0$. With the definition of the following constants

[TABLE]

which are formally equivalent to those defined in section III, and retaining only the leading order in $\varepsilon$ at eq.(49), we obtain

[TABLE]

with $\Delta_{i}$ the discrete Laplace operator. Dividing through by $\varepsilon$ and rescaling time as $T:=\mu t$ we finally obtain

[TABLE]

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 $\varepsilon\to 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 fisher1943relation . 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:

[TABLE]

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 hubbell2001unified ; o2010field , because of the low number of free parameters (only one) and of its straightforward mathematical derivation. Nonetheless, many limitations have emerged azaele2016statistical ; o2018cross 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 between $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, $\langle n\rangle$, $\bar{\lambda}$ and $\bar{\rho}$). $\langle n\rangle$ has been obtained directly from the empirical data, whereas $\bar{\lambda}$ and $\bar{\rho}$ 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 $\approx~{}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.