Generalized Lotka-Volterra model with hierarchical interactions

I. INTRODUCTION

The study of complex ecosystems has been an active research area in theoretical ecology since the 1970s, when evidence emerged of a sharp transition from stability to instability in model ecosystems of increasing complexity . By suggesting fundamental limits on the size, connectedness, and interaction variability of a stable ecosystem, these results appeared to contradict the prevailing ecological view of the time. Many ecological networks are both large and highly interconnected , leading to an a priori expectation that a more complex and well-connected ecosystem ought to be more stable than its simpler, more sparsely connected counterpart . Such an apparent contradiction has led to an increasingly detailed and nuanced search for consistent definitions of ecological stability and complexity across theory and experiment .

One powerful theoretical tool for analyzing which factors contribute to the stability of a large system of many interacting constituents, such as a complex ecosystem, is random matrix theory (RMT). This is the method that was employed by May in his seminal work . May started from the Jacobian of a hypothetical ecosystem, linearized about an “equilibrium.” He assumed that the entries of this Jacobian matrix were independent and identically distributed random variables. This allowed for the deduction of a stability criterion using a result from RMT: Girko's circular law for the distribution of eigenvalues of large independent and identically distributed random matrices . May's initial and somewhat austere model has been extended in recent years, accounting for additional features of ecosystems such as food web structure , spatial dispersal , alternative interpretations of “interaction strength” , and varying self-regulation .

Despite providing a great deal of insight, the RMT approaches to determining the stability of complex ecological communities suffer from one crucial drawback: There is no guarantee that the random matrix under consideration corresponds to the Jacobian matrix of any real dynamics linearized about a feasible equilibrium , that is, an equilibrium at which all species abundances are non-negative. To remedy this, recent works have instead studied the stability of complex ecosystems using the generalized Lotka-Volterra equations, which produce feasible equilibria by construction. In such models, the interactions between an initial pool of species are chosen randomly, but certain species are allowed to go extinct as the system evolves towards a stable equilibrium. The stability of the surviving community is dependent on which species survive and their equilibrium abundances and cannot be determined from the eigenvalues of the initial random interaction matrix. Therefore, RMT alone cannot determine the stability of a surviving feasible community. We must use alternative methods to find the properties of the surviving community and subsequently its stability.

In this work we extend a recent study , which used RMT to study the stability of a linear model (in the style of May) with hierarchical interactions (referred to in as the cascade model ). Here we incorporate such hierarchical interactions in the generalized Lotka-Volterra equations. We obtain analytical results for both stability and the composition of surviving communities using techniques from statistical physics and the theory of disordered systems, specifically dynamic mean-field theory . Our approach has several advantages. Foremost, as mentioned above, the equilibria we study are feasible by construction. Further, we are able to study the effects of hierarchical interactions not only on stability, but also on the emergent properties of the community, such as species abundances and survival probabilities. Finally, we are also able to identify not only when the ecosystem becomes unstable, but also the nature of the instability.

More specifically, we show that increasing the severity of hierarchical interaction in the community is a stabilizing force, as is increasing the proportion of interactions which are of predator-prey type (in agreement with, for example, ). We also demonstrate that larger asymmetry in the variance of species' interactions decreases stability. Further, we look at the properties of stable equilibria produced in such a hierarchical community. We find that the presence of hierarchical interactions leads to communities with non-Gaussian species abundance distributions and that species lower in the hierarchy (i.e., species that benefit less from interactions) are both less likely to survive and less abundant.

This remainder of the paper is structured as follows. In Sec. we describe the generalized Lotka-Volterra model with hierarchical interactions. In Sec. we outline how dynamic mean-field theory is used to trade the set of coupled differential equations with random interaction coefficients constituting the generalized Lotka-Volterra model for a smaller number of statistically equivalent coupled stochastic differential equations. We then analyze the fixed-point solution found in Sec. , noting that the fixed-point equations were also obtained with the cavity method in Refs. [Eq. ], but no stability analysis was performed in these references. We then discuss the effect of hierarchical interactions on the distribution of species abundances and the fraction of surviving species in Sec. . In Sec. we study the stability of the fixed-point solution. We summarize in Sec. .

II. MODEL DEFINITION

A. General block-structured interactions

Consider a community of $𝑁$ species, partitioned into $𝐵$ distinct subcommunities, indexed with $𝑎 = 1, \dots, 𝐵$ . Subcommunity $𝑎$ contains $𝑁 𝑎$ species, so $𝑁 = \sum 𝑎 𝑁 𝑎$ . We first describe a more general model before specifically introducing hierarchical interactions.

We denote the abundance of the $𝑖 t h$ species in group $𝑎$ by $𝑥 𝑎 𝑖 (𝑡)$ , so $𝑖 = 1, \dots, 𝑁 𝑎$ when we are referring to species in group $𝑎$ . Species abundances evolve according to generalized Lotka-Volterra (GLV) dynamics

̇ 𝑥 𝑎 𝑖 (𝑡) = 𝑥 𝑎 𝑖 (𝑡) [1 - 𝑥 𝑎 𝑖 (𝑡) + 𝐵 \sum 𝑏 = 1 𝑁 𝑏 \sum 𝑗 = 1 𝐴 𝑎 𝑏 𝑖 𝑗 𝑥 𝑏 𝑗 (𝑡)],

(1)

where the block-structured interaction matrix $(𝐴 𝑎 𝑏) 𝑖 𝑗$ dictates the influence of the $𝑗 t h$ species in group $𝑏$ on the $𝑖 t h$ species in group $𝑎$ . The coefficients $𝐴 𝑎 𝑏 𝑖 𝑗$ are correlated random variables with statistics to be specified (which may vary between between blocks). The diagonal elements $𝐴 𝑎 𝑎 𝑖 𝑖$ are set to zero. This general setup is illustrated in Fig. .

FIG. 1.

Illustration of the block-structured interaction matrix $𝐴 𝑎 𝑏 𝑖 𝑗$ . The figure shows the mean value (indicated by color) of interaction matrix elements $𝐴 𝑎 𝑏 𝑖 𝑗$ in three cases. (a) General block-structured interaction matrix with $𝐵 = 4$ subcommunities. (b) Example of an interaction matrix with hierarchical statistics, again with $𝐵 = 4$ . All blocks above the diagonal blocks have the same mean and variance, and similarly for blocks below the diagonal. (c) Model with an infinite number of subcommunities ( $𝐵 \to \infty$ ) and hierarchical statistics [see Eq. ].

