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 formulae 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.

Popular Summary

A wide variety of biological systems—such as flocks of birds, an ensemble of neurons, and even entire ecosystems—can be described using the tools of statistical mechanics and stochastic processes. In recent years, researchers have found that living systems seem to occupy a special point in their parameter space known as a critical point, similar to critical points seen in phase transitions. Using the idea as a starting point, we have come up with a framework for predicting which spatial patterns can arise in a simple model of population dynamics.

We find that there is a class of spatial stochastic processes for which it is possible to calculate explicit expressions of important community patterns in the close-to-critical regime. Applying our theoretical framework to the spatial distribution of hundreds of tree species in a tropical rainforest in Malaysia, we also show that a model with simple birth, death, immigration, and spatial movement is able to predict the species’ abundances across multiple spatial scales when the dynamics is close to criticality.

This work constitutes a significant theoretical advancement in the analysis of spatial stochastic processes, provides tools for analyzing more realistic models of real-world communities of living organisms, and suggests underpinning mechanisms that emerge in macroscopic spatial patterns.

Figure 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 $\gamma$ or hop onto one of the $2d$ nearest neighbors with probability $1-\gamma$. 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).

Figure 2

Stationary distribution as obtained from the phenomenological algorithm ($d=1$). (a)–(c) present a comparison between the simulated model as obtained from the phenomenological algorithm (histograms) outlined in Sec. 6 and the stationary solution as calculated with Eq. (36) (blue solid line). The lattice comprises 500 sites in total, and we carry out 50 000 independent realizations in $d=1$ with periodic boundary conditions, where parameters are $D=30$, $b_0=0.5$, $\mu =0.01$, and $\sigma =10$, and, hence, $\overline{\lambda }\approx 55$. Panels show results for segments of different lengths which include 10, 20, and 60 adjacent sites, respectively. The size of error bars (black lines) is twice as much the standard deviation, while the red solid line represents the mean-field solution of the system [i.e., Eq. (36), where $\mathrm{\Sigma }={\sigma }^2/\mu$]. (d) presents the comparison between the simulated and the analytic pair correlation function at stationarity. Red dots are from simulations, while the blue solid line is the analytic solution given by Eq. (7).

Figure 3

Stationary distribution as obtained from the phenomenological algorithm ($d=2$). (a)–(c) present a comparison between the simulated model as obtained from the phenomenological algorithm (histograms) outlined in Sec. 6 and the stationary solution as calculated with Eq. (36) (blue solid line). The square lattice comprises $200\times 200$ sites. We carry out 50 000 independent realizations in $d=2$ with periodic boundary conditions, with parameters $D=20$, $b_0=0.1$, $\mu =0.1$, and $\sigma =10$, and, hence, $\overline{\lambda }\approx 15$. (a)–(c) show results for different areas with radii of length 5, 10, and 20, respectively. The size of error bars (black lines) is twice as much the standard deviation, while the red solid line represents the mean-field solution of the system [i.e., Eq. (36), where $\mathrm{\Sigma }={\sigma }^2/\mu$]. (d) presents the comparison between the simulated and the analytic pair correlation function at stationarity. Red dots are from simulations, while the blue solid line is the analytic solution given by Eq. (7).

Figure 4

Species abundance distribution from a lowland tropical forest and prediction from the model. (a)–(c) present the comparison between the empirical data of the distribution of species’ abundances from the lowland tropical forest inventory of the Pasoh natural reserve (Malaysia) and 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 m as reported on the corresponding panels. For the statistical analysis, see Appendix pp7. In (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 $\overline{\lambda }\approx 2.5\times {10}^3\mathrm{m}$, with $\overline{\rho }\approx 8.9\times {10}^3$ and $⟨n⟩\approx 6.1\times {10}^{-4}$ trees per square meter for each species. The total number of species in the whole 50 Ha forest stand is approximately 900.

Figure 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 Ref. [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-\gamma )$. In the clockwise direction (red arrows), the total rate is $[b\gamma n+Dm+b_0][b\gamma m+D(n+1)+b_0][r(n+1)][r(m+1)]$. In the anticlockwise direction (blue arrows), the total rate is $[b\gamma m+Dn+b_0][b\gamma n+D(m+1)+b_0][r(m+1)][r(n+1)]$. These two total rates must be equal for the detailed balance to hold. However, for any arbitrary configuration ($m\ne n$), this equality is true only when $D(D-b\gamma )r=0$, i.e., for $b=0$, $r=0$, and $\gamma =1,1/2$.

Figure 6

Behavior of the function $\psi (x)$. The figure shows the behavior of $\psi (x)$ as defined in Eq. (27) in dim $d=1$ (blue solid curve), $d=2$ (yellow solid curve), $d=3$ (green solid curve), and $d=4$ (red solid curve).

Figure 7

Time evolution with the Doob-Gillespie algorithm. (a)–(c) present a comparison between the simulated model as obtained from the Doob-Gillespie algorithm (histograms), the analytic prediction as defined in Eq. (38) (blue solid line), and the mean-field solution (red line). Simulations are carried out in $d=1$ with periodic boundary conditions, where each site initially contains exactly 100 individuals, with a total number of 200 lattice sites; $L$ comprises ten adjacent sites, and time $t$ is expressed in units of ${\mu }^{-1}$. The size of error bars are twice as much the standard deviation. Results in (d) are from the same set of simulated data (red dots, with error bars) and compare simulations to the analytic curve (blue line) of the pair correlation function at large times ($t=20$). Parameters in these simulations are $b=600$ and $r=601$, with $\gamma =0.5$ and $b_0=5$ ($\lambda \approx 12$).

Figure 8

Time comparison for different areas and different times in one dimension. Simulations are carried out in $d=1$ with periodic boundary conditions, where each site initially contains 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 ${\mu }^{-1}$). Histograms represent data from simulations, solid lines are the analytic predictions, and error bars have length twice as much as 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.

Figure 9

Comparison between Gillespie and phenomenological simulation scheme. The four panels show the conditional PDF $P(N|L)$ simulated with the same set of parameters and on lattices of the same size but using different algorithms at stationarity and in one dimension. Blue dots represent data from the phenomenological scheme (and blue lines are the respective error bars), while red dots are from the Doob-Gillespie algorithm (with respective error bars). Black solid lines are the analytic predictions from Eq. (36). The parameters are $b=600$, $d=601$, $\gamma =0.5$, and $b_0=5$ ($D=150$), and, hence, $\lambda \approx 12$. The lattices have 200 total sites with periodic boundary conditions.