The persistence of bipartite ecological communities with Lotka–Volterra dynamics

The assembly and persistence of ecological communities can be understood as the result of the interaction and migration of species. Here we study a single community subject to migration from a species pool in which inter-specific interactions are organised according to a bipartite network. Considering the dynamics of species abundances to be governed by generalised Lotka–Volterra equations, we extend work on unipartite networks to we derive exact results for the phase diagram of this model. Focusing on antagonistic interactions, we describe factors that influence the persistence of the two guilds, locate transitions to multiple-attractor and unbounded phases, as well as identifying a region of parameter space in which consumers are essentially absent in the local community.


I. INTRODUCTION
Understanding patterns in the composition of ecological communities is one of the fundamental goals in ecology [1][2][3][4][5][6].A popular modelling framework for this problem considers in detail a single community in a local habitat embedded within a wider ecosystem portrayed as a species pool from which the local community can be invaded [7].For mathematical analysis, this is then supplemented with several further elements [5].The first is a dynamical model of the species abundances, and here a generalised Lotka-Volterra approach is typical but other possibilities exist [8,9].The second is a model of interspecific interactions, and in this regard a random-matrix model is often used, having a long history of shedding light on ecological questions [1,10] as well as acting as a baseline scenario against which more detailed and ecologically-motivated studies can be compared [11].
Despite the complexity of the resultant community-assembly model, analytical progress has been made [5,12,13].In particular, the dynamical cavity method (DCM), a method originating in the physics of disordered systems [14] but since adapted to a number of ecological problems [15][16][17][18][19][20][21][22], has brought significant insight into this model [23][24][25].Ref. 23 has given a comprehensive analysis of the phase diagram of this model and show that the DCM solution corresponds to a unique-fixed-point (UFP) phase, in which there exists a unique persistent community that is resistant to invasion.The DCM solution also gives the boundaries to multiple-attractor (MA), and unbounded phases.
As is typical for the DCM, the model of interactions studied in the above works is statistically homogeneous, i.e. the interaction between all species in the model is described by a single random matrix.As such, there is no a priori differentiation between the species, and the result is a single abundance distribution for the entire community.However, we know from the study of ecological networks [26][27][28][29] that interspecific interactions are anything but homogeneous and that ecological networks possess significant structure, such as trophic levels [30], nestedness [31], or modularity [32,33].One of the most common structures encountered in the network representation of ecological communities is that of the bipar-tite network, which depicts the interactions between two groups or guilds of species.The interactions described in these bipartite networks are typically either mutualistic, such as in plant-pollinator networks [34][35][36], or antagonistic such as in host-parasitoid [37][38][39][40] or trophic networks [37,[41][42][43].Bipartite networks are also found as natural components of multi-partite networks, e.g.Ref. 44 and 45, or multitrophic foodwebs e,g, Ref. 46.
In this paper, we apply the DCM to a species pool and hence local community in which the interactions are structured as a bipartite network but are otherwise random.We show the applicability of the DCM to this kind of structured scenario and derive analytic results for abundance distributions and persistence probabilities of each of the two guilds.Focusing on trophic bipartite networks, we describe the phase diagram of the consumer-resource community and show that the UFP, MA and unbounded phases of the unstructured, unipartite model still occur, but with phase boundaries that exhibit non-trivial scaling behaviour.Furthermore, we report the existence of a region in parameter space in which consumers are effectively absent from the persistent community.

II. BIPARTITE COMMUNITY-ASSEMBLY MODEL
Our species pool consists of two guilds of species in a bipartite ecological network, i. e. with interactions only occurring between species in different guilds.Let S (1) and S (2) be the number of species of each guild in the species pool, S = S (1) + S (2) be the total species number, and ρ (i) = S (i) /S be the corresponding ratios.We define N (i) α as the abundance of species α in guild i in the local community, and r α as its growth rate and carrying capacity respectively.These letter two quantities we define as positive, with their sign given by coefficient t (i) ∈ {−1, +1} that depends on whether the species in guild i grow or die out in absence of interaction.We then posit that dynamics of the abundances is described by the generalised Lotka-Volterra equations for species 1 ≤ α ≤ S (i) and guild i = 1, 2 where we adopt a periodic labelling convention that maps i = 3 onto guild i = 1.In Eq. (1), coefficients a represent the strength of interaction experienced by species α in guild i due to species β in guild i + 1.We arrange these interaction elements into the S (i) × S (i+1) matrices A (i,i+1) .We take all interactions between guild i and i+1 to be of the same the type (antagonistic, mutualistic or competitive), and correspondingly set the matrix elements as non-negative a (i,i+1) α,β ≥ 0, with the signs of the interaction provided by the interguild interaction signs c (i) ∈ {−1, +1}.
We then set the scaling of the matrix elements of A (i,i+1) with pool size such that ) are fixed as S → ∞ with ρ (i) constant.This is a natural choice when a are chosen from a non-negative distribution such as the half-normal distribution for which the standard deviation is proportional to the mean [a similar scaling was adopted in Ref. 47].We can thus rewrite matrices A (i,i+1) in terms of the centered, normalised matrix B (i,i+1) with elements We express correlations in the different interaction directions by choosing B (i,i+1) such that they have the property: with parameter γ ∈ [0, 1].We restrict ourselves to positive correlations here to avoid any ambiguity in the implied sign assignments of a . In these terms, the interaction blocks become: where J (i,i+1) is the S (i) × S (i+1) matrix of ones.Without loss of generality, we set the mean carrying capacity to be one, K (i) = 1, and then further parameterise the carryingcapacities such that κ (i) 2 = S (i) • Var K (i) is also fixed as S → ∞.This choice is justified a posteriori as being consistent with the interaction scaling.An alternate to this scaling scheme is discussed later.
In the following, the main quantity of interest will be the fraction of pool species in guild i that persist in equilibrium in which Θ is the Heaviside function.