The statistics of the interactions between species depend only on their respective subcommunities. We write

𝐴 𝑎 𝑏 𝑖 𝑗 = 𝜇 𝑎 𝑏 𝑁 + 𝜎 𝑎 𝑏 \sqrt 𝑁 𝑤 𝑎 𝑏 𝑖 𝑗,

(2)

where the $𝑤 𝑎 𝑏 𝑖 𝑗$ are random variables with the first and second moments

‾ 𝑤 𝑎 𝑏 𝑖 𝑗 = 0, ‾ (𝑤 𝑎 𝑏 𝑖 𝑗) 2 = 1, ‾ 𝑤 𝑎 𝑏 𝑖 𝑗 𝑤 𝑏 𝑎 𝑗 𝑖 = 𝛾 𝑎 𝑏 .

(3)

We have indicated the average over realizations of the matrix $𝑤 𝑎 𝑏 𝑖 𝑗$ with an overbar.

The model is thus fully specified by the parameters $𝜇 𝑎 𝑏, 𝜎 𝑎 𝑏$ , and $𝛾 𝑎 𝑏$ . More specifically, $𝜇 𝑎 𝑏$ is the average influence of species in group $𝑏$ on those in group $𝑎, (𝜎 𝑎 𝑏) 2$ is the variance of those interactions, and $𝛾 𝑎 𝑏 \in [- 1, 1]$ is a correlation coefficient controlling the proportion of interactions between species in subcommunities $𝑎$ and $𝑏$ which are of a predator-prey type (i.e., $𝐴 𝑎 𝑏 𝑖 𝑗 𝐴 𝑏 𝑎 𝑗 𝑖 < 0$ ). The factors of $1 / 𝑁$ in Eq. ensure a sensible large- $𝑁$ limit (see, e.g., ).

B. Hierarchical interactions

Hierarchical interactions are obtained through a specific choice of the model parameters $𝜇 𝑎 𝑏, 𝜎 𝑎 𝑏$ , and $𝛾 𝑎 𝑏$ in Eqs. and . We imagine a ranking of the $𝐵$ subcommunities in which, on average, species with higher rank gain more from lower ranked species than vice versa. Hence species higher up in the hierarchy are increasingly biased towards success.

Similar to Ref. , we achieve this with the choice

