Demographic noise in complex ecological communities

We introduce an individual-based model of a complex ecological community with random interactions. The model contains a large number of species, each with a finite population of individuals, subject to discrete reproduction and death events. The interaction coefficients determining the rates of these events is chosen from an ensemble of random matrices, and is kept fixed in time. The set-up is such that the model reduces to the known generalised Lotka-Volterra equations with random interaction coefficients in the limit of an infinite population for each species. Demographic noise in the individual-based model means that species which would survive in the Lotka-Volterra model can become extinct. These noise-driven extinctions are the focus of the paper. We find that, for increasing complexity of interactions, ecological communities generally become less prone to extinctions induced by demographic noise. An exception are systems composed entirely of predator-prey pairs. These systems are known to be stable in deterministic Lotka-Volterra models with random interactions, but, as we show, they are nevertheless particularly vulnerable to fluctuations.


Introduction
Ecological communities are inherently complex systems, in which species interact with each other, generating intricate equilibria, patterns and dynamics. The factors that promote or hinder stability are of great interest, together with the question of how much complexity an ecosystem can sustain at stable equilibria. This is at the heart of the so-called 'stability-diversity' or 'stability-complexity' debate, started fifty years ago by Robert May [1,2,3,4]. May started from a hypothetical equilibrium of an ecosystem with many species, and then modelled the Jacobian at this equilibrium as a random matrix (we will refer to ecosystems with random interactions as 'complex'). He then used knowledge from random matrix theory (i.e., the eigenvalue spectra of random matrices) to conclude that an increasing number of species ought to make an ecological equilibrium unstable in random ecosystems. This lead him to speculate that there must be some 'devious strategy' nature is using to maintain complex equilibria [2]. May's model has been subject to a variety of criticism (for an overview see e.g. [5]), and has been extended and refined in a number of different ways. This includes for example generalisations to spatial systems or to more structured random matrices [6,7,8].
Much of this work does not specify an actual dynamics for the ecosystem. A parallel stream of work therefore focuses on generalised Lotka-Volterra equations or related differential equations, with a random interaction structure similar to the hypothetical Jacobian matrices in May's approach [9,10,11,12,13,14,15,16,17,18,19]. The analysis of such systems is not straightforward, the nonlinearity of the dynamics combined with the random interactions make an analytical solution for individual instances impossible. The problem can however be recognised as analogous to questions in the theory of spin glasses and disordered systems, originally devised to describe certain materials at low temperatures [20]. Analytical tools are available from statistical physics to address dynamical systems with quenched (i.e., time-independent) random interaction coefficients [20,21,22,23]. These methods were first applied to random-replicator models in [9,10], these are systems that are structurally very similar to generalised Lotka-Volterra equations (gLVE). Subsequent work addresses gLVEs directly, and has established when such systems are stable or unstable respectively, and characterised the equilibria in the stable phase. This work is based on generating functionals [16] or the so-called cavity method [13,14,15]. For systems in which the interactions between species are symmetric the celebrated replica method can also be applied [18,19].
The starting point of these existing analyses is set by the random replicator or generalised Lotka-Volterra dynamics. While the interaction coefficients are drawn at random at the beginning the subsequent dynamics is assumed to be deterministic. Analytical work then focuses on describing where in parameter space stable fixed points are found, and what the statistics of these fixed points are. Crucially though, the effects of so-called demographic (or intrinsic) noise are often left out. The deterministic gLVE describe continuous abundances of species, and discard the inherent randomness of birth and death events. Neglecting this stochasticity is valid when the number of individuals for each species is large, any death or birth then only changes abundances by a small amount, and the stochasticity in the birth and death events 'averages out'. However, complex interaction models are also used to describe microbial populations [24,25] and in game theory [26,27], and it is known that demographic noise can significantly alter the behaviour of finite evolving populations in these areas [28,29]. In this work we therefore set out to study some of the effects of demographic (or intrinsic) noise on the dynamics of complex Lotka-Volterra dynamics. We note that the work in [18,19] makes first steps towards including demographic noise in complex Lotka-Volterra systems. However intrinsic stochasticity is there modelled by adding a Gaussian noise term to the gLVE.
In the present work we take a different approach. We introduce an individual-based Lotka-Volterra system, describing a population in which there are finitely many discrete individuals for each species. These individuals undergo discrete birth and death events, with rates governed by their interactions. The model is set up such that there is a control parameter setting the typical number of individuals belonging to a species (the parameter can be thought of the area or volume in which the ecosystem resides). In this way we recover the known random gLVE by sending this parameter to infinity. Our work focuses specifically on the effects of intrinsic noise on the stable fixed points of the limiting deterministic system. We study how this stochasticity can lead to additional extinctions and how the distribution of species abundances changes as a consequence. We ask what species in particular are at risk of such noise-driven extinctions. We note that individual-based models or other stochastic dynamics have been used for Lotka-Volterra dynamics for example in [30,32,33,34,35,36,37,38,39]. With the exception of [30] this previous literature does not focus on complex ecosystems however (that is ecosystems with random interaction coefficients). Our model contains two types of stochasticity. One is the quenched randomness of the interactions which chosen before the dynamics starts, but then remain fixed. The other type of noise describes randomness as the dynamics of the system unfolds. The combination of the two is where the main novelty of our work lies.
The remainder of the paper is organised as follows: In Sec. 2 we briefly summarise the deterministic random gLVE equations, and then introduce the individual-based model and describe how it can be simulated. To set the scene, Sec. 3 then contains a brief account of the known phenomenology of the deterministic gLVE model. This serves as a basis to compare the individual-based model against. Our main results are the contained in Sec. 4. We investigate the effects of demographic noise and show how stochasticity changes species abundances relative to the deterministic model. Sec. 5 finally contains our conclusions and a discussion of our results.

