Incorporating Heterogeneous Interactions for Ecological Biodiversity

Introduction.— The introduction of the field-theoretic approach Martin et al. (1973); Janssen (1976); De Dominicis and Peliti (1978); Täuber (2014) for stochastic processes has led to its application across various dynamical systems, including neural networks Crisanti and Sompolinsky (2018); Sompolinsky et al. (1988); Helias and Dahem (2022), statistical learning Mannelli et al. (2020); Poole et al. (2016); Segadlo et al. (2022), and game theory Galla and Zhang (2009); Baron (2021); Opper and Diederich (1992). Dynamical mean-field theory (DMFT), initially devised for investigating spin-glass dynamics Sompolinsky and Zippelius (1981); Crisanti et al. (1993), has found success in ecological modeling, particularly in solving generalized random Lotka-Volterra (GRLV) equations of species abundance Tikhonov and Monasson (2017); Yang and Park (2023); Galla (2018); Bunin (2017). This breakthrough sheds light on how complexity impedes the stability and diversity of large ecological communities Galla (2018); Bunin (2017); Roy et al. (2019); Biroli et al. (2018), aligning with earlier findings from the random matrix theory May (1972); Allesina and Tang (2012); Baron et al. (2023).

While the current DMFT method is predominantly constrained to fully-interacting or well-mixed communities, real-world data reveal complicated and heterogeneous structures in ecological networks Jordano et al. (2003); Krause et al. (2003); Dunne et al. (2002); Guimarães Jr. (2020); Montoya and Solé (2002), making them archetypical examples of complex networks in network science. From the perspective of network science, understanding the impact of structural heterogeneity on dynamics stands as a central issue, given its universality across disciplines Newman (2018); Albert and Barabási (2002); Barrat et al. (2008). In response, theoretical tools such as the heterogeneous mean-field theory (HMFT) have been developed Dorogovtsev et al. (2008); Pastor-Satorras and Vespignani (2001). The HMFT, based on the expectation that agents with the same number of interacting partners exhibit identical dynamic behaviors, has successfully explained the critical phenomena in spin systems of a general network structure, offering insights distinct from the conventional mean-field results Kim et al. (2001); Dorogovtsev et al. (2002); Iglói and Turban (2002); Kim et al. (2005); Lee et al. (2009). Moreover, the HMFT yielded fruitful insights into epidemics Castellano and Pastor-Satorras (2010); Costa et al. (2022) and synchronization Lee (2005); Rodrigues et al. (2016) across various networks.

In Ref. Barbier et al. (2018), the DMFT has been tested for its performance in explaining ecological systems on complex networks by extracting the interaction mean and variance at the system level. Although this approach effectively describes ecological systems with relatively homogeneous interaction structures, it fails to provide accurate approximations for heterogeneous cases. As an endeavor to incorporate structural heterogeneity into the theoretical framework for the dynamics of ecological systems, solvable models beyond the well-mixed structures have been investigated Lee et al. (2022); Marcus et al. (2022); Poley et al. (2023). However, developing a general framework that can account for both structural heterogeneity and randomness in interaction strength remains a fundamental challenge. Such a framework is essential for comprehensively understanding the diversity and stability of ecological systems, but it has yet to be done.

In this Letter, we propose a generic theoretical framework, heterogeneous dynamical mean-field theory (HDMFT), combining the two theoretical methods, the DMFT from the field-theoretic approach and the HMFT from network science, and apply it to the GRLV dynamics on general ecological networks and obtain the species abundance distribution and survival probability at the system and individual level. The most remarkable finding is that the number of surviving species diminishes as the community becomes cooperative in heterogeneous interaction structures, which is counterintuitive and never seen in well-mixed systems. The origin lies in the differentiated abundance and survival of individual species with their numbers of interacting partners, and our detailed analysis reveals the interplay between the heterogeneous structure and random strength of interaction. This combined framework deepens the understanding of heterogeneous ecological systems and can be utilized in various interacting systems beyond the scope of natural ecosystems.

Model.— We first construct a GRLV system with $S$ species as follows:

where $x_{i}$ is the abundance of $i-$th species with a self-growth rate $\lambda$, and $\sum_{j\backslash i}$ denotes the summation over all species except $i$. The last term represents the interspecific interaction with two factors: (i) the strength of interaction $J_{ij}$ depicting how much a species $j$ affects another species $i$ positively or negatively and (ii) the structure of interaction (the adjacency matrix in the context of networks) $A_{ij}=1$ or $0$ indicating the existence or absence of interaction (see Fig. 1 for a graphical illustration). A positive (negative) value of $J_{ij}$ indicates the suppression (promotion) of the growth of $i$ by $j$. We set the adjacency matrix to be symmetric, i.e., $A_{ij}=A_{ji}$ while we consider $J_{ij}$ and $J_{ji}$ as independent or correlated to some extent controlled by a parameter.

These interspecific interactions, $J_{ij}$ and $A_{ij}$, do not need to be uniform in nature, and thus we assume the randomness in both $J_{ij}$ and $A_{ij}$. As in previous studies Bunin (2017); Galla (2018); Roy et al. (2019), for analytic treatment, we set the interaction strength $J_{ij}$ as a Gaussian random variable with the first few moments given by

where $K=S^{-1}\sum_{i,j}A_{ij}$. Notations ${\langle\cdots\rangle}_{\mathbf{J}}$ and ${\llangle\cdots\rrangle}_{\mathbf{J}}$ represent moment and cumulant, respectively, averaged over an ensemble of interaction strength realizations $\{\mathbf{J}_{1},\mathbf{J}_{2},\cdots\}$. The sign of $\mu$ represents whether the considered community is overall competitive $(\mu>0)$ or cooperative $(\mu<0)$. The $\sigma$ characterizes the randomness of interaction strength. The reciprocity $-1\leq r\leq 1$ is related to the type of pairwise interactions, e.g., $r=-1$ with predator-prey interactions $(J_{ij}=-J_{ji})$.

To represent a heterogeneous interaction structure, we consider an ensemble of adjacency matrices $\{\mathbf{A}_{1},\mathbf{A}_{2},\cdots\}$ with each element $A_{ij}$ taking $0$ and $1$ with probability $1-p_{ij}$ and $p_{ij}$ respectively, and thus satisfying

The number of interacting partners of species $i$ is given by $k_{i}=\sum_{j}A_{ij}$ called degree and its ensemble average $\langle k_{i}\rangle=\sum_{j\backslash i}p_{ij}$ can be heterogeneous if the connection probabilities $p_{ij}$ are not identical. We use $p_{ij}$ proposed in the static model Goh et al. (2001) to generate networks with the power-law degree distributions $P(k)\sim k^{-\gamma}$, where $\gamma$ is the degree exponent, or the Poisson distribution $P(k)=\mathrm{Pois}(k;K)=K^{k}e^{-K}/{k!}$ Erdős and Rényi (1959); Gilbert (1959). We here approximate the connection probability by $p_{ij}\approx k_{i}k_{j}/(SK)$, which is often called the annealed approximation and valid in uncorrelated networks without a significant degree-degree correlation Newman (2002); Dorogovtsev et al. (2008); Park and Newman (2004). Under this approximation, the statistical property of $A_{ij}$ is solely determined by the degree sequence $\{k_{i}\}$, leading to the dynamical equations for the HMFT.

HDMFT.— Using the DMFT technique with the HMFT for Eq. (1) and setting $r=0$ for simplicity, one can obtain the abundance dynamics for a species with degree $k$ as follows:

The outcome can be roughly understood by approximating $\sum_{j\backslash i}J_{ij}A_{ij}x_{j}(t)$ in Eq. (1) as the sum of $k_{i}$ independent and identically distributed random variables, resulting in $\mu\frac{k}{K}m(t)$ along with the Gaussian noise $\eta(t)$, where $\langle\eta(t)\rangle=0$ and $\langle\eta(t)\eta(t^{\prime})\rangle=q(t,t^{\prime})$. Here, $m(t)$ and $q(t,t^{\prime})$ are

Equations (3) and (4) are the main results from the HDMFT and their rigorous derivation with moment-generating functional for general $r$ is given in Supplemental Material (SM) Sec. I.A. sup . Note that $\langle\cdots\rangle$ is the average over both species and ensembles of $(\mathbf{J},\mathbf{A})$.

The equilibrium abundance $x(k)\equiv\lim_{t\to\infty}x(t;k)$ follows a truncated Gaussian distribution $\rho(x|k)$,

where $z$ is a random variable sampled from the standard Gaussian distribution $\mathcal{N}(0,1)$. Compared with the well-mixed system Bunin (2017); Galla (2018), differences are found in the terms proportional to $k$ and $\sqrt{k}$ in Eq. (5), which bring the degree dependence and fundamental changes to the abundance distribution and survival probability.