𝜇 𝑎 𝑏 = {𝜇 - 𝜈, 𝑎 < 𝑏 𝜇, 𝑎 = 𝑏 𝜇 + 𝜈, 𝑎 > 𝑏, 𝜎 𝑎 𝑏 = {𝜎 / 𝜌, 𝑎 < 𝑏 𝜎, 𝑎 = 𝑏 𝜎 𝜌, 𝑎 > 𝑏, 𝛾 𝑎 𝑏 = 𝛾,

(4)

where we mostly focus on the case $𝜈 > 0$ . In this case species $𝑎 = 1$ is lowest in the hierarchy and species $𝑎 = 𝐵$ is highest. If we make the transformation $𝜈 \to - 𝜈$ and $𝜌 \to 1 / 𝜌$ , the ranking of species is reversed and we obtain identical results. We also require $𝜌 > 0$ , ensuring the variance of all interactions remains positive. An example of the structure of the interaction matrix resulting from this choice is illustrated in Fig. .

The model parameters $𝜇, 𝜈, 𝜎, 𝜌$ , and $𝛾$ can be separated into two sets: $𝜈$ and $𝜌$ describe the hierarchy of species, whereas $𝜇, 𝜎$ , and $𝛾$ describe the overall statistics. The parameter $𝜈$ is a measure of the strength of the hierarchy. A larger value of $𝜈$ increases the average benefit to species higher in the hierarchy and correspondingly decreases the average benefit to species lower in the hierarchy. The parameter $𝜌$ is a measure of the disparity between the variances of interactions $(𝜎 𝑎 𝑏) 2$ and $(𝜎 𝑏 𝑎) 2$ for $𝑎 \neq 𝑏$ . The parameter $𝜇$ characterizes the mean interaction strength across all pairs of species in the community, $𝜎$ is a measure of the overall variability of interactions, and $𝛾$ can be related to the proportion of interactions that are of predator-prey type $𝑝$ via $𝛾 = c o s (𝜋 𝑝) + 𝑂 (1 / \sqrt 𝑁)$ [see Sec. S2 of the Supplemental Material for details]. If $𝜈 = 0$ and $𝜌 = 1$ , there is no distinction between any two positions in the hierarchy and all species are statistically equivalent; this is the case considered in .

It is important to note that we distinguish between the notions of hierarchy, controlled by $𝜈$ , and the proportion of predator-prey pairs in the model, controlled by $𝛾$ . Hierarchy is present so long as $𝜈 \neq 0$ , in which case there is a natural ordering to the subgroups of species based on the relative average benefit of species higher up in the hierarchy. The proportion of predator-prey pairs, on the other hand, is not affected by hierarchy (i.e., a nonzero value of $𝜈$ ) to leading order in $𝑁$ (see Sec. S2 in ). This is due to the $1 / 𝑁$ scaling of the average value of $𝐴 𝑎 𝑏 𝑖 𝑗$ in Eq. compared to the $1 / \sqrt 𝑁$ scaling of the random term. Species higher in the hierarchy are therefore no more likely to predate on species lower down than species lower down in the hierarchy are to predate on species higher up. However, species higher in the hierarchy have higher average interaction coefficients than those lower in the hierarchy.

III. DYNAMIC MEAN-FIELD THEORY

We use dynamic mean-field theory in order to analyze the stability and community properties of the GLV system in Eq. . Ultimately our analysis will focus on models with hierarchical interactions, determined by the parameters $𝜇, 𝜈, 𝜎, 𝜌$ , and $𝛾$ . However, for now, we return to the more general block-structured general case discussed in Sec. [Fig. ].

The analysis involves taking the limit $𝑁 \to \infty$ while keeping the ratios $𝑛 𝑎 \equiv 𝑁 𝑎 / 𝑁$ constant such that

\sum 𝑎 𝑛 𝑎 = 1 .

(5)

That is to say, each subcommunity contains a large (formally infinite) number of species, but the proportions of species in each group remain fixed as $𝑁$ is varied.

The calculation closely follows the lines of , with modifications made to account for the block structure of the interactions in the community. Details are given in Sec. S3 in .

The dynamic mean-field approach results in the reduction of the initial set of coupled ordinary equations with random coefficients [given in Eq. and with $𝑁 \to \infty$ ] to a set of $𝐵$ stochastic integro-differential equations, one for each subcommunity. These describe the “typical” time evolution for the abundance of species in the different groups $𝑥 𝑎 (𝑡)$ . Carrying out the steps in Sec. S3 in , we arrive at the effective dynamics for species in subcommunity $𝑎$ ,

̇ 𝑥 𝑎 (𝑡) = 𝑥 𝑎 (𝑡) (1 - 𝑥 𝑎 (𝑡) + 𝐵 \sum 𝑏 = 1 𝑛 𝑏 𝜇 𝑎 𝑏 𝑀 𝑏 (𝑡) + 𝐵 \sum 𝑏 = 1 𝑛 𝑏 𝛾 𝑎 𝑏 𝜎 𝑎 𝑏 𝜎 𝑏 𝑎 \int 𝑡 0 𝑑 𝑡' 𝐺 𝑏 (𝑡, 𝑡') 𝑥 𝑎 (𝑡') + 𝜂 𝑎 (𝑡)) .

(6)

The variables ${𝜂 𝑎 (𝑡)}$ are colored Gaussian noise terms with the statistics

⟨ 𝜂 𝑎 (𝑡) ⟩ = 0, ⟨ 𝜂 𝑎 (𝑡) 𝜂 𝑏 (𝑡') ⟩ = 𝛿 𝑎 𝑏 𝐵 \sum 𝑐 = 1 (𝜎 𝑎 𝑐) 2 𝑛 𝑐 𝐶 𝑐 (𝑡, 𝑡'),

(7)

where we have used $⟨ \dots ⟩$ to denote averages over realizations of the effective dynamics in Eq. , i.e., over realizations of the noise ${𝜂 𝑎 (𝑡)}$ .

The macroscopic statistics of community $𝑎$ are given by

𝑀 𝑎 (𝑡) \equiv ⟨ 𝑥 𝑎 (𝑡) ⟩, 𝐶 𝑎 (𝑡, 𝑡') \equiv ⟨ 𝑥 𝑎 (𝑡) 𝑥 𝑎 (𝑡') ⟩, 𝐺 𝑎 (𝑡, 𝑡') \equiv ⟨ 𝛿 𝑥 𝑎 (𝑡) 𝛿 𝜂 𝑎 (𝑡') ⟩ .

(8)

The above quantities describe the average abundance of a species in subcommunity $𝑎$ , the autocorrelations (in time) of a species abundance in the community, and the response to perturbations, respectively. The solution of Eqs. determines the macroscopic statistics $𝑀 𝑎 (𝑡), 𝐶 𝑎 (𝑡, 𝑡')$ , and $𝐺 𝑎 (𝑡, 𝑡')$ self-consistently.

IV. FIXED-POINT EQUATIONS

A. General block-structured interactions

We now assume that the system reaches a fixed point such that $𝑥 𝑎 (𝑡) \to 𝑥 𝑎 *$ as $𝑡 \to \infty$ , where $𝑥 𝑎 *$ is a static random variable. We write

𝑀 𝑎 \equiv ⟨ 𝑥 𝑎 * ⟩, 𝑞 𝑎 \equiv ⟨ (𝑥 𝑎 *) 2 ⟩

(9)

for the first and second moments of the asymptotic abundances in each group of species $𝑎$ . The noise term $𝜂 𝑎 (𝑡) \to 𝜂 𝑎 *$ also asymptotically loses its time dependence and we write

𝜂 𝑎 * = 𝑧 𝑎 \sqrt \sum 𝑏 (𝜎 𝑎 𝑏) 2 𝑛 𝑏 𝑞 𝑏,

(10)

where the $𝑧 𝑎$ are independent zero-mean static Gaussian random variables with unit variance. In this fixed-point regime the response function $𝐺 𝑎 (𝑡, 𝑡') \to 𝐺 𝑎 (𝑡 - 𝑡')$ only depends on time differences $𝜏 = 𝑡 - 𝑡'$ and causality dictates that $𝐺 𝑎 (𝜏) = 0$ for $𝜏 < 0$ . We then define the integrated response function for subcommunity $𝑎$ ,

𝜒 𝑎 \equiv \int \infty 0 𝑑 𝜏 𝐺 𝑎 (𝜏) .

(11)

With Eqs. in mind, we find that the nontrivial fixed point of Eq. is described by

𝑥 𝑎 * = m a x ⎛ ⎜ ⎝ 0, 1 + \sum 𝑏 𝜇 𝑎 𝑏 𝑛 𝑏 𝑀 𝑏 + 𝑧 𝑎 \sqrt \sum 𝑏 (𝜎 𝑎 𝑏) 2 𝑛 𝑏 𝑞 𝑏 1 - \sum 𝑏 𝛾 𝑎 𝑏 𝜎 𝑎 𝑏 𝜎 𝑏 𝑎 𝑛 𝑏 𝜒 𝑏 ⎞ ⎟ ⎠ .

(12)

A similar expression was found for models without hierarchical structure in . We note that, depending on the value that the random variable $𝑧 𝑎$ takes, some species will be extinct at the fixed point ( $𝑥 𝑎 * = 0$ ), whereas others will have positive abundance.

The average over realizations of the effective process (denoted by angular brackets) is now an average over the static Gaussian random variables $𝑧 𝑎$ at the fixed point. This allows us to obtain self-consistent conditions for the statistics of $𝑥 𝑎 *$ . We find (see Sec. S4.1 in )

𝜒 𝑎 𝑢 𝑎 = 𝑤 0 (Δ 𝑎), 𝑀 𝑎 𝑢 𝑎 = 𝑤 1 (Δ 𝑎) \sqrt \sum 𝑏 (𝜎 𝑎 𝑏) 2 𝑛 𝑏 𝑞 𝑏, 𝑞 𝑎 (𝑢 𝑎) 2 = 𝑤 2 (Δ 𝑎) \sum 𝑏 (𝜎 𝑎 𝑏) 2 𝑛 𝑏 𝑞 𝑏,

(13)

where we have abbreviated

𝑢 𝑎 = 1 - \sum 𝑏 𝛾 𝑎 𝑏 𝜎 𝑎 𝑏 𝜎 𝑏 𝑎 𝑛 𝑏 𝜒 𝑏, Δ 𝑎 = 1 + \sum 𝑏 𝜇 𝑎 𝑏 𝑛 𝑏 𝑀 𝑏 \sqrt \sum 𝑏 (𝜎 𝑎 𝑏) 2 𝑛 𝑏 𝑞 𝑏

(14)

and where we have defined the functions

𝑤 𝑘 (Δ) = 1 \sqrt 2 𝜋 \int Δ - \infty 𝑑 𝑧 𝑒 - 𝑧 2 / 2 (𝑧 - Δ) 𝑘

(15)

for $𝑘 \in {0, 1, 2}$ . Equations are equivalent to those in Sec. IV.1 of the Supplementary Information of Ref. , derived with the cavity method. They also reduce to the fixed-point equations for the model without hierarchy found in Refs. if one takes $𝜇 𝑎 𝑏 = 𝜇 / 𝐵, 𝜎 𝑎 𝑏 = 𝜎 / \sqrt 𝐵$ , and $𝛾 𝑎 𝑏 = 𝛾$ .

Equations can be solved numerically for the quantities $𝜒 𝑎, 𝑀 𝑎$ , and $𝑞 𝑎$ as functions of the model parameters $𝜇 𝑎 𝑏, 𝜎 𝑎 𝑏, 𝛾 𝑎 𝑏$ , and $𝑛 𝑎$ . We note that the fraction of surviving species in subcommunity $𝑎$ , written $𝜙 𝑎$ , is given by

𝜙 𝑎 = 𝑤 0 (Δ 𝑎) .

(16)

Hence the solution of Eq. also provides the fraction of surviving species in each subcommunity.

For further analysis, it is useful to introduce the average over communities

〚 𝑋 〛 𝐵 = 𝐵 \sum 𝑎 = 1 𝑛 𝑎 𝑋 𝑎

(17)

for a quantity $𝑋 𝑎$ defined in each community. The overall fraction of surviving species in the system is then $〚 𝜙 〛 𝐵$ and the average abundance per species is $〚 𝑀 〛 𝐵$ .

B. Hierarchical interactions

We now analyze the fixed-point solution for a community with hierarchical interactions, for which the parameters $𝜇 𝑎 𝑏, 𝜎 𝑎 𝑏$ , and $𝛾 𝑎 𝑏$ are as in Eq. . It is convenient to introduce $𝛼 = 𝑎 / 𝐵$ and to work in the limit $𝐵 \to \infty$ ; the resulting interaction matrix is illustrated in Fig. . The index $𝛼 \in [0, 1)$ is now continuous and the constraint $\sum 𝐵 𝑏 = 1 𝑛 𝑏 = 1$ becomes $\int 10 𝑑 𝛼 𝑛 (𝛼) = 1$ (see also Sec. S4.1 in ). Quantities such as $𝑀 𝑎, 𝜙 𝑎$ , and $Δ 𝑎$ are now functions of $𝛼$ and averages over the index $𝑎$ are integrals, which we write as

〚 𝑋 〛 \equiv l i m 𝐵 \to \infty 〚 𝑋 〛 𝐵 = \int 10 𝑑 𝛼 𝑛 (𝛼) 𝑋 (𝛼),

(18)

dropping the subscript $𝐵$ in the limit. For example, the first relation in Eq. is now

𝑢 (𝛼) = 1 - 𝛾 𝜎 2 〚 𝜒 〛,

(19)

informing us that $𝑢$ has no dependence on $𝛼$ in the hierarchical model when $𝐵 \to \infty$ .

Recalling the definition of $Δ 𝑎$ in the second relation in Eq. , we also introduce the quantities

Δ 0 \equiv Δ (𝛼 = 0), Δ 1 \equiv Δ (𝛼 = 1) .

(20)

To proceed, we first find expressions for the quantities $𝑢, Δ 0$ , and $Δ 1$ , in terms of $𝜈, 𝜎, 𝜌$ , and $𝛾$ using Eq. ; these quantities are defined as a mathematical convenience. They will allow us to write explicit expressions for the macroscopic quantities of the system in which we are ultimately interested. Then we will find expressions for $〚 𝑀 〛$ and $〚 𝜙 〛$ [we also discuss expressions for the full functions $𝑀 (𝛼)$ and $𝜙 (𝛼)$ in Sec. ] in terms of $𝑢, Δ 0$ , and $Δ 1$ . In turn, this will enable us to calculate the macroscopic statistics of the system for any $𝜇, 𝜈, 𝜎, 𝜌$ , and $𝛾$ . Details of all remaining calculations in this section can be found in Sec. S4 in . Manipulation of Eq. reveals that $𝑢$ satisfies [see Eqs. (S49) and (S57) in Sec. S4 in ]

𝑢 2 = ℓ 𝜎 2 \int 10 𝑑 𝛼 𝑛 (𝛼) 𝑤 2 [Δ (𝛼)], 𝑢 (1 - 𝑢) = 𝛾 𝜎 2 \int 10 𝑑 𝛼 𝑛 (𝛼) 𝑤 0 [Δ (𝛼)],

(21)

where $ℓ$ is the logarithmic mean of $𝜌 2$ and $1 / 𝜌 2$ ,

ℓ \equiv 𝜌 2 - 1 / 𝜌 2 l n 𝜌 2 - l n 1 / 𝜌 2,

(22)

which satisfies $ℓ \geq 1$ , with equality only if $𝜌 = 1$ . We now show that the integrals in Eq. can be written explicitly in terms of only $𝑢, Δ 0$ , and $Δ 1$ . We achieve this by finding the separable differential equation for $Δ (𝛼)$ [Eq. (S67) in Sec. S4 in ],

𝑑 Δ 𝑑 𝛼 = 𝑛 (𝛼) 𝜃 (Δ, 𝑢),

(23)

with initial condition $Δ (0) = Δ 0$ and with

𝜃 (Δ, 𝑢) \equiv 𝑢 2 2 𝑢 𝜈 𝑤 1 (Δ) - ℓ 𝜎 2 l n 𝜌 2 Δ 𝑤 2 (Δ) .

(24)

Self-consistently, the solution to Eq. will further have to satisfy the constraint $Δ (1) = Δ 1$ .

One then uses Eq. to change variables in the integrals in Eq. , finding

\int 10 𝑑 𝛼 𝑛 (𝛼) 𝑤 2 [Δ (𝛼)] = \int Δ 1 Δ 0 𝑑 Δ 𝜃 (Δ, 𝑢) 𝑤 2 (Δ),

(25)

and similarly for the integral over $𝑤 0 [Δ (𝛼)]$ . If we also perform the same change of variables on the condition $\int 10 𝑑 𝛼 𝑛 (𝛼) = 1$ , then we arrive at the simultaneous equations for $𝑢, Δ 0$ , and $Δ 1$ ,

𝑢 2 = ℓ 𝜎 2 \int Δ 1 Δ 0 𝑑 Δ 𝜃 (Δ, 𝑢) 𝑤 2 (Δ), 𝑢 (1 - 𝑢) = 𝛾 𝜎 2 \int Δ 1 Δ 0 𝑑 Δ 𝜃 (Δ, 𝑢) 𝑤 0 (Δ), 1 = \int Δ 1 Δ 0 𝑑 Δ 𝜃 (Δ, 𝑢) .

(26)

Given the parameters $𝜈, 𝜎, 𝜌$ , and $𝛾$ , Eq. can be solved numerically to obtain $Δ 0, Δ 1$ , and $𝑢$ . As the parameter $𝜇$ and the function $𝑛 (𝛼)$ do not appear in Eq. , we conclude that $Δ 0, Δ 1$ , and $𝑢$ are independent of them. That is, $Δ 0, Δ 1$ , and $𝑢$ are functions of $𝜈, 𝜎, 𝜌$ , and $𝛾$ only.

Once $Δ 0, Δ 1$ , and $𝑢$ are determined, we calculate the average abundance and fraction of surviving species in the community with

1 〚 𝑀 〛 = Δ 1 𝜌 2 + Δ 0 Δ 1 𝜌 2 - Δ 0 𝜈 - 𝜇,

(27a)

〚 𝜙 〛 = 𝑢 (1 - 𝑢) 𝛾 𝜎 2 .

(27b)

Equation ceases to apply when

𝜈 = 0

, and a limit must be taken. Similar care has to be taken when finding

〚 𝜙 〛

in the limit

𝛾 \to 0

. Both limits are found in Sec. S4.5.3 in .

Inspection of Eqs. and reveals that, surprisingly, neither $〚 𝜙 〛$ nor $〚 𝑀 〛$ depends on $𝑛 (𝛼)$ , the relative sizes of the different subcommunities $𝛼$ . In fact, we will see in Secs. and that $𝑛 (𝛼)$ does not affect any of the properties of the community that we are interested in (see also Sec. S8 in ). Figure gives us insight into why this is the case. When the interactions are hierarchical, rather than a generally block structured as in Fig. , varying $𝑛 (𝛼)$ only affects the size of the diagonal blocks in the interaction matrix [Fig. ]. In the large- $𝐵$ limit, varying the sizes of the diagonal blocks (which are small compared to the size of the matrix) becomes insignificant [Fig. ]. We emphasize that the independence of our results from $𝑛 (𝛼)$ is a feature only of the present hierarchical model in the limit $𝐵 \to \infty$ and is not a general feature of the Lotka-Volterra model with block-structured interactions. One also sees that $〚 𝜙 〛$ is independent of $𝜇$ .

FIG. 2.

Variation of average abundance $〚 𝑀 〛$ and the fraction of surviving species $〚 𝜙 〛$ with $𝜈$ (the strength of the hierarchy). (a) Dependence of average abundance on $𝜈$ . The vertical dashed line indicates a divergence in average abundance when $𝛾 = 1$ and $𝜈 \approx 0.4$ . The system does not converge for $𝜈 ≲ 0.4$ . (b) Dependence of the fraction of surviving species on $𝜈$ . When abundances diverge (for $𝛾 = 1$ and $𝜈 \leq 0.4$ ), the fraction of surviving species cannot be meaningfully calculated and so the lower curve does not extend to the left of the vertical line. Markers are data from the numerical integration of Eq. , with $𝑁 = 𝐵 = 250$ , averaged over 40 runs, and lines are predictions from theory [Eqs. and , with $𝑢, Δ 0$ , and $Δ 1$ determined from Eq. ]. The model parameters used in both plots are $𝜎 = 0.3, 𝜌 = 1.5$ , and $𝜇 = 1.0$ .

Our theoretical predictions for the community average species abundance $〚 𝑀 〛$ and the community average survival fraction $〚 𝜙 〛$ are verified using the results of computer simulation of Eq. in Fig. . As an additional test of Eq. specifically, we also derive a condition under which only one species goes extinct in the nonhierarchical limit ( $𝜈 = 0$ and $𝜌 = 1$ ). A similar calculation was performed in Ref. , also in the absence of hierarchy. The calculation and the results can be found in Sec. S9 in .

V. FIXED-POINT DISTRIBUTIONS

A. Abundance distributions

At a stable equilibrium, we can calculate species abundance distributions (SADs) and rank abundance distributions (RADs). An SAD is obtained in simulations by binning species according to their abundances and producing a normalized histogram of the number of species in each bin [Fig. ]. An RAD is a plot of abundance against a ranking of species from 0 to 1, with the highest abundance species having a rank of 0 [Fig. ]. For reviews of both SADs and RADs see Refs. , for similar calculations without hierarchical interactions see Ref. , and for RADs derived from random replicator equations rather than the Lotka-Volterra equation see .

FIG. 3.

Abundance distributions with and without hierarchical interactions. All plots are for $𝜇 = 0.5, 𝜎 = 0.15$ , and $𝛾 = - 0.3$ . (a) Two SADs are plotted; the taller one is for the model without hierarchy ( $𝜈 = 0$ and $𝜌 = 1$ ) and is a clipped Gaussian. The flatter distribution in (a) is for $𝜈 = 1.5$ and $𝜌 = 23$ ; bars are from simulation. Also shown are the corresponding (b) RADs and (c) HADs. Square markers are from simulations for $𝜈 = 0$ and $𝜌 = 1$ and diamond markers are from simulations for $𝜈 = 1.5$ and $𝜌 = 23$ . Simulations are for $𝑁 = 𝐵 = 1000$ species, averaged over 200 runs. Lines are from the theory.

We also introduce hierarchical abundance distributions (HADs). Similar to an RAD, we rank species from 0 to 1, but this time such that a species with rank $𝑟$ is higher in the hierarchy than $𝑟 𝑁$ species and lower than $(1 - 𝑟) 𝑁$ species. An HAD is then a plot of abundance as a function of this hierarchical rank. The SADs and RADs are both derived from $Prob (𝑥 | 𝛼)$ , the probability that a species has abundance $𝑥$ , given that it sits at position $𝛼$ in the hierarchy. This is the probability density for $𝑥 𝑎 *$ in Eq. . One arrives at (see also Sec. S5 in ) the clipped Gaussian distribution

Prob (𝑥 | 𝛼) = [1 - 𝜙 (𝛼)] 𝛿 (𝑥) + 𝐻 (𝑥) 𝜑 (𝑥 | 𝛼),

(28)

where $𝐻 (\cdot)$ is the Heaviside step function [ $𝐻 (𝑥) = 1$ for $𝑥 \geq 0$ and $𝐻 (𝑥) = 0$ otherwise],and the delta function $𝛿 (𝑥)$ represents species that have gone extinct and hence have an abundance of $𝑥 = 0$ . The function $𝜑 (𝑥 | 𝛼)$ is a (normalized) Gaussian in $𝑥$ such that $\int \infty 0 𝑑 𝑥 𝜑 (𝑥 | 𝛼) = 𝜙 (𝛼)$ ,

𝜑 (𝑥 | 𝛼) = 𝑤 1 [Δ (𝛼)] 𝑀 (𝛼) \sqrt 2 𝜋 e x p [- 12 (𝑤 1 [Δ (𝛼)] 𝑀 (𝛼) 𝑥 - Δ (𝛼)) 2] .

(29)

The probability density function $Prob (𝑥 | 𝛼)$ is the SAD for species at position $𝛼$ in the hierarchy. We can therefore obtain the SAD for the whole community by integrating over all positions

〚 Prob (𝑥) 〛 \equiv \int 10 𝑑 𝛼 𝑛 (𝛼) Prob (𝑥 | 𝛼) .

(30)

To calculate $〚 Prob (𝑥) 〛$ explicitly, first Eq. are solved to obtain $Δ 0, Δ 1$ , and $𝑢$ . We then reexpress Eq. as a function of $Δ (𝛼)$ and use the same change of variables as in Eq. to calculate the integral , obtaining

〚 Prob (𝑥) 〛 = (1 - 〚 𝜙 〛) 𝛿 (𝑥) + 𝐻 (𝑥) \int Δ 1 Δ 0 𝑑 Δ 𝜃 (Δ, 𝑢) 𝜑 (𝑥 | Δ)

(31)

(details are found in Sec. S5 in ). Similarly to $〚 𝑀 〛$ and $〚 𝜙 〛$ , the distribution $〚 Prob (𝑥) 〛$ is independent of $𝑛 (𝛼)$ (see Sec. S8 in ).

In the case where there is no hierarchy ( $𝜈 = 0$ and $𝜌 = 1$ ), the underlying unclipped distribution $〚 𝜑 (𝑥) 〛 = \int 10 𝑑 𝛼 𝑛 (𝛼) 𝜑 (𝑥 | 𝛼)$ is itself a Gaussian distribution. However, as demonstrated in Fig. , this simple form is lost in a hierarchical community and $〚 𝜑 (𝑥) 〛$ is no longer Gaussian and is not even symmetric around its maximum. Broadly speaking, increasing the value of $𝜈$ lowers the modal abundance and increasing the value of $𝜌$ increases the spread in abundances.

Rank abundance distributions are also calculated from the distribution $〚 Prob (𝑥) 〛$ . We observe that if species are ranked on a scale from 0 to 1 by descending abundance, then the rank of a species with abundance $𝑥$ is $\int \infty 𝑥 〚 Prob (𝑥') 〛 𝑑 𝑥'$ . The plot in Fig. shows abundance on the vertical axis and rank on the horizontal axis. Rank abundance distributions are often the preferred representation of species abundances as they do not suffer from loss of information due to species binning, as SADs do .

To compute a hierarchical abundance distribution we rank species on a scale $𝑟 \in [0, 1]$ . Species with index $𝛼$ have rank

𝑟 (𝛼) = \int 𝛼 0 𝑑 𝛼 𝑛 (𝛼),

(32)

ensuring that $𝑟 𝑁$ species are lower in the hierarchy and $(1 - 𝑟) 𝑁$ species are higher in the hierarchy. One then produces an HAD, such as in Fig. , parametrically, with the horizontal axis equal to the rank and the vertical axis equal to the abundance $𝑀 (𝛼)$ , which is derived from Eq. in Sec. S5 in and given by

𝑀 (𝛼) = 𝐴 〚 𝑀 〛 𝑤 1 [Δ (𝛼)] \times e x p (ℓ 𝜎 2 l n 𝜌 𝑢 2 \int Δ (𝛼) Δ 0 𝑑 Δ' 𝜃 (Δ', 𝑢) 𝑤 2 (Δ')),

(33)

where the constant $𝐴$ ensures that $\int 10 𝑑 𝛼 𝑛 (𝛼) 𝑀 (𝛼) = 〚 𝑀 〛$ for $〚 𝑀 〛$ given by Eq. . Surprisingly, just like $〚 𝑀 〛, 〚 𝜙 〛$ , and $〚 Prob (𝑥) 〛$ , HADs are also independent of $𝑛 (𝛼)$ (this is shown in Sec. S8 in ).

B. Fraction of surviving species

A given species' survival probability as a function of the ranking $𝑟 (𝛼)$ can be produced in a similar way to the hierarchical abundance distributions discussed in Sec. . We simply replace $𝑀 (𝛼)$ by $𝜙 (𝛼) = 𝑤 0 [Δ (𝛼)]$ . By the same reasoning, we deduce that a plot of survival probability against rank $𝑟 (𝛼)$ is independent of both $𝜇$ and $𝑛 (𝛼)$ (see Sec. S8 in ). The distribution of survival probabilities is flat if $𝜈 = 0$ and $𝜌 = 1$ and is an increasing function of $𝑟$ in the presence of hierarchy (see Fig. ). This is an indication that the composition of a hierarchical community before the dynamics are run is vastly different from the resulting stable community; species lower in the hierarchy are both less abundant and less likely to survive. On further noting that the area under an HAD is $〚 𝑀 〛$ and the area under the survival curve is $〚 𝜙 〛$ , we see that introducing hierarchy produces smaller communities dominated by high ranked species. Further, from Figs. and , the mean abundance and the fraction of surviving species in communities with hierarchy are mostly lower than those in models without hierarchy (flat curve in both), demonstrating that few species benefit at all from hierarchical interactions.

FIG. 4.

Survival distributions with and without hierarchical interactions. The parameters used are $𝜎 = 0.75, 𝛾 = - 0.6$ , and, from top to bottom along the left edge of the plot, $(𝜈, 𝜌) = (0, 1), (0, 3), (3, 3), (3, 1), (3, 13)$ . Markers are the average fraction of survivors for every 40th species in simulations with $𝑁 = 𝐵 = 500$ , averaged over 400 runs.

VI. STABILITY

A. Stability conditions

In Sec. we found the statistics of the surviving species abundances by presuming a static solution to Eq. . In this section we discuss when this fixed-point solution is valid and thus under what conditions we have a stable and feasible equilibrium. As in Refs. , we find that instability can occur either through linear instability against small perturbations to the abundances or through a divergence in species abundances.

1. Linear instability

Following along the lines of (see Sec. S6 in ), we use a linear stability analysis to find that the system is unstable to perturbations in species abundances when

𝜎 2 \geq ℓ (ℓ + 𝛾) 2 1 〚 𝜙 〛,

(34)

where the average survival probability $〚 𝜙 〛$ is to be determined from Eq. and . In the limit of no hierarchy ( $𝜈 = 0$ and $𝜌 = 1$ ), Eq. reduces to the stability condition found in Refs. (Sec. S6.3.3 in ). Interestingly, the same criterion as in Eq. can be obtained with the machinery of random matrix theory, noting that $𝑁 * = 〚 𝜙 〛 𝑁$ species survive the dynamics asymptotically and reach a fixed point. The bulk of the eigenvalue spectrum of an $𝑁 * \times 𝑁 *$ random matrix with hierarchical statistics as in Eq. (cf. Fig. ) crosses the imaginary axis precisely when Eq. is satisfied . Similar observations were made in the case without hierarchy ( $𝜈 = 0$ and $𝜌 = 1$ ) in Ref. . A crucial difference between the dynamical and random matrix approaches is that the fraction of survivors $〚 𝜙 〛 = 𝑁 * / 𝑁$ is determined from Eq. in the dynamical theory, but $〚 𝜙 〛$ is an independent parameter of the model in the random matrix approach. In particular, Eq. would lead one to conclude that $𝜈$ has no effect on linear stability if a random matrix approach were used. However, as $〚 𝜙 〛$ is itself a function of the parameters $𝜈, 𝜎, 𝜌$ , and $𝛾$ , no such conclusion is drawn here.

2. Diverging abundances

To find the point at which species abundances diverge, we solve Eq. , together with Eq. , for the point at which $1 / 〚 𝑀 〛 = 0$ . We find that abundances diverge if

𝜇 \geq Δ 1 𝜌 2 + Δ 0 Δ 1 𝜌 2 - Δ 0 𝜈 .

(35)

On eliminating $Δ 0$ and $Δ 1$ using Eqs. , Eq. can be solved to yield the critical value of one of the parameters ( $𝜇, 𝜎, 𝜈, 𝜌$ , and $𝛾$ ) given the others. The resulting predictions for the point at which the mean abundance diverges is exact when the system is stable with respect to small perturbations, but is only approximate when the system becomes linearly unstable before the point of divergence is reached (see Sec. S7 in for numerical justification). We include both cases in Fig. . The stability conditions and provide a comprehensive analytical picture of stability in the hierarchical system. Similarly to Eq. , Eq. reduces to the corresponding stability criteria of Refs. in the nonhierarchical limit (Sec. S4.5.4 in ).

FIG. 5.

Hierarchical interactions and stability. (a) Similar to the phase diagrams in Ref. , we show the three possible dynamical behaviors of the system without hierarchy in the $𝜇 - 𝜎$ plane. We have $𝛾 = 0.8$ ( $20 %$ predator-prey interactions), $𝜈 = 0$ , and $𝜌 = 1$ . In (b) and (c) we choose the points (b) $𝜇 = - 0.5$ and $𝜎 = 0.7$ and (c) $𝜇 = 0.5$ and $𝜎 = 0.8$ and show stability in the $𝜈 - 𝜌$ plane. The left circle (yellow) in (a) and the circle in (b) are the same point in parameter space ( $𝜇 = - 0.5, 𝜈 = 0, 𝜎 = 0.7, 𝜌 = 1$ , and $𝛾 = 0.8$ ) and similarly for the right (red) circle in (a) and the circle at the center of panel (c) with $𝜇 = 0.5, 𝜈 = 0, 𝜎 = 0.8, 𝜌 = 1$ , and $𝛾 = 0.8$ .

B. Phase diagrams

The phase diagrams in Figs. and illustrate the effects of hierarchical interaction ( $𝜈$ and $𝜌$ ) on stability. When $𝜈$ is sufficiently large, increasing $𝜈$ typically decreases average abundances and the fraction of surviving species (see Fig. ), thereby pushing the system away from diverging abundances and from linear instability. Larger deviations of $𝜌$ from unity push the system both towards linear instability and generally towards infinite abundances. These effects are demonstrated in Figs. and , which indicate that a large enough value of $𝜌$ will always result in an unstable system and a large enough value of $𝜈$ will always result in a stable one. We also demonstrate the (separate) effects of changing $𝜈$ and $𝜌$ in Fig. . If our parameters are restricted such that either $𝜈 = 0$ or $𝜌 = 1$ [shown in Fig. ], then we can analytically demonstrate the effects of varying only one parameter (either $𝜈$ or $𝜌$ ) on linear instability. Details can be found in Sec. S6.3 in .

FIG. 6.

Phase diagrams (horizontal axis shows the proportion of predator-prey pairs). The onset of instability [(a) linear instability and (b) diverging abundances] are shown as black lines for $(𝜈, 𝜌) = (0, 1)$ , as well as for $(𝜈, 𝜌) = (3, 1)$ and $(𝜈, 𝜌) = (0, 3)$ as indicated. Stability is colored as labeled for the case $(𝜈, 𝜌) = (0, 1)$ . The lengths of arrows indicate the effect of increasing the values of $𝜈$ or $𝜌$ by one ( $𝜈 \to 𝜈 + 1$ or $𝜌 \to 𝜌 + 1$ , respectively).

Figures and also reveal how the precise combination of $𝜈$ and $𝜌$ can affect stability and community composition. Specifically, communities with $𝜈 > 0$ and $𝜌 > 1$ tend to have lower abundances than those with $𝜈 > 0$ and $𝜌 < 1$ and are further from the diverging abundances transition.

From Fig. we find that the influence of $𝜈$ on linear stability is relatively small compared with its influence on the point at which abundances diverge. This is reminiscent of the fact that the overall average interaction strength $𝜇$ has no effect on linear stability but affects the average abundance. Figure also reveals that $𝜌$ has a more dramatic effect on stability in communities with a large proportion of predator-prey pairs. In particular, any value of $𝜌 \neq 1$ removes the special behavior of the model without hierarchy and all interactions of the predator-prey type ( $𝜌 = 1, 𝜈 = 0$ , and $𝛾 = - 1$ ). In this special case, there is no linear instability , as indicated in Fig. . Any value $𝜌 \neq 1$ will however lead to the possibility of instability in the system with only predator-prey interactions. The stabilizing influence of $𝜈$ , on the contrary, is only mildly affected by the proportion of predator-prey pairs.

The influence of the remaining parameters $𝜇, 𝜎$ , and $𝛾$ on stability is relatively straightforward and similar to the model without block structure ( $𝜈 = 0$ and $𝜌 = 1$ ) previously studied in and in studies of complex ecosystems based on the spectra of random matrices . The parameter $𝜇$ has no effect on linear stability, as Eq. has no $𝜇$ dependence. Instead, from Eq. we see that $𝜇$ increases average abundances (noting that $Δ 0$ and $Δ 1$ are independent of $𝜇$ ), thereby moving the system towards instability via diverging abundances. This can be seen in Figs. and : A large enough value of $𝜇$ will always lead to diverging abundances in the community. Increasing $𝜎$ , on the other hand, pushes the system towards both linear instability and diverging abundances, as demonstrated in Figs. and . Increasing the correlation parameter $𝛾$ (i.e., decreasing the proportion of predator-prey interaction pairs) also pushes the system towards linear instability and diverging abundances, as demonstrated in Fig. .

VII. CONCLUSION

Our analysis of the generalized Lotka-Volterra model with hierarchical interactions has focused on the effect of hierarchy on both stability and structure in complex ecological communities. We have extended previous work on the stabilizing impact of predator-prey-like relationships by considering both the average severity of the hierarchy ( $𝜈$ ) and the proportion of interactions of predator-prey type ( $𝛾$ ) in a single dynamical model. We found that increases to both factors are stabilizing. We also found that increased heterogeneity in interaction variances ( $𝜌$ ) is a destabilizing force.

The dynamic mean-field theory approach, unlike an approach based on the spectra of random matrices, guarantees a feasible equilibrium and gives access to properties of the ecosystem other than stability. We found that communities with a strong hierarchy are dominated by species at the top, which are both more abundant and more likely to survive asymptotically. Further, hierarchy leads to more complex, non-Gaussian abundance distributions.

In order to find fixed-point equations for the hierarchical model with an infinite number of subcommunities, we first considered a related and more general community, divided into a finite number of subcommunities. Fixed-point equations were then obtained, resulting in an effective abundance for a representative species $𝑥 𝑎 (𝑡)$ in each community. This is in contrast to most similar studies employing dynamic mean-field theory , which do not have a structured population and accordingly only a single effective species. Our more general approach could allow for progress to be made investigating more heterogeneous interaction structures in and beyond ecology, such as other block-structured interaction matrices , metapopulation models on complex networks , or trophic levels .

Codes can be found in , as well as the method for simulating the community dynamics with Eq. , numerical solution procedures for solving the hierarchical fixed point equations , and the data and code for producing all figures.

ACKNOWLEDGMENTS

We acknowledge funding from the Spanish Ministry of Science, Innovation and Universities, the Agency AEI, and FEDER (EU) under the grant PACSS (No. RTI2018-093732-B-C22), the Maria de Maeztu program for Units of Excellence in R&D (Program No. MDM-2017-0711) funded by Grant No. MCIN/AEI/10.13039/501100011033, and the Engineering and Physical Sciences Research Council UK, Grant No. EP-T517823-1.

Generalized Lotka-Volterra model with hierarchical interactions

Abstract

Physics Subject Headings (PhySH)

Article Text

I. INTRODUCTION

II. MODEL DEFINITION

A. General block-structured interactions

B. Hierarchical interactions

III. DYNAMIC MEAN-FIELD THEORY

IV. FIXED-POINT EQUATIONS

A. General block-structured interactions

B. Hierarchical interactions

V. FIXED-POINT DISTRIBUTIONS

A. Abundance distributions

B. Fraction of surviving species

VI. STABILITY

A. Stability conditions

1. Linear instability

2. Diverging abundances

B. Phase diagrams

VII. CONCLUSION

ACKNOWLEDGMENTS

Supplemental Material

References (59)