Model definitions
In this section, we first summarise the conventional deterministic random Lotka-Volterra model, and then define an individual-based model of these dynamics.

Deterministic generalised Lotka-Volterra dynamics with random interactions
The model describes N interacting species, labelled i = 1, ..., N . The abundance of species i at time t is denoted by x i (t), we collect the abundances into the vector x = (x 1 , . . . , x N ). The generalised Lotka-Volterra equations are then The coefficients α ij describe the interactions between the species, representing the benefit or detriment that species i receives from species j. Consequently, a negative coefficient α ij indicates a detrimental effect of species j on species i. If α ij and α ji have opposite signs, then species i and j constitute a prey-predator pair (i.e., there is antagonistic interaction). If α ij and α ji are both negative, the two species compete with each other, and if both coefficients are positive then they have a mutually beneficial relation with one another. The quantity K i in Eq. (1) is the carrying capacity of species i in monoculture, this is the equilibrium value of the abundance of species i if there is no interaction with any other species (α ij = 0 ∀j = i). In principle heterogeneous K i can be considered, however we focus on the simpler case K i ≡ 1 for all i.
Following [13,14,18,16] we choose the interaction coefficients from a fixed distribution at the beginning, they then remain constant as the dynamics proceeds. More precisely, the off-diagonal coefficients α ij are drawn from a Gaussian distribution with the following mean, variance and correlations, where the overbar denotes an average over the ensemble of random interaction matrices. The scaling with N is chosen to guarantee a well-defined thermodynamic limit N → ∞ [20]. In our notation the diagonal coefficients α ii do not enter in Eq. (1), but we note the stabilising term −x i inside the bracket on the right-hand side, which could be absorbed in the interaction term by defining α ii = −1.
A positive value of µ means that species tend to benefit from each others presence. On the other hand, a negative µ indicates competition. The parameter σ 2 controls the variance of the interactions, and describes the degree of heterogeneity or complexity of the species interactions. Finally the symmetry parameter Γ accounts for the correlations between the interaction coefficients of any pair of species, and can take values between −1 ≤ Γ ≤ 1.
In the limit N → ∞, the fraction of predator-pray pairs in the system (i.e, pairs i, j with α ij α ji < 0) can be obtained performing a Gaussian integral over the two quadrants, {α ij ≥ 0, α ji ≤ 0} and {α ij ≤ 0, α ji ≥ 0}, of the joint distribution of α ij and α ji (see e.g. [40]). This results in a non-linear and decreasing dependence of the fraction of predatorpray pairs p on Γ in the limit of large N , p = 1 2 − 1 π sin −1 (Γ). In particular, for Γ = 0, α ij and α ji are uncorrelated and there are 50% of predator-prey pairs in the system. If Γ = 1, then α ij = α ji (fully correlated) and there are no predator-prey pairs. Finally, for Γ = −1 (full anti-correlation), all pairs in the system are of the predator-prey type.
Despite the fact that the interaction coefficients are chosen at random at the beginning, we will refer to in Eq. (1) as deterministic. This is because there is no further stochasticity during the actual dynamics.