If the degrees of individual species follow a degree distribution $P(k)$, the stationary abundance distribution of the whole community is given by $\rho(x)=\sum_{k}P(k)\rho(x|k)$. When all the species have the same degree $K$, i.e., $P(k)=\delta_{k,K}$, the truncated Gaussian distribution is recovered, $\rho(x)=\rho(x|K)$ as in fully-interacting systems Bunin (2017); Galla (2018). However, when species exhibit varying degrees, the coalition no longer guarantees a truncated Gaussian distribution for $\rho(x)$. In Fig. 2, we present $\rho(x)$ obtained with the Poisson distribution $\mathrm{Pois}(k;K)$ and power-law distributions $P(k)\sim k^{-\gamma}$ with $\gamma=3.5$ and $2.5$, ordered by distribution width (the second cumulant). The DMFT approach extracts only the mean and variance of the overall interactions, including zero components, and approximates $J_{ij}A_{ij}$ as a Gaussian random variable with those mean and variance Barbier et al. (2018). While this method performs well in predicting abundance distributions when degree heterogeneity is relatively small, it fails in cases of large degree heterogeneity. In contrast, our HMDFT method agrees well with simulation results, even with significant degree heterogeneity. Furthermore, our approach accurately describes $\rho(x|k)$, distinguishing between species with varying numbers of interacting species, a capability inaccessible to the DMFT method alone (See Sec. II.C. in SM sup ).

Average abundance and survival probability.—The system-level influence of structural heterogeneity is best manifested in the average abundance $\langle x\rangle$ and the survival probability $\phi=\langle\Theta(x)\rangle$ with the Heaviside step function $\Theta(x)$. As shown in Fig. 3, for a competitive community $(\mu>0)$, both $\langle x\rangle$ and $\phi$ decrease monotonically with $\mu$, i.e., competition curbs the community from growing. On the other hand, in the cooperative community $(\mu<0)$, as $|\mu|$ gets larger, $\langle x\rangle$ increases but $\phi$ decreases. It is counterintuitive, as this result indicates that the more species benefit from each other, the more they eventually vanish, whereas the community itself proliferates. This dramatic drop in diversity reflected by $\phi$ is absent in well-mixed systems such as fully-connected $(K=S-1)$ case or relatively homogeneous interactions $[P(k)=\mathrm{Pois}(k;K)]$, indicating the important role of interaction structures. Without the explicit consideration of degree heterogeneity, the DMFT alone can never predict this phenomenon.

To excavate the reason behind this seemingly counterintuitive behavior appearing in the cooperative community, let us consider the case with no randomness in strength, $\sigma=0$. Then all species survive with nonzero abundance $x(k)=\lambda+|\mu|mk/K$, trivially leading to $\phi=1$. Thus, cooperative communities are feasible with homogeneous interaction strength by prohibiting survival probability from $\phi$ falling off. Next, let us consider the case without structural heterogeneity, i.e., $k_{i}=K$ for all $i$. Then one can find that both $m$ and $q$ increase with increasing $|\mu|$ but their ratio $q/m^{2}$ remains constant and so does the survival probability $\phi$, independent of $\mu$ (See Sec. I.C. in SM sup ). For small $\sigma$, we found that $\gamma<4$ is necessary to observe such a diversity drop (See Sec. I.G. in SM sup ). Therefore, without either heterogeneity in strength or structure is there no diversity drop in $\phi$ for $\mu<0$.

The counterintuitive drop of diversity with cooperation can be intuitively understood in a star-like network, where numerous “peripheral” species interact with a small number of “hub” species. The influences of the peripheral species on the hub, when summed up, are likely to be cooperative on average with a small fluctuation due to the central limit theorem. In contrast, from the viewpoint of peripherical species, the effect from the hub fluctuates under $J_{ij}\sim\mathcal{N}(\mu/K,\sigma^{2}/K)$. Therefore, the hub is very likely to benefit from the peripheral species while hubs could help or suppress the growth of peripheral species by fluctuations. Moreover, the expectedly large abundance of the hub can threaten a peripheral species even if only a weak competitive interaction is present between them. This is highly contrasted to the case of no significant heterogeneity in degrees, in which all species, topologically similar, essentially impose an averaged effect on one another with similar abundances.