III. DYNAMICAL CAVITY METHOD
For an overview of the DCM in the ecological unipartite context, we refer the reader to the tutorial article of Ref. 24, as well as to the work of Ref. 23.A full account of our derivation of the bipartite case is given in Appendix A, but the essence of the method is that an equilibrium configuration is considered to which a new species from each guild is added.The action of the pre-existing community on the added species is treated exactly, but the reciprocal action of the added species on the community is small in the large-S limit and treated in linear response.Since the added species are identical with other species from the same guild, this leads to a closed system of equations that can be solved self-consistently.A discussion of our solution technique for the equations is given in Appendix B and a discussion of the validity given in Appendix C.
The central result in this analysis is that the abundances of species within a guild are each distributed according to truncated Gaussians [5,12,13,23,25] with interdependent parameters.The key properties of the distributions are described by two quantities ∆ (i) ; i = 1, 2. In particular, the fraction of species in guild i that persist in equilibrium is given by in which w k are a set of functions defined via In Appendix A we give details of the equations that determine ∆ (i) in the most general case.
A. Phase diagram Our main focus will be on antagonistic interactions and for concreteness we will use the language of trophic interactions In this setting, we identify guild 1 with the resource species and guild 2 with consumers.The corresponding choices of sign are t (1) = −t (2) = +1, such that in the absence of interactions, the resource-species abundances grow to carrying capacity and the consumers die out; and c (1) = −c (2) = −1 such that the interaction is beneficial to the consumers and detrimental to the resources.
As in the unipartite case, the bipartite model is found to exhibit three phases: the UFP described by the cavity solution outlined above, plus the unbounded and MA phases.In the unbounded phase, one or more of the species abundances diverge, such that one or both of lim t→∞ N (i) → ∞.For the abundances in the cavity solution to remain bounded, we require that both sides of Eq. ( 8) (in the κ (i) = γ = 0 case) are greater than zero.As shown in Appendix D, the implication of this is that the location of phase boundary is asymptotically given by which holds for γ = 0 but arbitrary κ (i) .As example, let us assume that the interaction strengths are distributed according to a half-normal distribution for which mean and standard deviation are related as µ (i) = σ (i) 2/(π − 2).Then, assuming that the couplings are symmetric, σ (1) = σ (2) = σ, we find that the unbounded phase occurs when σ > σ c with critical interaction strength The key observation is that because 3 for large |∆|, the critical interaction strength σ c diverges exponentially with pool size S and becomes inaccessible.At finite S, the transition occurs at finite interaction strength.
Stability analysis of the cavity solution shows that, as in the unipartite case, it becomes unstable and gives way to an MA phase.Appendix E shows that the boundary to the MA phase obeys the equation (valid for γ = 0, arbitrary κ (i) ) To find the parameters of this boundary in the κ (i) = 0 case, we look for overlap of this curve with that described by Eq. ( 7).Since both curves as symmetric with respect to interchange of ∆ (i) , the boundary behaviour where the two curves just cease to overlap occurs when ∆ (1) = ∆ (2) = ∆.From this, we determine that the critical parameters occur when w 0 (∆) = w 2 (∆).This happens at ∆ = 0, at which point w 0 (∆) = 1 2 .The result is that the critical parameters for the MA transition obey For σ (1) = σ (2) , the MA phase therefore occurs when σ > σ MA = 4/(ρ (1) ρ (2) ) 1/4 .