Individual-based model
We now define an individual-based variant of the generalised random Lotka-Volterra dynamics.
2.2.1. Birth-death dynamics. In order to incorporate demographic stochasticity into the deterministic Lotka-Volterra model, we focus on a population of discrete individuals and consider the dynamics as a birth-death Markov process in continuous time, driven by discrete events. In each event, an individual of a species i reproduces or dies, increasing or reducing respectively the abundance of that species. The rates for these events are T + i and T − i respectively. Therefore, at each time step there are 2N possible reactions or events. Similar approaches have been used in [30,32,33] and especially in [41].
We write n i (t) for the number of individuals of species i at time t, and n(t) = [n 1 (t), . . . , n N (t)]. In order to define the rates, T ± i , we first introduce a parameter V , characterising the area or the volume in which the population resides. We will chose the birth and death rates such that the deterministic gLVE model is recovered in the limit V → ∞. We write x i (t) = n i (t)/V . The deterministic Lotka-Volterra equations (1) (with K i ≡ 1) can then be written as where we have separated the terms with positive α ij from those with negative interaction coefficients. These terms set the rates at which individuals of species i reproduce and die in the stochastic model. We define, There is therefore a baseline per capita reproduction rate of one [first term in T + i (n)], and reproduction rates are enhanced by positive interaction with other species. The first term in the death rate T − i (n) describes intra-species competition, and the second term (increasing the death rate) represents detrimental effects on species i by other species. Given that the n i scale linearly in V , the reaction rates are of order O(V ), indicating that the total number of events in the system per unit time is proportional to V .
The stochastic process of n defined by the rates T ± i (n) can formally be described by a master equation. Using standard methods [42,43,44] one can then show that the moments n i of the process follow the Lotka-Volterra dynamics in Eq. (3) in the limit of large V (in this limit, the distribution of n at time t becomes sharply peaked around its mean). In other words, the quantities x i = lim V →∞ n i /V follow Eq. (1).
We stress that the choice of rates T ± i (n) in Eq. (4) is not the only one leading to Eq. (1) in the deterministic limit. For example, one could consider combined birthdeath events, in which an individual of one species is replaced by that of another. We do however expect the principal effects we observe in the stochastic model to be robust against variation of these details.

Simulation method.
We use the well-known exact Gillespie algorithm to simulate the invididual-based model [45,46]. This generates a statistically faithful ensemble of trajectories of the stochastic process. The main principles of the simulation are an exponentially distributed time increment until the next event, followed by a stochastic selection of the type of event. Both the parameter of the exponential and the probabilities with which the different possible events are chosen are set by the T ± i (n). In detail, the simulation proceeds as follows: 1. Draw a random interaction matrix α ij with the statistics in Eq. (2). Set the initial conditions for the abundances n i (i = 1, . . . , N ) at t = 0.

Background: behaviour of the deterministic model
In this section we briefly describe the known phenomenology of the deterministic Lotka-Volterra system with random interactions in Eq. (1). This serves as a baseline for our study of effects induced by demographic noise in Sec. 4.

Generating functional analysis
The gLVE system with random interaction matrices has been studied analytically in the limit N → ∞ with tools from the theory of disordered systems, in particular using pathintegral analyses [16,17] or the so-called cavity method [13,14,15]. Static approaches to the symmetric system (Γ = 1) use the celebrated replica method, see e.g. [18,19].
In the essence these approaches involve carrying out the disorder-average (often via a static or dynamic partition function), and then deriving an effective dynamics or an effective potential for a representative species. The fixed points of this dynamics (or the minima of the potential) can then be characterised further. Ultimately, these approaches all lead to the same results of course when they can be applied. However, generating-functional and cavity methods are more general in that they do not require any restrictions on Γ. The replica method on the other hand allows for a more in-depth analysis of the 'energy' landscape of the model and of the resulting analogies to spin glasses [18,19].
We here do not describe the analysis in detail, instead we focus on the equations for the key order parameters describing the stable fixed points of the gLVE.
The gLVE system exhibits a stable phase, in which a given realisation of the random interaction matrix leads to one unique stable fixed point with overwhelming probability as N → ∞. This phase is found at sufficiently low values of the paramaters µ and σ, see Sec. 3.2 for further details.
The key order parameters in this fixed point phase are the mean abundance per species at the fixed point, M , the second moment of the fixed point abundances, q, and the so-called susceptibility χ, which indicates how a fixed point is shifted in abundance space in response to a permanent external perturbation [13,14].
From the path-integral analysis and a subsequent fixed-point ansatz (or alternatively, from the cavity approach), one derives the following relations for these order parameters [13,14,16], From the solution of these equations, the species abundance distribution is found as where G(x) is a Gaussian distribution with mean (1 + µM )/(1 − Γσ 2 χ) and variance qσ 2 /(1 − Γσ 2 χ) 2 , clipped at x = 0 (abundances only take non-negative values). The coefficient φ = w 0 (∆) is the fraction of surviving species. Examples of these species abundance distributions are shown further below in Fig. 7 (a)-(c).

