Linking complex microbial interactions and dysbiosis through a disordered Lotka–Volterra model

We thus collect all the parameters estimated from the data in a vector $\mathit{\pi }=(h,q_d,q_0,K)$. As we will better detail in the Methods, we develop a moment matching inference algorithm to infer the model parameters $\mathit{\theta }$, as depicted in Figure 1. The idea of the method is to introduce a cost function $\mathcal{C}(\mathit{\theta }|\mathit{\pi })$, representing a total relative error for some self-consistent equations. If the parameters $\mathit{\theta }$ are such that the right part of the self-consistent equation equals the left part, the problem is considered solved. Because the landscape associated with this cost function presents several minima, we perform multiple optimization procedures to collect an ensemble of possible solutions, from which we retain the top 30. First, we find that different solutions ${\mathit{\theta }}^{\ast }$ of the optimization problem provide ecological insights into the underlying microbiome populations.

As originally predicted in Altieri et al., 2021, among all the parameters that define Equation 1, the only ones relevant for reproducing the theoretical phase diagram are the amplitude of demographic noise and the heterogeneity of interactions. The mean interaction strength, provided it is sufficiently positive, does not play a significant role. This prediction is fully confirmed by the inference procedure applied to the two microbiome datasets, allowing us to identify a universal signature that distinguishes healthy from unhealthy states. Figure 3a shows, indeed, that inferred noise ($T$) and interaction heterogeneity strength ($\sigma$) for healthy and diseased microbiomes are clustered in the two-dimensional plane.

In particular, the SAD for the healthy cohort is robust among the different solutions of the inference procedure, as depicted by the superposition of the different curves in the inset of Figure 3a. On the other hand, SADs inferred from unhealthy patients have high sensitivity to different solutions. In particular, some of them display a mode for high-abundance species (light red lines in Figure 3a), a signature of dominant strain in the gut. Consistently, the distribution of the interactions $P({\alpha }_{i,j})$ generated through the inferred parameters μ and σ is different between healthy and diseased cohorts, giving a distinct pattern of interactions (see Figure 3b), a result that is compatible with that found by Bashan et al., 2016. Remarkably, we find that dysbiosis reduces the heterogeneity of interaction strengths, a result also observed when taking correlations as a proxy for interactions (Seppi et al., 2023).

We then assess how close the inferred σ and $\beta =1/T$ (a.k.a. inverse temperature in a statistical physics approach) are to the critical RSB line of the dgLV ($\mathcal{R}=0$), evaluated by keeping all the other parameters constant (see Methods). We find again that the replicon values $\mathcal{R}$ corresponding to each solution of our optimization protocol are significantly different for the two investigated microbiome phenotypes (see Figure 4a). In particular, diseased microbiomes are closer to marginal stability within the RS ansatz (Altieri et al., 2021; Mézard et al., 1987; de Almeida and Thouless, 1978). Furthermore, by investigating the shape of the SAD given by Equation 2, we can estimate the ratio between niche (represented by species interaction) and neutral (represented by birth/death and immigration) ecological forces, which can be captured by the quantity $\psi$ (Wu et al., 2021). It detects the emergence of peaks in the SAD as a hallmark of niche processes (see Appendix 2).

Inspired by field-theory arguments (see Methods and Appendix, Section S2), we can call the mass of the theory the $m$ parameter, as defined above in Equation 2. In classical and quantum field theory, the particle–particle interaction embedded in the quadratic term is typically referred to as a mass source. In our context, $m=1-\beta {\sigma }^2(q_d-q_0)$ captures quadratic fluctuations of species abundances, as also appearing in the expression of the leading eigenvalue of the stability matrix. When $m\to 0$, the analytical order parameters diverge and the system enters the unphysical regime of unbounded growth. As such, the mass term can be considered a complementary stability measure, capable of capturing the transition to the unbounded growth regime.

In the model, two kinds of effects compete to shape the community structure. On the one hand, we have niche effects, encoded in disordered interactions and thus tracked by the parameters μ, σ, and $K$. Their overall effect is selective and tends to concentrate the SAD around the typical abundance value. On the other hand, we have neutral effects encoded in the stochastic dynamics and immigration, governing the low-abundance regime of the SAD. When the demographic noise amplitude is stronger than immigration ($\nu <1$, as in our case), the SAD exhibits a low-abundance integrable divergence. In the opposite scenario, for $\nu >1$, there is no divergence, and the SAD is modal. Since interactions are random, the probability of observing an internal mode can be estimated as the fraction of SADs realizations having non-trivial solutions to the stationary point equation. Such a quantity, dubbed as the niche–neutral ratio, can be analytically evaluated:

(5)

${\displaystyle {\displaystyle {\displaystyle \psi =\frac{1}{2}\text{Erfc}\left(\frac{{\zeta }^{\ast }+\overline{\zeta }}{\sqrt{2}{\sigma }_{\zeta }}\right)+\frac{1}{2}\text{Erfc}\left(\frac{{\zeta }^{\ast }-\overline{\zeta }}{\sqrt{2}{\sigma }_{\zeta }}\right),}}}$

where ${\zeta }^{\ast }=\sqrt{\frac{4(1-\nu )m}{\beta }}$ and $\overline{\zeta }=K-\mu h$. When $\psi \approx 1$, niche and neutral forces give comparable contributions to the dynamics, as both low-abundance divergence and a finite abundance mode coexist in the SAD. Finally, if the typical abundance diverges, we enter the unbounded growth phase, which means that the mass $m$ and the niche–neutral ratio $\psi$ are not independent, as suggested by the analytical expression for $\psi$. For an exhaustive derivation of this result, see Appendix 2. With the obtained model parameters, we are able to evaluate $m$ and $\psi$ for healthy and diseased microbiomes. Also, in this case, healthy and diseased microbiomes are visibly clustered, as shown in Figure 4. Unhealthy microbiomes turn out to be closer to the unbounded growth phase, and the niche–neutral ratio is larger by five orders of magnitude than the healthy case ${\psi }_D\approx {10}^5{\psi }_H$. This leads us to argue that selective pressure is way larger in diseased states, while in the healthy one, birth and death effects are the key drivers of the dynamics. These results are also confirmed by the SAD shapes in the inset of Figure 3 (panel a).

In summary, in the Results section, we show that (i) the inference pipeline robustly recovers demographic noise and interaction heterogeneity by calculating $h$, $q_0$, and $q_d$ from the data; and (ii) these parameters cluster according to health status, with diseased microbiomes lying closer to the replica-symmetry-breaking threshold, indicating reduced ecological resilience.