IV. RESULTS
In visualising the results of this calculations we reduce the number of independent parameters by choosing σ (1) = σ (2) = σ, and κ (1) = κ (2) = κ.Furthermore, although µ (i) and σ (i) are in general independent parameters, here we consider them to be derived from a half-normal distribution for which they are related as In Fig. 1 we plot the persistent fractions φ (i) as a function of interaction strength σ in the simplest case of κ = γ = 0. We show analytic results from the DCM described previously as well as results obtained from numerical simulations of the GLV equations (described in Appendix G).Overall agreement is good, and is seem to improve for larger values of the pool size S.As in the unipartite case, the DCM solution still gives a good account of the simulation results in the MA phase (to the right of the red dotted line in Fig. 1), despite the loss of stability of the cavity solution in this region.
Figure 2 shows that this agreement between simulations and analytics also extends to non-zero values of κ and γ, and thus that the DCM is able to address correlations in the interaction matrix and a distribution of carrying capacities in this structured context.
The inset of Fig. 1 gives numerical evidence of the location of the MA transition.Here we plot the standard deviation of the abundance values N (i) α sampled over initial conditions, and averaged over instances of interaction matrix and carrying capacities and over all species in both guilds.For presentation, this measure of the fluctuations in the final equilibrium state is normalised against the mean abundance over all runs.We see that, for an interaction strength lower than the critical value predicted by the DCM, the fluctuations are extremely small, but that around the critical interaction strength they start to rise.Clearly this is indicative of the unique fixed point below the transition giving way to the MA phase above it, as the multiple attractors will have different abundance distributions and hence finite fluctuations from instance to instance.
From Figs. 1 and 2 we see that, even though the parameters are chosen symmetrically, there are significant differences between the two persistence fractions.In particular, at smaller interaction strengths, φ (1) (resources) is larger than φ (2) (consumers) and at larger interaction, this trend is reversed.This is perhaps not surprising, given the different roles of the two guilds.However, in the asymptotic limit (S → ∞; see Appendix F) for κ = γ = 0 and symmetric parameters we find that ∆ (1) = ∆ (2) = ∆ with ∆ determined by w 2 (∆) = 2/σ 2 for any σ = 0.In this case, then, the abundance distribution of the two guilds becomes identical and the two φ (i) curves overlap (a result shown as an orange line in Fig. 1).
Apparent from these plots is that for a range of interaction strength from zero upwards, the persistence fraction of the consumers is suppressed to being close to zero, such that at φ (1) , S = 50 φ (2) , S = 50 φ (1) , S = 100 φ (2) , S = 100 φ (1) , S = 200 The fraction of persistent resources (φ (1) , dark blue) and consumers (φ (2) , light blue) as a function of interaction strength σ with parameters κ (i) = γ = 0, ρ (i) = 1  2 , and with interactions symmetric, σ (i) = σ, and drawn from a half-normal distribution.Good agreement is seen between analytic results (lines) and numerical simulations (markers) and this increases as the size of the species pool S. The red dashed line indicates the transition from a unique fixed point (UFP) to multiple-attractor (MA) phase.For small interaction strengths, the consumer fraction drops to almost zero.The orange line shows the asymptotic S → ∞ limit, in which both φ (i) are identical for these parameters.INSET: Relative fluctuations in equilibrium abundances of the numerical solution taken over realisations of the interaction matrix and carrying capacities.Marked increases occur around the critical interaction strength of the MA transition: σ = 2 for γ = κ = 0; σ = 2.39 for (γ = 0, κ = 1).For (γ = 1, κ = 0), the transition is outside this σ-range, and no increase in fluctuations in seen.
Figure 3 show different aspects of the phase diagram of the antagonistic bipartite model.We plot results for a pool size of S = 50, such that the consumer-suppression regions are easily visible on the same scale as the MA transition.Note that we do not show that unbounded transitions as, according to the arguments presented in the previous section, these occur as large values of interaction strength.
Figure 3 shows the persistent fractions φ (i) as a function of interaction strength σ and carrying capacity width κ.Increasing κ results in the persistent fractions dropping more quickly as σ increases, and also a noticeable drop in the peak number of consumers.This is a result of the wider distribution of K (i) α values giving resource species a carrying capacity closer to zero and therefore more likely to becomes locally extinct and consumers (unsigned) carrying capacities further from zero, and therefore more likely to die out rapidly.The red dashed line shows the point of the MA transition, which is seen to move to higher interaction strength with increasing κ.Finally, the green dashed line shows the boundary Eq. ( 14) of the suppressed-consumer region.This is a constant as a function of κ here, but as κ increases, the transition out of the suppression region becomes less sharp.
Figure 3(b) shows how the persistent fractions change as a function of σ and ρ (1) , the fraction of resource species in the species pool.In this case we see how the consumersuppression region depends on the composition of the species pool with the suppression becomes more extensive for smaller ρ (1) , i.e. fewer resource species in the pool.The opposite effect is also observed -for large values of the resource fraction ρ (1) , the consumer shows an extensive range of interaction strength for which φ (2) ≈ 1 and thus all consumers in the pool are supported in the local community.The MA transition line is symmetric about ρ (1) = 1/2, given the dependence on ρ (1) ρ (2) = ρ (1) (1 − ρ (1) ).
Finally, Fig. 3(c) looks at the role of the correlation parameter γ in determining community persistence.Increasing correlations (γ > 0) universally serves to increase both persistent fractions, with extensive regions where φ (i) > 0.5 for both consumers and resources.Transition to the MA phase is moved to higher σ for increasing γ and being in the MA phase is correlated with a drop in both persistent fraction values.