Instabilities and phase diagram
Two types of instability of the fixed points in the gLVE system have been established [13,14,47]. One is a linear instability, this sets in when Therefore for σ 2 < σ 2 c (and in absence of the other type of instability) one expects the Lotka-Volterra dynamics to have stable fixed points in the limit of large number of species, N . The fixed points are unique, i.e. independent of initial conditions, for a given interaction matrix. For σ 2 > σ 2 c , the system can show several different types of behaviour. For example, the dynamics can remain volatile and potentially chaotic, or there can be a very large number of fixed points (initial conditions then determine which one of these is reached). For further details see for example [17].
The second type of instability occurs when the mean abundance per species, M , diverges. This indicates unbounded growth of the ecosystem. No closed-form expression is available for the point in parameter space where this happens, however the onset of this instability can be obtained numerically from Eqs. (5), and imposing 1/M = 0 [13,14,16].
The resulting phase diagram in the (µ, σ 2 )-plane is shown for different values of Γ in Fig. 1. The system is stable and Eqs. (5) are valid in the region below the dashed lines (linear instability) and to the left of the solid lines (diverging abundance).

Effects of demographic noise
We now study the effects of demographic noise on the ecosystem. Our analysis focuses on parameters in the stable phase of the gLVE. The rationale for this focus is discussed in the conclusions in Sec. 5.

Effects on trajectories
We start by illustrating trajectories of the deterministic and stochastic systems in Fig. 2, for one particular realisation of the interaction matrix. Model parameters are chosen such that the deterministic system approaches a stable fixed point (indicated by the dashed lines). As seen in the figure the abundances of the stochastic system broadly fluctuate about the deterministic fixed-point values. By virtue of the central limit theorem the variance of these fluctuations can be expected to scale as V −1 .
As discussed in the previous section, some species become extinct at the deterministic fixed point. Those species are normally also found to die out in the system with demographic noise. Additionally, there can be species which survive in the absence of noise, but which are driven to extinction by fluctuations in the stochastic system. An example can be seen in the bottom right of Fig. 2.
While the deterministic system settles down to a fixed point in the stable regime, the dynamics in finite populations remains subject to fluctuations. As a consequence the mean number of surviving species in the stochastic system (averaged over realisations) shows a mild decay over time, see Fig. 3(a). There is therefore no true stationary state for the system with intrinsic noise [one can expect all species eventually to become extinct at extremely long times (scaling exponentially in V )]. For the purposes of our analysis we do however have to choose a time at which we study the stochastic system. We choose this time scale (t ≈ 100) such that the deterministic system has reached its fixed point. We have also verified that results for the stochastic system do not change qualitatively if a longer time scale is chosen (e.g. t ≈ 200). Further details will be described below.
We first study the magnitude of fluctuations more systematically. To this end we define (for a single realisation), where · · · t stands for an average over time after some equilibration period. More precisely we calculate this average in the stochastic system during the time interval from t = 60 to t = 100. Thus, h 2 rel describes the relative fluctuations of the species abundances in time. We also define As seen in Fig. 4(a) the absolute fluctuations of species abundances increase with σ 2 for fixed Γ > −1, showing a particularly steep rise as the linear instability in the deterministic model is approached. One possible reason for this is the growth of the species abundances when the deterministic system approaches the instability [as shown in Fig. 5(b)]. For Γ = −1 there is no instability and the predator-prey relations in the ecosystem prevents large abundances. As a result, the increase in absolute fluctuations is only moderate. (We remark that throughout the paper we often distinguish between the cases Γ = −1 and Γ > −1. This is a broad classification, we cannot exclude that the system behaves similarly to the fully anti-correlated case when Γ is near −1.) Surprisingly, relative fluctuations decrease with σ 2 for Γ > −1 [see Fig. 4(b)]. The relative effect of the intrinsic noise on the abundances is hence reduced when increasing the diversity in the interactions. However, this effect becomes less pronounced when interactions are anti-correlated (more negative values of Γ). Indeed, if all interactions are of a predator-prey type (Γ = −1) relative fluctuations around found to increase only mildly with σ 2 .