Abundance and survival of individual species.— Given such difference between hub and peripheral species, we investigate how the fate of individual species is differentiated by their degrees. Using the HDMFT solution in Eq. (5), we calculate the survival probability of a species with degree $k$, represented by a cumulative Gaussian distribution $\phi(k)=\int_{-\infty}^{\infty}dx\>\Theta(x)\rho(x|k)=(2\pi)^{-1/2}\int_{-\infty}^{\Delta(k)}dz\exp(-z^{2}/2)$ with

The two terms in Eq. (6) represent the self-growth and the expected influences of $k$ neighboring species, respectively, as seen in rewriting Eq. (5) as $x(k)=\sigma\sqrt{qk/K}\max(0,\Delta(k)-z)$. Each term dominates for $k\ll k_{*}$ and $k\gg k_{*}$, respectively, with the characteristic degree $k_{*}\equiv\lambda K/(|\mu|m)$.

For the competitive case with $\mu>0$, species with large $k$ may have negative $\Delta(k)$ and thereby have small $\phi(k)$. Specifically, the survival probability $\phi(k)$ sharply decreases in the vicinity of $k_{*}$ [see Fig. 4(a)]. Furthermore, a strong competition drags $k_{*}$ down, manifested as decreasing $\phi$ with respect to $\mu$ in Fig. 3(a). On the other hand, for the cooperative case with $\mu<0$, $\Delta(k)$ becomes positive for all $k$, resulting in $\phi(k)>1/2$. For the species with as small degree as $k\ll k_{*}$ or as large as $k\gg k_{*}$, the self-growth or the expected sum of the partners’ influences is likely to dominate the fluctuation of the latter, resulting in large $\Delta(k)$ and $\phi(k)$. If a species’ degree $k$ is comparable to $k_{*}$, the fluctuation can be so significant as to make $x(k)$ small or zero, reducing $\Delta(k)$ and $\phi(k)$. We also find that due to $k_{*}$ decreasing with $\mu$, low-degree species are more likely to be eliminated under stronger cooperation [see Fig. 4(b)].

Phase diagram.— For a comprehensive comparison with previous studies Bunin (2017); Galla (2018), we investigate how the degree heterogeneity affects the phase diagram consisting of the unique fixed point (UFP), unbounded growth (UG), and multiple attractors (MA) phases in the $(\mu,\sigma)-$plane by using both numerical solutions to Eq. (1) and the HDMFT prediction. The results in Fig. 5 first reveal that the UFP phase shrinks as the degree heterogeneity increases. The triple point is given by $(\mu_{t},\sigma_{t})=(0,\sqrt{2K^{2}/\langle k^{2}\rangle})$ from the HDMFT approach, which converges to $(0,0)$ for $P(k)\sim k^{-\gamma}$ with $2<\gamma<3$ due to the divergence of the second moment $\langle k^{2}\rangle$ (See Sec. I.F. in SM sup ). Consequently, the system cannot transit directly from the UFP to the UG phase under such strong heterogeneity. Note that the triple point in Fig. 5(c) does not precisely converge to $(0,0)$ due to finite-size effects.

Summary and discussion.— We have proposed a novel MF framework capable of addressing two types of heterogeneity in the random Lotka-Volterra systems: the interaction strength with $\mathbf{J}$ and the interaction structure with $\mathbf{A}$. The obtained MF dynamics reveals the differentiation of the individual species’ behaviors depending on their degrees, helping us to understand the presence of the characteristic degree at which the randomness of strength is the most dominant and the origin of nontrivial emergent features of the whole community when both types of heterogeneity are present.

As a final cautionary remark, our approach is based on the moment-generating functional expanded up to the leading order of $K^{-1}$. The higher-order terms such as $\sim K^{-2}$ and quartic interactions $\sim x^{4}$ should be considered to study the case of $K$ not large enough, the investigation of which will clarify the effects of interaction sparsity. Furthermore, it will allow us to investigate whether our results that degree heterogeneity is a prerequisite for causing diversity drops in cooperative communities still hold for even smaller $K$.

Our framework can be extended in diverse directions, including the incorporation of general distributions of interaction strength, beyond the Gaussian distribution. Another challenging work can be understanding the topological properties of the community of surviving species, including the size and number distributions of the possibly fragmented subcommunities. Finally, we would like to emphasize the broad applicability of our method beyond natural ecosystems, as long as a system is connected by networks and describable through simple equations like GRLV. This includes diverse domains such as economic or political ecosystems, demonstrating the relevance of our theoretical framework across various fields.