V. DISCUSSION
We have shown here that the DCM generalises to structured ecological models, and specifically to a bipartite structure with consumer-resource Lotka-Volterra dynamics.We have seen that this model exhibits phases analogous to those of the unipartite model, and that the DCM allows us to map the boundaries between them.The key feature of the bipartite model is the existence of two guilds, and we have seen that the composition of the persistent community in terms of these two guilds depends both on the strength of the interaction between them, the guild size in the species pool, as well as parameters such as γ and κ (i) .Interestingly, in the S → ∞ limit the ratio of the mean equilibrium abundances becomes which depends only on the ratio of mean interaction strengths.The ratio of the persistent fraction is not as simple.
One novel feature of our results is that, for weak interactions, the fraction of persistent consumers, φ (2) , is suppressed.This is ecologically reasonable, as it means that if trophic interactions are too weak, consumers can not be sustained in the community.What is perhaps surprising is that the transition to a sustained consumer presence is reasonably abrupt in parameter space.Furthermore, we have shown that the width of this suppressed region depends on the size of the species-pool, becoming narrower as S increases.For interaction strengths above the consumer onset, properties of the two guilds become similar with both φ (i) falling off with interaction strength.
The most obvious baseline for comparison is the unipartite model with γ = −1, as this has interactions arranged in consumer-resource pairs Ref. 48.The persistent fraction in that model shows a monotonic decrease with interaction strength, similar to the asymptotic behaviour of the bipartite model in Fig. 1.The unipartite model shows no suppression for small couplings because, although the interactions are all trophic, their orientations are random.In contrast, the bipartite structure enforces a consistent direction to the interactions and this paves the way for guild-level effects.
Whilst we have focused on trophic interactions, different choices of the sign factors t (i) and c (i) allow different interguild interactions to be studied.With c (i) = +1 we have a mutualistic bipartite network, either obligate (t (i) = −1) or facultative (t (i) = +1).Adapting the reasoning from the trophic case, we find that, irrespective of t (i) , the transition to the unbounded phase in this mutualistic model occurs for large ∆ (i) , and this allows to obtain the transition point as µ (1) Taking the case of symmetric coupling drawn from the half-normal distribution, the critical coupling is given by such that bounded phase is obtained when σ < σ c ∼ S −1/2 .The extent of this bounded phase therefore decreases with pool size and vanishes in the limit.And thus the characteristic behaviour of mutualistic bipartite interactions is towards non-persistence of the community.This might be reconciled with the manifold observation of bipartite mutualistic networks in nature in a number of ways.It might imply that mutualistic interactions are very weak, but this seems unlikely given the important role these interactions typically play in the lifecycles of the participants.It could also indicate limitations in the dynamical model -inclusion of saturating interactions is an obvious improvement that could be made.But perhaps the most interesting possibility is that the model suffers from studying mutualism in isolation, and in nature these mutualistic networks modules exist as modules in larger networks with interactions of various types.Finding the networks conditions which allows mutualisms to persist therefore becomes an important future question.The sign allocation t i) = +1 and c (i) = −1 gives a model in which the two guilds compete with one another.In this case the transition to unboundedness occurs when µ (1) µ (2)   (σ (1) With symmetric couplings and matrix elements from the half-normal distribution again, this translates into a critical coupling such that the bounded phase occurs for σ < σ c ∼ S 1/2 .The two guilds therefore both always persist in the large-pool limit.Although we were unable to find reports of bipartite competition networks in the literature (presumably due to the difficulty of observing such interactions), "negative non-trophic" interactions have been reported as part of larger multiinteraction networks [49,50] and often with a particular association with facilitation [51].
In both mutualistic and competitive cases, the symmetry in interaction sign, c (1) = c (2) , ensures that neither guild is significantly suppressed relative to the other.It is interesting to compare these results for the persistence of communities with different interaction types with the conclusions derived from a linear stability analysis of the relevant interaction matrices in which 1 is a unit matrix, and where we set γ = 0 for simplicity such that A (1,2) and A (2,1)   are independent.From e.g.Ref. 47 and references therein, we know that the asymptotic spectrum of A will consist of two parts: a bulk, and a pair of isolated "macroscopic" eigenvalues.In the uncorrelated case, the bulk spectrum will be a circle in the complex plane [52] with centre at (-1,0) and radius of √ σ (1) σ (2) ρ (1) ρ (2) 1/4 .The macroscopic eigenvalues are given by with the scaling ∼ √ S justifying the "macroscopic" moniker.In the consumer-resource case, c (1) c (2) = −1, and the macroscopic contribution to the spectrum is purely imaginary.Thus it is the bulk that determines the stability.In contrast, for both competitive and mutualistic interactions, we have c (1) c (2) = 1 and the stability properties of these two interaction types will be the same.In these cases we have λ macro,+ real and positive and therefore this eigenvalue dominates stability considerations.The trend, then, from local stability analysis is that both competitive and mutualistic bipartite interactions are unstable, whereas antagonistic interactions are stable.This stands in contrast with the DCM results which identifies the antagonistic and competitive structures as persisting, whilst the mutualistic one does not.
We now discuss the scaling of the interaction coefficients [53].The most obvious alternative to Eq. ( 3) is to set so that the mean scales like the variance (rather than the standard deviation).This scaling is like that in the unipartite model [23].The difference between these two choices is perhaps best appreciated from the spectrum of A of Eq. ( 19).In the alternative scaling, A has the same bulk spectrum as before but now the S-dependence of macroscopic eigenvalues removed.Thus, eigenvalues λ macro,± cease to be macroscopic, and the entire spectrum scales as S 0 .The DCM equations for this alternative scaling can be obtained from those presented here by scaling µ (i) → µ (i) S −1/2 and κ (i) → κ (i) √ S (i) .This results in a modification of the phase diagram.Confining ourselves to the trophic case, the transition to the unbounded phase now occurs at coupling strength ∼ S 0 , rather than ∼ S 1/2 , and the width of the suppressed zone becomes ∼ S 0 rather than ∼ S −1/2 .On the other hand, the MA transition remains in the same place, being determined by σ (i) and not µ (i) .Distributions that scale like Eq. ( 3) are straightforward to realise -the half-normal-distribution used here is a simple example.Not so for Eq. ( 21) as this requires a distribution defined on non-negative support a > 0 in which the ratio of mean to standard deviation scales like S −1/2 , this inevitably results in a heavily skewed distributionRef.54.Other than this seemingly extreme properties required of the distribution, the second problem with this scaling is that it limits the parameter values for which the DCM solution is accurate, since the skewness of the distribution compromises the Normality assumption (see Appendix C) unless we have µ (i) /σ (i) 1 .A second variation of the model is to drop the restriction that the interaction elements be non-negative.This obviously changes the intent of the model as, from a starting point in which all interactions are e.g. trophic interactions, negative values of a (i,j) α,β mix in some interactions that are mutualistic, some competitive, and some that remain antagonistic but opposite in direction.Nevertheless, if the majority of the interactions remain of the original type, it still makes sense to differentiate the two guilds along the original lines.Such a model might be appropriate for bipartite plant-microbe networks [55] where interactions are complicated and of different signs [56,57].Dropping this restriction does not change derivation of the DCM equations, but it does affect the validity argument presented in Appendix C. If we choice the scaled matrix elements b (i,j) α,β from a normal distribution (which necessarily permits negative values), then we remove any concerns about the Normality of the final fluctuations, and there are no limits on the validity of the DCM equations from a skewness point of view.This means that the scaling of Eq. ( 21) works just as well as the scaling of Eq. (3).However, preservation of guild interaction identity still requires µ (i) /σ (i)  1, and so the useful parameter regime of the model stays the same.Looking to the future, this work opens up the study of community assembly within other block-structured ecological networks.Here we think of two particular geometries: "hub and spokes" in which a central guild interacts with a number of further guild, as in Ref. 44, and a "ladder" such as a food web with perfect trophic coherence, i.e.where basal species are consumed exclusively by primary consumers, primary consumers are consumed exclusively by secondary consumers, and so forth [30].Tripartite ecological networks [47,[58][59][60][61] span both categories, and could either describe a single interaction type e.g.antagonism in a plant-pest-parasitoid network, or mixed interactions such as in plant-mutualist-parasitoid or herbivore-plant-mutualist networks.Following the approach set out here, the DCM should allow us to map the persistence and coexistence conditions across these diverse network structures which, although certainly just caricatures, represent important aspects of the organisation of interactions central to natural ecosystems.and, for purposes of derivation, added small perturbation ξ  α and to each guild we then add a species, labelled with index 0 and assumed to be statistically homogeneous with the other species in its guild.We denote the new equilibrium abundances as m (i) α .These abundances will obey the same equations as Eq.(A1) but with m α and with sums extended by one.This allows us to identify the perturbations as Since both m (i) 0 and b scale as S 0 , we have ξ α → 0 as S → ∞ such that the perturbations are weak.The abundances with and without the perturbation can therefore be related via linear response as where, in the last expression, we have set the off-diagonal response blocks v (i,i+1) to zero, a result which will be justified later.Under the assumption that m (i) 0 > 0, substitution of the previous result into the corresponding equilibrium condition gives 0 = λ This can be rearranged as As reasoned in Ref. 23, the denominator of this expression is a finite number with negligible fluctuations and thus be replaced by its expectation value.Evaluating the expectation value of the sum, we obtain which defines the mean response coefficient, ν (i) .The denominator of Eq. (A7) can therefore be written as Expressions for ν (i) can be obtained by differentiating Eq. (A7) and noting that we only obtain a non-zero contribution with a probability φ (i) : We also see from Eq. (A7) that ∂m = 0, consistent with the earlier assertion that v (i,i+1) α,β = 0.The final equation required to close the set is f (1) + f (2) = 1 by definition.

Appendix B: Solution of DCM equations
We can collapse the above set of equations into a pair of equations, solution of which gives parameters ∆ (i) as a function of model parameters.First we use Eq.(A19) to write Eq. (A9) as This maps onto a pair of independent quadratics with solution where have taken the positive root so that û(i) → u (i) for γ → 0.
We then turn to Eq. (A15) and use Eq.(A18) to obtain with the matrix There are then two different solution routes, depending on whether both κ (i) = 0 or not.Assuming at least one κ (i) is non-zero and that M is non-singular, we obtain a solution R (i) = gS −1/2 P (i) with where we take the positive root P (i) > 0 since R (i) > 0. This gives such that Eq. (A17) may be written as The two resultant equations can then be solved to give ∆ (i) as functions of model parameters.

(B14)
Appendix C: Validity of the DCM analysis The DCM solution involves replacing the numerator in Eq. (A7) with a normallydistributed random variable.The Berry-Esseen theorem [62,63] supplies an upper bound for the error in the cumulative-density function approximated in this way.Applying this to the case in hand, with κ (i) = 0 initially, we find that error bound to be where C is some constant.Since the first factor only contains moments of the fractional variables n The key point is that in both cases this value is independent of S. The result is that E (i) vanishes in the limit thanks to the explicit S −1/2 dependence of Eq. (C1).It seems plausible that this behaviour extends to any similar distribution with the same scaling.For κ (i) = 0, the fluctuations in λ also need to be taken into account.Considering these independently of the interactions and taking the distribution of K (i) to a gamma distribution, the error associated with the normal approximation will be proportional to Thus, with matrix elements as in Eq. ( 3), the error associated with the normality approximation scales like S −1/2 .This is consistent with other approximations in the derivation and means that in the S → ∞ limit, the DCM solution becomes exact.At finite S, the above limitations do have a bearing on the parameters for which we can expect agreement between numerics and DCM expressions.With distributions chosen as in Sec.G, close agreement requires √ S 1 and √ S κ (i) .The main effect of this is the restriction of the value of κ (i) used in the plots here to κ (i) ≤ 5.
In the alternative scaling scheme discussed in Eq. ( 21) with a (i,i+1) 0,α > 0, we obtain b (i,i+1) 0,α 3 ∼ √ S and therefore the overall error scales as S 0 .In this case we can no longer rely on the asymptotic limit to make our distributions converge to the Gaussian limit.Rather, the accuracy of the DCM solutions is dependent on other model parameters.
In particular, for interaction coefficients distributed according to a Gamma distribution, we would require σ (i) ρ (i+1) /µ (i) 1. Lifting the a (i,i+1) 0,α > 0 restriction changes this requirement.For example, we might then choose the elements themselves to be normally distributed, in which case b (i,i+1) 0,α 3 = 0 and any potential convergence issues vanish (in either scaling approach).
response to these perturbations can be captured by matrices v (i,j) with elements v ∀i, α.To derive the cavity equations, we consider an initial community with equilibrium abundances n

, 3 ≈ 1 .
it scales like S 0 .This leaves the scaling behaviour of the error term dependent on the second factor.To calculate b go back to Eq. (3) and consider the original distribution of the matrix elements a 72, i.e. a constant.For a gamma distribution we find b