Community properties
In Fig. 5 we study how the size of the surviving community (indicated by the fraction of survivors, φ) and the mean abundance per species, M , vary with σ 2 , for different values of the correlation parameter Γ. We note again that these quantities vary in time in the  stochastic model (e.g. φ decreases with time). The results in the figure are therefore only valid at the time they were measured (t = 100). We have however checked that the qualitative behaviour as a function of σ 2 and Γ does not change at later times.
Unsurprisingly, the fraction of survivors is systematically lower in the stochastic model (markers) than in the deterministic one (solid lines). This is because fluctuations cause some species to eventually become extinct in the stochastic model while the same species survive in the absence of demographic noise. The average abundance M is measured across all species, including those that go extinct. As a consequence M is also lower in the stochastic model than in the deterministic gLVE.
The data for the fraction of surviving species in Fig. 5(a) indicates that the difference between the communities resulting from the stochastic model and that of the deterministic dynamics becomes more pronounced as the fraction of predator-prey pairs in the interactions increases (i.e., as Γ is lowered to approach Γ = −1). One possible explanation is that, although communities in the deterministic model are more stable for predator-prey interactions (and there are fewer deterministic extinctions), species abundances are generally smaller [ Fig. 5(b)]. Therefore extinctions caused by fluctuations in the stochastic model become more prevalent as Γ becomes more negative.
In Fig. 6 we analyse in more detail how the difference between the size of the surviving communities in the two models, measured by |φ det − φ stoch |, changes with the complexity of the interaction strengths σ 2 . We observe virtually no difference in the number of survivors in the absence of interactions (σ 2 = 0). At the fixed point of the deterministic system all species abundances then take the value x i = 1, and extinctions in the stochastic model are rare. As σ 2 is increased from zero we observe an initial increase of |φ det − φ stoch | with σ 2 for all values of the correlation parameter Γ that we have tested.
For fully anti-correlated interactions, Γ = −1, this increase continues to larger values of σ 2 . In this situation, the relative magnitude of fluctuations of abundances in the stochastic model grows with σ 2 [Fig. 4(b)]. This leads to an increased number of extinctions in the stochastic model, caused by intrinsic noise.
When interactions are not fully anti-correlated (Γ > −1), the difference |φ det −φ stoch | becomes smaller as σ 2 increases (except for the region near σ 2 ≈ 0 discussed above). Hence the deterministic and stochastic models become more similar to one another, and extinctions occur mostly due to interactions between species (as opposed to demographic noise). This is in line with our earlier observation that relative fluctuations become weaker with increasing σ 2 when Γ > −1 [ Fig. 4(b)], and that the mean abundance per species M increases [ Fig. 5(b)]. Extinctions due to fluctuations hence become more rare with increased variance of the interactions.

