## The deterministic model

We start from the classic Bever deterministic model (Bever 2003) that couples the population dynamics of two competing plant species, A and B, to their soil, *S**A* and *S**B* = 1 - *S**A*. The system of coupled ordinary differential equations (ODE) that describes this model is (Eqs. 1):

$$\left\{\begin{array}{c}\frac{dA}{dt}=A ({r}_{A}-{m}_{A}A-{c}_{BA }B+{\alpha }_{A}{S}_{a}+{\beta }_{A}(1-{S}_{A})\\ \frac{dB}{dt}=B ({r}_{B}-{m}_{B}B-{c}_{AB }A+{\alpha }_{B}{S}_{a}+{\beta }_{B}(1-{S}_{A})\\ \frac{d{S}_{A}}{dt}={S}_{A}\left(1-{S}_{A}\right)\left[\frac{A}{A+B}-v\frac{B}{A+B}\right] \end{array}\right.$$

Where A and B are the population size of species A and B, and *S**A* is the soil axis variable that reaches is maximum value when the soil communities equals the one associated to a monoculture of plant *A*. The parameter *r* and *m* are the intrinsic growth rate and logistic term of the classic logistic population model, and the terms “*c*” are the competition coefficients. Note that while *r*, *m* and *c* are strictly positive, the parameters *α**a*, *α**b*, *β**a* and *β**b* are either positive or negative. The two alphas parameters (*α**a* and *α**b*) describe the net effect of soil of A on plant A and B, while the beta parameters (*β**a* and *β**b*) describe the net effect of the soil of B on plant A and B. The soil state variable thus can fluctuate between the soil typical of A and that typical of B. Finally, the parameter *v* describes the impact of B on the soil of A. The behaviour of this model is well known (Bever, 2003), and in the supplementary materials (RScript_Bever_Model_Stochastic.R) we provide an R script that reproduces the behaviour of the model with an example based on the same parameters originally presented by Bever (2003). If the system starts away from the equilibrium values, and the combination of parameters is such that the system allows plant coexistence, the trajectories of the system will be characterised by sinusoidal oscillations that quickly dampen until the system reaches the equilibrium, which is stable (Fig. 1a).

## Stochastic model: formulation

We introduce environmental stochasticity to the ODE system of Eqs. 1 through the Itô calculus formalism for continuous stochastic differential equations (SDE). SDEs for interacting populations are constructed by combining the deterministic component (the ODE of Eqs. 1 for us), which represents the “drift” term, and either a continuous time Markov chain describing probabilistic transitions within the population due demographic stochasticity (Allen 2010) or an external source of random variability due to environmental fluctuations (Allen 2010; Dobrow 2016), or a combination of both demographic and environmental stochasticity (Lande et al. 2003). Here, we are mostly interested in broad scale population dynamics, where the impact of demographic stochasticity is likely to be minor given the size of the population under consideration and considering that, typically, demographic stochasticity scales as the inverse of population size (Lande et al. 2003). We thus consider a general model with environmental stochasticity only, which can in the future be expanded to include also demographic stochasticity. The general SDE model then is:

$${d}{N}={\mu }\left({N}\right) dt+{\Sigma } {d}{W}$$

Where the bold notation indicates vectors, and so N is a vector with population sizes A, B and the soil state variable, dW is a vector of independent Wiener processes (often referred to as white noise in ecology), which we use as the source of environmental stochasticity (Lande et al. 2003; Dobrow 2016), and the matrix Σ is the diffusion matrix, which quantifies the intensity of stochasticity for each state variable. The drift term **µ**, which is a function of population size, corresponds to the deterministic structure of our ODEs, that is the classic Bever model (Bever 2003). If plant A, B and Soil were affected by independent environmental stochasticity, the diffusion matrix would simply be

$${\Sigma }=\left[\begin{array}{ccc}{\sigma }_{A}& 0& 0\\ 0& {\sigma }_{B}& 0\\ 0& 0& {\sigma }_{s}\end{array}\right]$$

And the corresponding SDEs then is

$$\left\{\begin{array}{c}dA=[A \left({r}_{A}-{m}_{A}A-{c}_{BA }B+{\alpha }_{A}{S}_{a}+{\beta }_{A}\left(1-{S}_{A}\right)\right]dt+{\sigma }_{A}d{W}_{1} \\ dB=[B \left({r}_{B}-{m}_{B}B-{c}_{AB }A+{\alpha }_{B}{S}_{a}+{\beta }_{B}\left(1-{S}_{A}\right)\right]dt+{\sigma }_{B}d{W}_{2} \\ d{S}_{A}=\left[{S}_{A}\left(1-{S}_{A}\right)\left[\frac{A}{A+B}-v\frac{B}{A+B}\right]\right]dt+ {\sigma }_{S}d{W}_{3} \end{array}\right.$$

This model can be made more realistic by introducing correlations in the source of environmental stochasticity, and so to the responses of plants and soil to environmental fluctuations, which are most likely correlated. For example, a flood event will modify water availability, but also produce physical disturbance and produce anoxic conditions, at the same time. The individual effects of all these factors on the different plants and soil may differ, but the resulting responses of the state variables might be correlated. The more general model thus accounts for correlation in the three stochastic terms via a covariance matrix

$${\text{K}}_{ij}=\left[\begin{array}{ccc}{{\sigma }_{A}}^{2}& {\sigma }_{A}{\sigma }_{B}{\rho }_{BA}& {\sigma }_{A}{\sigma }_{S}{\rho }_{SA}\\ {\sigma }_{A}{\sigma }_{B}{\rho }_{BA}& {{\sigma }_{B}}^{2}& {\sigma }_{S}{\sigma }_{B}{\rho }_{SB}\\ {\sigma }_{A}{\sigma }_{S}{\rho }_{SA}& {\sigma }_{S}{\sigma }_{B}{\rho }_{SB}& {{\sigma }_{s}}^{2}\end{array}\right]$$

where the variance terms *σ**2* for plant A and B and soil *S* reflects the strength of environmental stochasticity, and the terms *ρ* are the correlations. Note that this matrix is symmetric, simply because the covariance between variable *i* and *j* is the same as the covariance between *j* and *i.* We used Cholesky factorization (Allen 2010; Kloeden et al. 2012) of the covariance matrix K*ij* to define the diffusion matrix that accounted for correlated environmental fluctuations in our matrix. The diffusion matrix now becomes:\({\Sigma }=\left[\begin{array}{ccc}{\sigma }_{A}& 0& 0\\ {\sigma }_{B}{\rho }_{BA}& {\sigma }_{B}\sqrt{1-{{\rho }_{BA}}^{2}}& 0\\ {\sigma }_{S}{\rho }_{SA}& \frac{{\sigma }_{S}({\rho }_{BS}-{\rho }_{AS})}{\sqrt{1-{{\rho }_{BA}}^{2}}}& {\sigma }_{s}\sqrt{1-{{\rho }_{AS}}^{2}-\frac{{({\rho }_{BS}-{\rho }_{AS})}^{2}}{1-{{\rho }_{BA}}^{2}}}\end{array}\right]\)

Note that if the correlation terms *ρ* are set to zero, the matrix returns to the special case of independent Wiener processes, with variance terms in the diagonal equal to the three *σ* terms for *A*, *B* and *S* and all the off-diagonal elements equal to zero.

## Computation

The key goal of the analysis of a stochastic process is the characterization of the probability distribution that describes how the population changes with time (Karlin and Taylor 1975). The variation of the population is not deterministic, meaning that given the same initial conditions, the population can take many different paths or trajectories. The population of size *N* is thus best described by a probability distribution *P**(N, t)* that quantifies the key features of the collective behaviour of the paths. The mean and variance of *N* over time *t* are some of these key features. Another important feature is the probability that the population goes below a certain value, for example below zero (i.e. extinction). Or also, the probability for the population to be at a certain value (for example two times the initial population size, or zero) after a certain amount of time (i.e. hitting time). For simple processes, the probability distribution *P**(N, t)* has an analytically close form that can be derived by the Fokker-Planck equation that corresponds to the Itô SDE. Our system, however, has a complex form, and we could not resolve it analytically. We thus simulated it numerically, using the Euler-Maruyama method (Allen 2010; Dobrow 2016), and we provide the full R script that implemented this method for our system (RScript_Bever_Model_Stochastic.R).

We analysed multiple scenarios, starting from the same parameters as in (Bever 2003), apart from equating the carrying capacity of the two plant species (which we did just for convenience, and with no effect on the final results). First, we compared a scenario in which the stochastic components of the process are not correlated with a scenario in which they are correlated and, within each of these two scenarios, we also compared a scenario where the variances of plant A and B were relatively high with a scenario in which the variances where low (results in Fig. 2). Also, we compared a scenario in which only the variance of A was high (or low) with the opposite scenario (results in Fig. 3). Finally, we also investigated how variance controls the probability of extinction at different level of soil feedback *I**s* (results of Fig. 4). We replicated the simulation of each scenario 1000 times, and estimated the probability of a certain outcome (e.g. unfeasible system, that is at least one species goes extinct) as the frequency of that outcome over 1000 independent replicates of the same process. In the scenarios in which we changed the soil feedback *I**s* = *α**a* - *α**b* - *β**a* + *β**b*, we set to negative the PSF by keeping three parameters (*α**b*, *β**a*, *β**b*) and decreasing the other parameters (*α**a*). But PSF being the same, we expected that the probability of extinction changed with different combination of parameters. We thus also explored scenarios in which we kept the PSF positive, constant, and compatible with coexistence at low level of variance, and then progressively increased the difference between parameters *α**b* and *β**a*.