Species abundance distribution
We next look at the effects of intrinsic noise on species abundance distributions (SAD). In deterministic model the mean abundance increases with σ 2 , and the SAD becomes wider. Species with low abundance in the deterministic model are generally more likely to go extinct in the face of noise than those with high abundance. So we expect that any deviation of the SAD in the stochastic model from that of the determinisitic model will result mostly at small abundances x. This is confirmed in Fig. 7.
Broadly, the picture is here that the gap between the SADs of the deterministic and stochastic models is localised mostly at x < x c , with x c some critical abundance, and that for abundances x > x c the SAD is mostly unaffected by intrinsic fluctuations. (This is not a strict classification, we are not suggesting that the SADs of the two models are identical for x > x c . But we think the notion of a value of x below which most of the deviations are seen is useful.) The numerical value of x c will depend on the model parameters. In particular, the data in Fig. 7 suggests that x c grows with σ 2 . This is because absolute fluctuations of species abundances become larger with σ 2 , as previously shown in Fig. 4(a). Consequently, as σ 2 is increased, extinctions due to demographic noise can occur for species with higher abundances.
In Fig. 8, we show the probability for a species to go extinct if it has attained an abundance x in the deterministic model. To measure this, we first fixed an interaction matrix, integrated the gLVE numerically, and determined the abundances of the surviving species in the resulting equilibrium. We then ran multiple simulations of the stochastic model, with the same interaction matrix, and determined how likely the different species were to survive up to time t = 100. This was then repeated for different realisations of the interactions, for a given choice of σ 2 and Γ. We observe in Fig. 8 that species with higher abundance in the deterministic model also become prone to extinction as σ 2 increases. Further, the data demonstrates that making the interactions more anti-correlated (increasingly negative values of the correlation parameter Γ), the probability to go extinct for species with high to moderate x decreases [at fixed σ 2 , compare panels (a), (d), (g) in Fig. 8]. Indeed, reducing Γ at fixed σ 2 reduces absolute fluctuations as previously shown in Fig. 4(a).
We note that the extinction probabilities in Fig. 8 are, in isolation, not immediately a statement about the total number of extinctions in the system due to fluctuations. What is shown in Fig. 8 is the probability that a species becomes extinct in the stochastic system provided it attains a given abundance in the deterministic system. The total number of extinctions due to noise results from a combination of these probabilities and the abundance distribution in the deterministic system. For Γ > −1, these abundances move towards higher values as the variance of interactions, σ 2 , is increased [see Fig. 5(b)]. There are therefore two competing effects acting on total extinctions as σ 2 is increased, increasing abundances (reducing the propensity for extinction) and the increasing probability to go extinct for any fixed deterministic abundance, x.

Abundance dispersion
As a further point we now study of the differences between the abundances a species attains in the stochastic and the deterministic systems with the same interaction matrix. Generally we find that the abundances in the stochastic system, although fluctuating in time, are closely related to the deterministic fixed point abundances for systems that are far away from any instability. However, as we will discuss, this is not the case when parameters are such that the deterministic system approaches the unstable phase.
To study how abundances in the two different models differ from one another, we produce scatter plots of these two quantities. Each point on the graph is a species, and the coordinates of the points are the species' abundances in the deterministic and the stochastic systems with the same interaction matrix. If the abundances in both models are fairly similar, one would expect to see a dense region of points near the diagonal. Conversely, points away from the diagonal indicate differences between the abundances in the stochastic and deterministic models. We use the temporal average value of the abundances measured in the steady state of the deterministic system. We average the abundances in the stochastic system from t = 60 to t = 100. For the deterministic abundances, we average from t = 90 to t = 100, this is sufficient to capture the deterministic fixed point of the dynamics. For all simulations we use N = 50, V = 100, and µ = −1. Data is shown in Fig. 9. Each of the scatter plots in the different panels is from an aggregate of data from 10 different interaction matrices.
The figure demonstrates that the abundances of the two models are quite similar for small σ 2 , but differ more substantially when approaching the unstable regime, therefore appearing more dispersed in the graph. While we only show data for one particular choice of the correlation parameter Γ, we have tested other values as well, and found similar behaviour.
To capture these findings more compactly we define the average absolute difference between abundances as where we restrict the sum over i to species that survive both in the deterministic and in the stochastic model (as indicated by the set S). Numerically, we identify these as species with x i det t > 0.001 and x i stoch t > 0. The normalising factor N S in Eq. (10) is the number of such surviving species. Including only the set of survivors avoids that the outcome is dominated by species that go extinct in the stochastic model, but survive in the absence of noise.
Indeed, as shown in Fig. 10, δx grows with σ 2 , with quicker increase for positive Γ and only slowly for systems with anti-correlated interactions (Γ ≈ −1), in-line with (We remark that the density of points in the plots reduces as we increase σ 2 . This is because the species abundances become larger and because we fix the range of abundances in the different panels for better comparison.) the behaviour of the absolute fluctuations of abundances shown in Fig. 4 (a). Similar behaviour is found if all species are included.

Summary of stochastic effects
We find that the fraction of surviving species, φ, is systematically lower in the individualbased model than in the deterministic gLVE. This is not surprising, and mainly serves to confirm that additional extinctions occur in the stochastic model, driven by demographic noise. The mean abundance per species, M , is also found to be lower in the stochastic model.
We also find, that at a fixed value of the correlation parameter Γ, absolute fluctuations of species abundances in the individual-based model grow with the variance σ 2 of the interaction coefficients. Hence we observe a more dispersed pattern in scatter plots of abundances in the stochastic models versus those in the absence of noise. Differences between the species abundance distributions of the deterministic and stochastic models extend to larger and larger abundances as the complexity of interactions is increased and absolute fluctuations grow. Species with higher abundances in the deterministic model hence become more vulnerable to fluctuations when the complexity of interactions is high.
It is known that abundances in deterministic model grow with the complexity σ 2 . This growth, combined with that of absolute fluctuations, creates competing forces on the relative magnitude of abundance fluctuations. We find that relative fluctuations decrease as the interactions become more varied when the system is not primarily composed of predator-prey pairs (Γ > −1). The total number of extinctions due to noise then reduces with σ 2 , despite the fact that species with higher deterministic abundances become affected by noise. This is because the increase of abundances in the deterministic model with σ 2 more than compensates for the increase in absolute fluctuations.
For systems composed only of predator-prey pairs (Γ = −1), relative abundance fluctuations increase with complexity. The number of noise-driven extinctions then also grows. The difference between the communities in the stochastic and determinstic systems becomes larger as the proportion of predator-prey pairs increases. Systems with Γ = −1 are more stable deterministically, and overall there are fewer extinctions, but there are also more species with lower abundances, hence the community of species is more vulnerable to demographic noise.

Conclusions and discussion
In summary we have used an individual-based model to study the effects of demographic noise in Lotka-Volterra systems with complex random interactions. As we have shown, intrinsic stochasticity leads to additional extinctions, and changes the composition of the surviving community. Species which have a low abundance in the absence of noise are particularly prone to the effects of fluctuations. How strong these effects are depends on the detailed circumstances, in particular on the variance of the interactions and on the number of predator-prey relations in the interaction matrix.
Our work complements existing literature approaching complex ecological communities with deterministic differential equations (with random, but fixed interaction coefficients), and a large body of work assessing stability of equilibria based on the spectra of random matrices. Neither of these approaches explicitly account for granular populations of individuals, subject to discrete birth and death events. The effects of the resulting demographic noise are particularly relevant when the number of individuals per species is small enough that abundances cannot be approximated by continuous variables. This may be relevant for example in the case of microbial communities or for models in game theory, where effects of intrinsic noise are known to be able to significantly change the outcome.
There are several limitations to our work. First, we focus on parameter choices for which the deterministic gLVE system has a unique stable fixed point. In parameter regimes where abundances diverge in the deterministic system the number of Gillespie steps required per unit of simulated time of the individual-based model also diverges. This makes simulations up to a sensible time horizon difficult. At the same time we expect that the relative effect of intrinsic fluctuations is then also rather limited, similar to what we have seen in Fig. 4(b). One aspect one could consider in future work is to characterise the effects of stochasticity just above the linear instability lines in Fig. 1, when the divergence of abundances has not yet set in [13,14].
Secondly, our study is based on simulations, no analytical theory is available at present to systematically study extinctions due to intrinsic stochasticity. Finally, we restrict our analysis to relatively simple interaction matrices, with no particular substructure. There are therefore multiple avenues for future work. This can include expansions of the Kramers-Moyal or van-Kampen type (in the limit of large, but finite V ), which would lead to sets of stochastic differential equations which capture the two types of randomness -that of the interaction matrix, and intrinsic noise in finite populations. The latter would be approximated as dynamic Gaussian in those series expansions. Initial steps have been taken in [30]. It would be interesting to see if further theoretical analysis and characterisation of the intrinsic fluctuations is possible, potentially making further simplifications such as the linear-noise approximation. One could also attempt to connect such approaches with existing replica analyses at nonzero temperature [19]. Finally, our work can be extended to more intricate interaction structures (e.g. those in [6]). The role of competition, mutualism, and nestedness could be studied, and it would also be interesting to see what the effects of intrinsic noise are in complex ecosystems with spatial structure. Further work could also focus on the possibility of noise-induced patterns in individual-based models of spatial complex ecosystems [7,8].