# Spatial heterogeneity of interaction strength has contrasting effects on synchrony and stability in trophic metacommunities

Pierre Quévreux, Bart Haegeman, Michel Loreau

► To cite this version:

Pierre Quévreux, Bart Haegeman, Michel Loreau. Spatial heterogeneity of interaction strength has contrasting effects on synchrony and stability in trophic metacommunities. Ecology Letters, 2023, 10.1111/ele.14297. hal-03829838v3

HAL Id: hal-03829838

https://hal.science/hal-03829838v3

Submitted on 27 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spatial heterogeneity of interaction strength has contrasting effects on synchrony and stability in trophic metacommunities

Pierre Quévreux¹, Bart Haegeman¹, and Michel Loreau¹

¹Theoretical and Experimental Ecology Station, UAR 2029, CNRS, 09200 Moulis, France

## Abstract

Spatial heterogeneity is a fundamental feature of ecosystems, and ecologists have identified it as a factor promoting the stability of population dynamics. In particular, differences in interaction strengths and resource supply between patches generates an asymmetry of biomass turnover with a fast and a slow patch. The coupling of these two energy channels by mobile predators has been identified to increase stability at different scales by promoting the asynchrony of population dynamics between each patch. Here, we demonstrate that asymmetry has a contrasting effect on the stability of metacommunities receiving localised perturbations. We built a model of an asymmetric metacommunity with two patches linked by the dispersal of predators and in which prey receive stochastic perturbations only in one patch. Perturbing prey in the fast patch synchronises the dynamics of prey biomass between the two patches and destabilises predator dynamics by increasing their temporal variability. Conversely, perturbing prey in the slow patch decreases the synchrony of their dynamics and stabilises predator dynamics. This discrepancy between the responses is due to the asymmetric transmission of perturbations caused by the different distributions of biomass between the fast and the slow patch. Consequently, the fast patch drives the dynamics of the metacommunity and imposes synchrony while the slow patch does not. Therefore, local perturbations can have opposite consequences at the regional scale depending on the characteristics of the perturbed patch. Our results have strong implications for conservation ecology and suggest reinforcing protection policies in fast patches to dampen the effects of perturbations and promote the stability of population dynamics at the regional scale.

## Key words

source-sink, stochastic perturbations, food chain, dispersal, asymmetry, conservation

## Introduction

Since May (1972) demonstrated that stability was not an inherent property of ecological interaction networks, ecologists have been relentlessly looking for the mechanisms ensuring ecosystem stability. Spatial heterogeneity has long been identified as one of the main factors promoting the mechanisms underlying the maintenance of biodiversity and the stability of ecosystems. For instance, in competitive metacommunity models, spatial heterogeneity provides local favourable conditions to each species of the regional pool (Holt, 1984; Cheson, 2000; Amarasekare and Nisbet, 2001), which in turn ensures species persistence in less favourable patches by source-sink dynamics (Mouquet and Loreau, 2002, 2003; Loreau et al., 2003). The stability of the temporal dynamics of species biomass is ensured by the asynchrony of the dynamics between patches, which leads to compensatory dynamics

(Loreau et al., 2003; Loreau and de Mazancourt, 2008). In trophic metacommunities, spatial heterogeneity has also been identified as a stabilising factor (Steele, 1974; Hastings, 1977, 1978), but the underlying mechanisms are more complex due to the interplay between trophic and spatial dynamics.

Inspired by the description of fast and slow energy channels by soil ecologists (i.e., in terms of biomass turnover), Rooney et al. (2006) noted the stabilising effect of the asymmetry of energy flows in ecosystems with a food web model consisting of one mobile predator feeding on two energy channels. In their model, the asymmetry of energy flow is generated by different interaction strengths between predators and prey (i.e. increased attack rate in one energy channel compared to the other one, see Figure 1) and different consumption rates of a common resource by the two basal species, which in turn promotes the asynchrony of prey biomass dynamics in response to perturbations. Although synchrony patterns are tightly linked to stability patterns, because the asynchrony of local population dynamics leads to more stable dynamics (low biomass
<table><caption>Table 1: Approximative relative increase in predation risk between low-risk and high-risk environments (equivalent to the asymmetry of interaction strength $\gamma$ in Figure 1). See Gorini et al. (2012) for an extended review and more references.</caption><thead><tr><th>Predator</th><th>Prey</th><th>γ</th><th>Reference</th></tr></thead><tbody><tr><td>American marten</td><td>Vole species</td><td>1.6</td><td>Andruskiw et al., 2008</td></tr><tr><td>Wolf</td><td>Moose</td><td>14-100</td><td>Gervasi et al., 2013</td></tr><tr><td>Wolf</td><td>Roe deer</td><td>2.5-8</td><td>Gervasi et al., 2013</td></tr><tr><td>Wolf</td><td>Elk</td><td>10</td><td>Kauffman et al., 2007</td></tr><tr><td>Savannah predators</td><td>Savannah ungulates</td><td>1.5-4.5</td><td>Thaker et al., 2011</td></tr><tr><td>Artificial gecko</td><td>Australian predators</td><td>2.8</td><td>Hansen et al., 2019</td></tr><tr><td>Lynx</td><td>Roe deer</td><td>2</td><td>Gehr et al., 2020</td></tr><tr><td>Puma</td><td>Vicuña</td><td>1.6</td><td>Donadio and Buskirk, 2016</td></tr></tbody></table>

CV) at higher scales due to compensatory dynamics (Loreau et al., 2003; Gonzalez and Loreau, 2008; Loreau and de Mazancourt, 2013; Wilcox et al., 2017), subsequent studies suggested that increased asymmetry does not necessarily leads to increased stability. For example Ruokolainen et al. (2011) presented a model in which biomass fluctuations can become more variable with increasing asynchrony. Hence, the relationship between asymmetry and stability is not trivial and the mechanisms governing asynchrony through the difference in energy flow between the fast and slow channels are not well understood. To fill this gap, we propose to consider the effects of asymmetry from the metacommunity perspective since recent theoretical studies were able to accurately explain the synchrony and stability patterns in metacommunities (Quévreux et al., 2021a,b).

Metacommunities embody the spatial dimension of interaction networks: they consist of distant patches connected by the dispersal of the organisms living in each patch (Leibold et al., 2004; Leibold and Chase, 2017). The metacommunity framework is particularly suitable to represent the spatial heterogeneity observed in ecosystems because each community has its own characteristics such as biomass turnover. Following Rooney et al.'s (2006) model, many studies implemented spatial heterogeneity through the asymmetry of interaction strength and/or resource supply to manipulate the difference in biomass turnover between the energy channels hosted by each patch (Goldwyn and Hastings, 2009; Ruokolainen et al., 2011; Anderson and Fahimipour, 2021). In particular, interaction strength is key in community dynamics because it governs food web structure, stability (Neutel et al., 2002) and biomass distribution (Barbier and Loreau, 2019) by simultaneously determining predator growth and prey mortality. Therefore, its significant variations observed in nature, often reported as predation risk by prey in field studies (Table 1), should lead to dramatic variations in community functioning across space.

In addition to the asymmetry of interaction strength, Rooney et al. (2006) highlighted the im-

portance of mobile predators coupling two different energy channels, a keystone role in ecosystem functioning largely reported by empirical studies (Schindler and Scheuerell, 2002; Vadeboncoeur et al., 2005; Schmitz, 2004; Olff et al., 2009; Dolson et al., 2009; Schmitz et al., 2010). In Rooney et al.'s (2006) model, the perturbation of the mobile predator leads to an asynchronous response of prey, which stabilises the food web. However, Quévreux et al. (2021a) showed that the perturbation and dispersal of particular trophic levels govern synchrony and stability in symmetric metacommunities. In asymmetric metacommunities, the perturbation of particular patches should lead to different synchrony and stability patterns at the metacommunity scale because of the different dynamics in each patch in response to perturbations. In parallel to the keystone role of mobile predators, keystone communities (sensus Mouquet et al. (2013), which are equivalent to keystone patches), should have a major influence on synchrony and stability patterns. Therefore, we expect that asymmetry is not a generic stabilising factor, as claimed by Rooney et al. (2006), but strongly depends on which patch is perturbed according to its characteristics. To explore this statement, we consider a simple metacommunity model of two patches hosting a predator-prey couple and with asymmetric interaction strength and resource supply. The stability of the metacommunity is assessed by the response at different scales (e.g. CV of the biomass of a species at the local and regional scales) when prey receive stochastic perturbations in one of the two patches. We show contrasting effects of asymmetry: perturbing prey in the fast patch (equivalent to the fast channel defined by Rooney and McCann, 2012) promotes prey synchrony and decreases predator stability at the metapopulation scale while perturbing the slow patch has the opposite effects.
<figure><img src="image_2.png" /><figcaption>Figure 1: The metacommunity model consists of two patches, each sustaining a predator-prey couple linked by predators, which disperse at a very high scaled rate $d_2$. Prey grow at a rate $g/D$ and have a positive effect $\epsilon a$ on predators, while predators have a negative effect $ma$ on prey. Each species population is also limited by self-regulation $D$ (negative intraspecific interactions). Spatial heterogeneity is embodied by the asymmetry of resource supply (green) and the interaction strength (purple), which are higher in patch #1 by factors $\omega$ and $\gamma$, respectively. Consistent with Rooney et al. (2006), patch #1 is called the fast patch, and patch #2 is called the slow patch. Prey receive stochastic perturbations either in patch #1 or in patch #2 (red arrows).</figcaption></figure>

$$ \frac{1}{D} \frac{dB_1^{(1)}}{dt} = B_1^{(1)} \left( \omega \frac{g}{D} - B_1^{(1)} - \gamma ma B_2^{(1)} \right) \qquad (1a) $$

$$ \frac{1}{mD} \frac{dB_2^{(1)}}{dt} = \underbrace{B_2^{(1)} (-B_2^{(1)} + \gamma \epsilon a B_1^{(1)})}_{\text{intra-patch dynamics}} + \underbrace{d_2 (B_2^{(2)} - B_2^{(1)})}_{\text{dispersal}} \qquad (1b) $$

# Methods

## Metacommunity model

We use the model proposed by Quévreux et al. (2021a) based on the food chain model developed by Barbier and Loreau (2019). The model consists of two patches that each sustain a food chain with Lotka-Volterra predator-prey interactions (equations (1a) and (1b)). $B_1^{(k)}$ and $B_2^{(k)}$ are the biomasses of prey and predators, respectively, in patch $k$. Prey have a positive effect $\epsilon a$ on predators ($\epsilon$ is the conversion efficiency and $a$ is the interspecific interaction rate relative to intraspecific interactions), and predators have a negative effect $ma$ on prey ($m$ is the predator to prey metabolic rate ratio)(Figure 1). The time scale of the system is rescaled by the metabolic rate of prey, and biologic rates of each species $i$ are rescaled by its intraspecific interaction rate $D_i$. Therefore, we obtain the relative growth rate $g/D$ and the scaled dispersal rate $d_i$. Considering scaled parameters and aggregated parameters ($\epsilon a$ and $ma$) enables us to explore a wide range of ecological situations. We refer to Appendix S1-1 for a detailed description of the food chain model and analysis methods. All the necessary information to fully understand our results are in the main text and the supplementary

information only serves to give additional technical elements to fully reproduce our work and proofs of the robustness of our results. Parameters and their values are summarised in Table 2.

We reproduce the two main features of Rooney et al.’s (2006) model. First, predators disperse at a very high scaled rate $d_2 = 10^6$, while prey are immobile ($d_1 = 0$), and strongly couple the two patches. Slightly mobile prey ($0 < d_1 \ll d_2$) should not change the results because Quévreux et al. (2021a) showed that the species for which dispersal has the strongest influence drives the coupling between the two patches. Second, resource supply and interaction strength are asymmetric between patches since they are higher in patch #1 by factors $\gamma$ and $\omega$ respectively (Figure 1). Patch #1 corresponds to the fast energy channel, in which biomass has a high turnover, while patch #2 corresponds to the slow channel. Therefore, we call patch #1 the fast patch and patch #2 the slow patch. We set $\gamma = \omega$ to ensure species persistence over the entire range of parameters (see Figure S2-6 in the supporting information) but varying them independently does not qualitatively change the results (see Figure S2-14 in the supporting information). In the following, we only refer to $\gamma$ for the sake of simplicity and only consider $\gamma \ge 1$ because $\gamma \le 1$ just swaps the roles of patches #1 and #2.
<table><caption>Table 2: Table of parameters. $\sigma_i$ is set very small to keep the system in the vicinity of equilibrium. More combinations of $\varepsilon a$ and $ma$ are tested in the supporting information. $d_2$ is set very high to emphasise the high mobility of predators and their ability to couple prey populations. $\omega$ is set equal to $\gamma$. $r = 0$ removes the energetic limitations of the food chain and makes interactions the dominant factors determining biomass distribution and stability patterns, as in Barbier and Loreau (2019).</caption><thead><tr><th>Parameter</th><th>Interpretation</th><th>Value</th></tr></thead><tbody><tr><td>σ<sub>i</sub></td><td>standard deviation of stochastic noise</td><td>10<sup>-3</sup></td></tr><tr><td>g</td><td>net growth rate of prey</td><td>1</td></tr><tr><td>r</td><td>death rate of predators</td><td>0</td></tr><tr><td>D</td><td>self-regulation</td><td>1</td></tr><tr><td>ε</td><td>conversion efficiency</td><td>0.65</td></tr><tr><td>m</td><td>predator/prey metabolic rate ratio</td><td>0.65</td></tr><tr><td>a</td><td>attack rate</td><td>1.54</td></tr><tr><td>εa</td><td>positive effect of prey on predators</td><td>1</td></tr><tr><td>ma</td><td>negative effect of predators on prey</td><td>1</td></tr><tr><td>d<sub>2</sub></td><td>scaled dispersal rate of predators</td><td>10<sup>6</sup></td></tr><tr><td>ω</td><td>asymmetry of resource supply</td><td>[1,10]</td></tr><tr><td>γ</td><td>asymmetry of interaction strength</td><td>[1,10]</td></tr></tbody></table>

## Response to stochastic perturbations

We use the same methods as Quévreux et al. (2021a) to study the response of metacommunities to stochastic perturbations. Indeed, recent studies advocate for the use of the temporal variability of biomass (Haegeman et al., 2016; Arnoldi et al., 2018), which is measured by the coefficient of variation (CV), and can be easily measured experimentally. In addition, Wang and Loreau (2014, 2016), Wang et al. (2019), and Jarillo et al. (2022) showed that CVs scale up from local populations to community, regional and metacommunity levels, therefore providing a comparison of stability at different scales. Here, we provide only a brief description of the main concepts, but a thorough description is available in Appendix S1.

Prey in the fast or slow channel receive stochastic perturbations that are represented by equation (2).

$$dB_i = \underbrace{f_i(B_1, \dots, B_S)dt}_{\text{Deterministic}} + \underbrace{\sigma_i \sqrt{B_i^*} dW_i}_{\text{Perturbation}} \quad (2)$$

$f_i(B_1, \dots, B_S)$ represents the deterministic part of the dynamics of species $i$, as described by equations (1a) and (1b)). Stochastic perturbations are defined by their standard deviation $\sigma_i$ and $dW_i$, a white noise term with a mean of 0 and variance of 1. Perturbations also scale with the square root of the biomass at equilibrium $B_i^*$ of the perturbed population. Such scaling makes the perturbations similar to demographic stochasticity (from birth-death processes) that evenly affect each species regardless of abundance (Arnoldi et al., 2019). In other words, the ratio of mean species biomass variance to perturbation variance is roughly independent of biomass, which disentangles the effect of asymmetry on perturbation transmission from its effect on species abundance. Therefore, for different perturbations

affecting different species with the same value of standard deviation $\sigma_i$, we generate a similar variance at the metacommunity scale regardless the abundance of the perturbed species and excite the entire metacommunity with the same intensity (see Figure S2-3 in the supporting information).

In the following, we assess the temporal variability of the biomass of each population induced by stochastic perturbations affecting the metacommunity. Therefore, we linearise the system in the vicinity of equilibrium to obtain equation (3) where $X_i = B_i - B_i^*$ is the deviation from equilibrium.

$$\frac{d\vec{X}}{dt} = J\vec{X} + T\vec{E} \quad (3)$$

$J$ is the Jacobian matrix, which represents the linearised direct effects between populations in the vicinity of equilibrium, and $T$ defines how the perturbations $E_i = \sigma_i dW_i$ apply to the system (i.e., which species they affect and how they scale with biomass, where $T$ is a diagonal matrix whose terms are $T_{ii} = \sqrt{B_i^*}$).

Because the system is at steady state, the stationary variance-covariance matrix $C^*$ of species biomasses (variance-covariance matrix of $\vec{X}$, see the demonstration in Appendix S1-5) can be calculated from the variance-covariance matrix of perturbations $V_E$ (variance-covariance matrix of $\vec{E}$) by solving the Lyapunov equation (4) (Arnold, 1974; Wang et al., 2015; Arnoldi et al., 2016; Quévreux et al., 2021a).

$$JC^* + C^*J^\top + TV_E T^\top = 0 \quad (4)$$

The expressions for $V_E$ and $T$ and the method to solve the Lyapunov equation are detailed in Appendix S1-6. From the variance-covariance matrix $C^*$, we compute the coefficient of correlation of the biomass dynamics between the two populations of
each species (see equation (22) in the supporting information) and we measure the stability with the coefficient of variation (CV) of the biomass. In addition, biomass CVs can be measured at different scales: population scale (e.g., biomass CV of prey in patch #1), metapopulation scale (e.g., CV of the total biomass of prey) and metacommunity scale (e.g., CV of the total biomass of predator and prey put together) to assess the effects of asymmetry at local and regional scales (Figure 3A and see Appendix S1-7). Finally, we quantify the synchrony of the dynamics of the different populations with the coefficient of correlation, which is also computed from the variance-covariance matrix $C^*$ (Appendix S1-7).

# Results

## Effects on stability

We describe how the asymmetry of interaction strength $\gamma$ shapes metacommunity stability at different scales. Since predators have a very high scaled dispersal rate ($d_2 = 10^6$), their populations are perfectly correlated and display the same dynamics. Our main result is that prey become more correlated when they are perturbed in patch #1 (fast channel in which $\gamma > 1$), while they become more anticorrelated when they are perturbed in patch #2 (Figure 2). Increasing $\gamma$ amplifies the difference in correlation between these two scenarios, and this pattern qualitatively holds for various combinations of the physiological and ecological parameters $\varepsilon a$ and $ma$ (see Figure S2-9 in the supporting information).

Increasing the asymmetry of interaction strength $\gamma$ has contrasting effects on biomass CV at different scales as well (Figure 3A). At the population scale, it increases the biomass CV of each population when prey are perturbed in the fast patch (Figure 3B). When prey are perturbed in the slow patch, increasing $\gamma$ slightly alters the biomass CV of prey in patch #1, increases the biomass CV of prey in patch #2 and decreases the biomass CV of predators. This discrepancy can be attributed to the strong effect of $\gamma$ on prey biomass in patch #2 (Figure 4): prey biomass strongly decreases with $\gamma$ in patch #2, which increases their biomass CV.

At the metapopulation scale, the asymmetry of interaction strength $\gamma$ increases the biomass CV of prey in both scenarios of perturbation (Figure 3C). However, this result is not true for all values of $\varepsilon a$ and $ma$ (Figure S2-10A in the supporting information) because of the various responses of prey biomass to $\gamma$ among patches (Figure S2-8A in the supporting information). The biomass CV of predators is higher when prey are perturbed in the fast patch (patch #1) compared to the case in which prey are perturbed in the slow patch (#2) (Figure 3C), which is consistent for all values of $\varepsilon a$ and

<figure><img src="image_3.png" /></figure>

Figure 2: Spatial correlation between the populations of each species depending on asymmetry of interaction strength $\gamma$ when predators disperse and prey are perturbed in patch #1 or #2. The dashed line emphasises the value of the correlation of prey populations without asymmetry ($\gamma = 1$). Note that the curves for predators overlap because their high dispersal that perfectly correlates their dynamics regardless of the perturbed patch.

$ma$ (see Figure S2-10A in the supporting information).

Finally, stability at the metacommunity scale depends on the distribution of biomass and CV among species. In our particular case ($\varepsilon a = 1$ and $ma = 1$), predators have the largest total biomass (Figure 4) and drive the biomass CV at the metacommunity scale for low values of asymmetry of interaction strength $\gamma$ (Figure 3D). For high values of $\gamma$, when prey are perturbed in patch #2, the CV of total biomass increases with $\gamma$ because it is driven by prey in patch #2, whose biomass CV is much higher than the biomass CV of predators, which compensates for their lower biomass. Other values of $\varepsilon a$ and $ma$ lead to other distributions of biomass and CV among species, which can make prey to drive the stability at the metacommunity scale (see Figures S2-8 and S2-10 in the supporting information).

## Underlying mechanisms

To unveil the mechanisms governing the stability of heterogeneous metacommunities, we look deeper into the dynamics after a pulse perturbation (Figure 5A) and explain them with the direct effects between species quantified by the Jacobian matrix (see equation (3)). When the perturbation of prey occurs in patch #1, the strong direct effect of prey on predators (and vice versa) in patch #1 due to $\gamma$ (Figure 5B) leads to a strong response of predators ①, which in turn drives the response of the two prey populations ②. In detail, predator biomass in patch #1 first increases because of the abundance of prey. Then, predators deplete prey biomass in
A

<figure><img src="image_4.png" /></figure>

B

<figure><img src="image_5.png" /></figure>

C

<figure><img src="image_6.png" /></figure>

D

<figure><img src="image_7.png" /></figure>

Figure 3: Stability at different scales depending on asymmetry of interaction strength $\gamma$ when predators disperse and prey are perturbed in patch #1 or #2. A) The temporal variability in the metacommunity is assessed by the coefficient of variation (CV) of biomass at different scales: population scale (biomass CV of one species in one patch), metapopulation scale (CV of the total biomass of one species across patches) and metacommunity scale (CV of the total biomass of the entire metacommunity). B) Biomass CV at the population scale. Note that the curves for predators overlap because their high dispersal perfectly balances their biomass distribution between the two patches. C) Biomass CV at the metapopulation scale. D) Biomass CV at the metacommunity scale.

<figure><img src="image_8.png" /></figure>

Figure 4: Distribution of the biomass of each species among patches depending on the asymmetry of interaction strength $\gamma$. Note that the curves for predators overlap because their high dispersal that perfectly balances their biomass between the two patches.

In turn, the very low direct effect of predators on prey in patch #2 does not allow perturbations to ripple back to patch #2 where prey slowly respond from the initial perturbation (4) (Figure 5B). This slow response is emphasised by the source-sink dynamics in the metacommunity (Figure 4 and Figure S2-5 in the supporting information), which leads to a lower biomass of prey in patch #2 compared to a metacommunity without dispersal, therefore decreasing biomass flows in patch #2 and its response speed. This difference in response speed between patches #1 and #2 leads to the anti-correlation of prey populations because it increases the time interval in which they have opposite variations: an increase in the biomass of prey in patch #1, which follows the initial decrease due to predation, and a slow decrease in prey biomass in patch #2.

both patches and correlate their dynamics, which explains why asymmetry of interaction strength $\gamma$ increases prey correlation when prey are perturbed in patch #1.

When the perturbation of prey occurs in patch

## Discussion

We have shown that the asymmetry of interaction strength has contrasting effects on stability depending on which patch is perturbed. Perturbing prey in the fast patch (in which interaction strength
<figure><img src="image_9.png" /><figcaption>Figure 5: **A)** Time series of biomasses rescaled by their value at equilibrium after an increase of 20% in prey biomass in patch #1 (left panel) or patch #2 (right panel) for a value of interaction strength asymmetry $\gamma = 3$. **B)** Direct effect of prey on predator (and vice versa) depending on interaction strength asymmetry $\gamma$. Direct effects correspond to the terms of the Jacobian matrix and the dashed line emphasises the null direct effects. On the central schema, the numbers are the numeric values of the terms of the Jacobian matrix corresponding to each arrow for $\gamma = 3$.</figcaption></figure>

is the highest) tends to synchronise the dynamics of prey biomass and increases the temporal variability of predator dynamics at the metapopulation scale, while perturbing prey in the slow patch decreases the synchrony of prey dynamics and the temporal variability of predator dynamics. This discrepancy between the responses is due to asymmetric transmission of perturbations within each patch, itself caused by different biomass distributions and interaction strengths between the two patches. Perturbations are strongly transmitted from the fast patch to the slow patch, while the reverse transmission is weak. Consequently, the fast patch drives the dynamics of the metacommunity and synchronises prey dynamics, while the slow patch does not, and the almost independent respond from perturbation in each patch decreases the synchrony of prey dynamics.

## Stability in a heterogeneous world

Our results show that spatial heterogeneity, which is represented by the asymmetry of interaction strength and resource supply as in Rooney et al. (2006), generates mechanisms that alter local and regional dynamics, which deeply changes the synchrony of population dynamics and the stability of the metacommunity at different scales. Quévreux et al. (2021a) showed that, in a homogeneous metacommunity, the spatial correlations between patches can be obtained from the within-patch correlations, the dispersing species making the link between the two (see Figure S2-28 for a summary of the results of Quévreux et al., 2021a). In other words, knowledge of the dynamics at the local scale is enough to understand the stability pattern at the metacommunity scale. In a heterogeneous metacommunity, a similar approach does not
work because patches do not contribute equally to the dynamics. In a heterogeneous metacommunity, a similar approach does not work because patches do not contribute equally to the dynamics. In particular, a patch with fast energy flow can have an overwhelming impact (see Figure S2-7 in the supporting information). Clearly, the dynamics at the metacommunity scale cannot be assessed by the dynamics at the local scale, as in Quévreux et al. (2021a), and they are an emergent property resulting from the tight interplay between the strength of perturbation transmission in each patch.

Rooney et al. (2006) verbally explained that each patch has dynamics with different speeds: the fast channel (with higher interaction rates and resource supplies) enables a quick response after a perturbation while the slow channel dampens the dynamics in the long term and prevents the system from overshooting. By considering the stability at different scales, our results contrast with this explanation. On the one hand, the asynchrony of prey dynamics, when they are perturbed in patch #2 (Figure 2), stabilises the dynamics of predators because their resource supplies are asynchronous. On the other hand, the dynamics of prey at the metapopulation scale are not stabilised by their asynchrony (Figure 3C) because of the low local stability in patch #2 (Figure 3B), which decreases the overall stability of prey. The potential stabilising effect of asymmetry depends both on the perturbed patch and the considered trophic level. Therefore, the overall stability at the metacommunity scale is governed by the relative contributions of the various populations in response to local perturbations, and asymmetry *per se* does not have a stabilising effect.

Our description of the mechanisms underlying the apparent stabilising effects of spatial heterogeneity should enlighten the results of previous theoretical studies. Goldwyn and Hastings (2009) and Ruokolainen et al. (2011) found that the asymmetry of interaction rate leads to asynchrony by generating out-of-phase dynamics in a system with endogenous oscillations. In particular, Ruokolainen et al. (2011) found a U-shaped relationship: for moderate asymmetry, the spatial asynchrony of predator and prey populations is maximal, which leads to optimum stability at the metacommunity scale. Our results suggest that moderate asymmetry would alter the phase of the oscillations in each patch while keeping the amplitude of oscillations equivalent, therefore promoting asynchrony. Conversely, a strong asymmetry would increase the imbalance between oscillation amplitude and enable the fast patch to take over the slow patch, which would bring back synchrony. However, their results rely on phase-locking (Jansen, 1999; Lloyd and May, 1999; Goldwyn and Hastings, 2008; Vasseur and Fox, 2009), which is the coupling of the phase of oscillators embodied by predator-prey pairs in each patch. Although our results provide interesting insight into metacommunity dynamics, they can-

not grasp the fine mechanisms underlying nonlinear phenomena such as phase-locking and further studies are needed to identify these mechanisms.

## Generality of the effects of asymmetry on stability

Our main results is that asymmetry is stabilising when the slow patch is perturbed, while it is destabilising when the fast patch is perturbed. This result is strikingly robust to several deviations from the original model we have described. First, we show that the described mechanisms are valid for a wide range of ecological and physiological parameters leading to various distributions of biomass among predators and prey (see Figures S2-8 and S2-9 in the supporting information). In addition, we observe the same results for longer food chains as long as prey populations are directly coupled by the dispersing predator (see Figures S2-18 and S2-19 in the supporting information). Currently, we do not identify a clear pattern for species lower in the food chain over a wide range of ecological and physiological parameters but further studies are needed to investigate the potential indirect effects propagating across the food chain. Second, the mechanisms are not restricted to prey populations coupled by a mobile predator but also apply to predator populations coupled by a mobile prey (see Figure S2-22 in the supporting information). Therefore, we anticipate that mobile predators are not the only major drivers of synchrony and stability in ecosystems (Schindler and Scheuerell, 2002; Vadeboncoeur et al., 2005; Dolson et al., 2009; Olff et al., 2009; Rooney and McCann, 2012), and resource species may also have an equivalent impact. Taken together, these two points strongly suggest that the mechanisms underlying stability and synchrony in response to perturbations should be general to metacommunities regardless of the ecological parameters, biomass distribution and dispersal among species.

Spatial heterogeneity has often been presented as a generic condition generating mechanisms ensuring stability, but our results contradict this statement. The models focusing on the asymmetric feeding of consumers on different energy channels or different patches showed that it promotes the existence of stable equilibria (McCann et al., 1998), greater asymptotic resilience (Rooney et al., 2006), asynchrony of prey in response to predator perturbation (Rooney et al., 2006) and out-of-phase limit cycles (Goldwyn and Hastings, 2009; Ruokolainen et al., 2011). All these studies considered measures of stability aiming to capture the general stability properties of metacommunities and miss the targeted effects of perturbations as we explained earlier. Although asymmetry does not necessarily promote stability, our results show that general mechanisms drive the response of metacommunities to localised perturbations, therefore providing a valuable framework to assess the response of ecosystems
to localised perturbations due to human activity. Additionally, these mechanisms enable us to understand the effect of environmental perturbations affecting all patches. As demonstrated by Arnoldi et al. (2019), environmental perturbations affect abundant populations the most, which is the prey population in the fast patch in our case (see Figure S2-25 in the supporting information). Therefore, we anticipate that the fast patch will govern the dynamics of metacommunities in which all populations are perturbed (see Figure S2-27 in the supporting information).

## Implications for conservation

The metacommunity framework has long been used in conservation ecology (Johnson et al., 2013; Schiesari et al., 2019; Patrick et al., 2021). Conservation efforts are usually concentrated on particular locations and useful management must consider the ecological processes acting at the landscape scale (Van Teeffelen et al., 2012; Chase et al., 2020). For instance, spatial heterogeneity is key to ensuring species coexistence and diversity at the regional scale, which ultimately provides important ecosystem services in agricultural landscapes (Bennett et al., 2006). A large corpus of theoretical studies explored the local response of communities in a landscape receiving perturbations (Mouquet et al., 2011; Economo, 2011; Holyoak et al., 2020; Jacquet et al., 2022). However, these studies focus on extinction events recovered by dispersal events in a patch dynamics framework, and little is known about the effect of moderate or small perturbations. In this context, the present study provides valuable insight into fine-scale dynamics in response to perturbations.

Our results show that species interactions are a major driver of synchrony in heterogeneous metacommunities. Even if the species of interest does not disperse significantly, the synchrony of the dynamics of its different populations can strongly depend on the interactions with another species with a higher dispersal across the landscape. For instance, Howeth and Leibold (2013) showed that predatory fish promote the asynchrony of oscillating populations of zooplankton in a mesocosm experiment. Therefore, species endorsing this role are called "mobile link organisms" (Lundberg and Moberg, 2003) and are particularly targeted by conservation policy because they have major impacts on community dynamics and ecosystem functioning (Soulé et al., 2005; Brodie et al., 2018). Such a species can be considered keystone species (Mills and Doak, 1993) and must be clearly identified to properly manage the conservation of the other interaction species. However, our results show that mobile link organisms are not the only driver of metacommunity stability, and the patch being perturbed also has a major impact. The concept of a keystone community, defined by Mouquet et al. (2013) for

communities whose destruction causes species extinction or a decrease in biomass production, can be applied to better assess the stability of metacommunities. Keystone communities are usually identified as those patches that are strongly connected to other patches in the spatial network (Resetarits et al., 2018), but our results suggest that the dynamical properties of each patch can be important as well. For instance, the fast patch can be identified as a keystone patch because of its ability to synchronise the dynamics of the other patches. Therefore, identifying the communities living in fast and slow patches should be key for conservation management aiming to mitigate the effects of perturbations.

According to our results, mitigating the effects of perturbations affecting the patch in which interaction strength is the highest is critical to avoid the synchrony of prey dynamics (Figure 2) and ensure predator stability (Figure 3C). Then, the patch in which the interaction strength between the species of interest and the mobile link organism is the highest must be identified. Conservation policies usually target preserved areas because they are characterised by high species richness but identifying them as fast or slow patches is not trivial. Urban ecology is a relevant example because many species dwell in cities and less anthropised ecosystems (e.g., agricultural and natural landscapes). Urban areas can be considered fast patches because of the abundance of resources (parameters $\omega$ in our model) for opportunistic species, but they can also be considered slow patches because of the reduced predation pressure (parameter $\gamma$ in our model), cities acting as safe spaces (see Shochat et al. (2006) and Shochat et al. (2010) for review). Typically, birds and rodents can find plenty of food due to human wastes, public parks and feeding while experiencing less predation (Rebolo-Ifrán et al., 2017). Therefore, focusing conservation efforts on urban areas to mitigate the perturbations affecting their ecosystem may be as important as protecting wild areas to protect species at the metapopulation scale.

## Conclusion

Asymmetry of interaction strength, and spatial heterogeneity in general, is not stabilising factor *per se* because perturbing prey in the fast patch leads to the synchrony of the dynamics of prey populations and increases the temporal variability of the mobile predator linking the two patches. Therefore, the response of metacommunities to perturbations is strongly context dependent, i.e., a good knowledge of the characteristics of each patch relative to each other is required to assess stability at the metacommunity scale. Based on our findings, we advocate for conservation efforts to target key patches not only according to species richness or biomass density but also according to the distribution of interaction strength across the metacommunity.
Acknowledgements

This work was supported by the TULIP Labora-
tory of Excellence (ANR-10-LABX-41) and by the
BIOSTASES Advanced Grant, funded by the Euro-
pean Research Council under the European Union's
Horizon 2020 research and innovation programme
(666971).

This paper has been peer-reviewed and rec-
ommended by PCI Ecology and Werner Ulrich
(10.24072/pci.ecology.100512). We thank the three
reviewers Phillip P.A. Staniczenko, Ludek Berec
and Diogo Provete for their constructive comments.

Data accessibility

The R codes to reproduce the results
and the figures are available on GitHub
(https://github.com/PierreQuevreux/model_
metacommunity_spatial_heterogeneity).

References

Amarasekare, P., & Nisbet, R. M. (2001). Spatial heterogeneity, source-sink dynamics, and the local coexistence of competing species. *The American Naturalist*, **158**(6), 572-584. https://doi.org/10.1086/323586

Anderson, K. E., & Fahimipour, A. K. (2021). Body size dependent dispersal influences stability in heterogeneous metacommunities. *Scientific Reports*, *11*(1), 17410. https://doi.org/10.1038/s41598-021-96629-5

Andruskiw, M., Fryxell, J. M., Thompson, I. D., & Baker, J. A. (2008). Habitat-mediated variation in predation risk by the american marten. *Ecology*, *89*(8), 2273-2280. https://doi.org/10.1890/07-1428.1

Arnold, L. (1974). *Stochastic differential equations: Theory and applications*. Wiley.

Arnoldi, J.-F., Bideault, A., Loreau, M., & Haegeman, B. (2018). How ecosystems recover from pulse perturbations: A theory of short- to long-term responses. *Journal of Theoretical Biology*, **436**, 79-92. https://doi.org/10.1016/j.jtbi.2017.10.003

Arnoldi, J.-F., Loreau, M., & Haegeman, B. (2016). Resilience, reactivity and variability: A mathematical comparison of ecological stability measures. *Journal of Theoretical Biology*, **389**, 47-59. https://doi.org/10.1016/j.jtbi.2015.10.012

Arnoldi, J.-F., Loreau, M., & Haegeman, B. (2019). The inherent multidimensionality of temporal variability: How common and rare species shape stability patterns (J. Chase, Ed.). *Ecology Letters*, **22**(10), 1557-1567. https://doi.org/10.1111/ele.13345

Barbier, M., & Loreau, M. (2019). Pyramids and cascades: A synthesis of food chain functioning and stability. *Ecology Letters*, **22**(2), 405-419. https://doi.org/10.1111/ele.13196

Bennett, A. F., Radford, J. Q., & Haslem, A. (2006). Properties of land mosaics: Implications for nature conservation in agricultural environments. *Biological Conservation*, **133**(2), 250–264. https://doi.org/10.1016/j.bioccon.2006.06.008

Brodie, J. F., Redford, K. H., & Doak, D. F. (2018). Ecological function analysis: Incorporating species roles into conservation. *Trends in Ecology & Evolution*, **33**(11), 840–850. https://doi.org/10.1016/j.tree.2018.08.013

Chase, J. M., Jeliazkov, A., Ladouceur, E., & Viana, D. S. (2020). Biodiversity conservation through the lens of metacommunity ecology. *Annals of the New York Academy of Sciences*, **1469**(1), 86–104. https://doi.org/10.1111/nyas.14378

Chesson, P. (2000). Mechanisms of maintenance of species diversity. *Annual Review of Ecology and Systematics*, **31**(1), 343–366. https://doi.org/10.1146/annurev.ecolsys.31.1.343

Dolson, R., McCann, K., Rooney, N., & Ridgway, M. (2009). Lake morphometry predicts the degree of habitat coupling by a mobile predator. *Oikos*, **118**(8), 1230–1238. https://doi.org/10.1111/j.1600-0706.2009.17351.x

Donadio, E., & Buskirk, S. W. (2016). Linking predation risk, ungulate antipredator responses, and patterns of vegetation in the high Andes. *Journal of Mammalogy*, **97**(3), 966–977. https://doi.org/10.1093/jmammal/gyw020

Economo, E. P. (2011). Biodiversity conservation in meta-
community networks: Linking pattern and persis-
tence. *The American Naturalist*, **177**(6), E167-
E180. https://doi.org/10.1086/659946

Gehr, B., Bonnot, N. C., Heurich, M., Cagnacci, F., Ciuti, S., Hewison, A. J. M., Gaillard, J.-M., Ranc, N., Premier, J., Vogt, K., Hofer, E., Ryser, A., Vimercati, E., & Keller, L. (2020). Stay home, stay safe—Site familiarity reduces predation risk in a large herbivore in two contrasting study sites (L. Prugh, Ed.). *Journal of Animal Ecology*, **89**(6), 1329–1339. https://doi.org/10.1111/1365-2656.13202

Gervasi, V., Sand, H., Zimmermann, B., Mattisson, J., Wabakken, P., & Linnell, J. D. C. (2013). Decomposing risk: Landscape structure and wolf behavior generate different predation patterns in two sympatric ungulates. *Ecological Applications*, **23**(7), 1722–1734. https://doi.org/10.1890/12-1615.1

Goldwyn, E. E., & Hastings, A. (2008). When can dispersal synchronize populations? *Theoretical Population Biology*, **73**(3), 395–402. https://doi.org/10.1016/j.tpb.2007.11.012

Goldwyn, E. E., & Hastings, A. (2009). Small heterogeneity has large effects on synchronization of ecological oscillators. *Bulletin of Mathematical Biology*, **71**(1), 130–144. https://doi.org/10.1007/s11538-008-9355-9

Gonzalez, A., & Loreau, M. (2008). The causes and consequences of compensatory dynamics in ecological communities. *Annual Review of Ecology, Evolution, and Systematics*, **40**(1), 393–414. https://doi.org/10.1146/annurev.ecolsys.39.110707.173349

Gorini, L., Linnell, J. D. C., May, R., Panzacchi, M., Boitani, L., Odden, M., & Nilsen, E. B. (2012). Habitat heterogeneity and mammalian predator-prey interactions: Predator-prey interactions in a spatial world. *Mammal Review*, **42**(1), 55–77. https://doi.org/10.1111/j.1365-2907.2011.00189.x

Haegeman, B., Arnoldi, J.-F., Wang, S., de Mazancourt, C., Montoya, J. M., & Loreau, M. (2016). Resilience, invariability, and ecological stability across levels of organization. *bioRxiv*. https://doi.org/10.1101/085852

Hansen, N. A., Sato, C. F., Michael, D. R., Lindenmayer, D. B., & Driscoll, D. A. (2019). Predation risk for reptiles is highest at remnant edges in agricultural landscapes (T. M. Lee, Ed.). *Journal of Applied Ecology*, **56**(1), 31–43. https://doi.org/10.1111/1365-2664.13269

Hastings, A. (1977). Spatial heterogeneity and the stability of predator-prey systems. *Theoretical Population Biology*, **12**(1), 37–48. https://doi.org/10.1016/0040-5809(77)90034-X

Hastings, A. (1978). Spatial heterogeneity and the stability of predator-prey systems: Predator-mediated coexistence. *Theoretical Population Biology*, **14**(3), 380–
395. https://doi.org/10.1016/0040-5809(78)90015-1

Holt, R. D. (1984). Spatial heterogeneity, indirect interactions, and the coexistence of prey species. *The American Naturalist*, **124**(3), 377-406. https://doi.org/10.1086/284280

Holyoak, M., Caspi, T., & Redosh, L. W. (2020). Integrating disturbance, seasonality, multi-year temporal dynamics, and dormancy into the dynamics and conservation of metacommunities. *Frontiers in Ecology and Evolution, 8*, 571130. https://doi.org/10.3389/fevo.2020.571130

Howeth, J. G., & Leibold, M. A. (2013). Predation inhibits the positive effect of dispersal on intraspecific and interspecific synchrony in pond metacommunities. *Ecology*, **94**(10), 2220-2228. https://doi.org/10.1890/12-2066.1

Jacquet, C., Munoz, F., Bonada, N., Datry, T., Heino, J., & Jabot, F. (2022). Temporal variation of patch connectivity determines biodiversity recovery from recurrent disturbances. *bioRxiv*. https://doi.org/10.1101/2022.01.02.474736

Jansen, V. A. A. (1999). Phase locking: Another cause of synchronicity in predator–prey systems. *Trends in Ecology & Evolution*, **14**(7), 278–279. https://doi.org/10.1016/S0169-5347(99)01654-7

Jarillo, J., Cao-García, F. J., & De Laender, F. (2022). Spatial and ecological scaling of stability in spatial community networks. *Frontiers in Ecology and Evolution, 10*, 861537. https://doi.org/10.3389/fevo.2022.861537

Johnson, P. T. J., Hoverman, J. T., McKenzie, V. J., Blaustein, A. R., & Richgels, K. L. D. (2013). Urbanization and wetland communities: Applying metacommunity theory to understand the local and landscape effects (M. Cadotte, Ed.). *Journal of Applied Ecology*, **50**(1), 34–42. https://doi.org/10.1111/1365-2664.12022

Kauffman, M. J., Varley, N., Smith, D. W., Stahler, D. R., MacNulty, D. R., & Boyce, M. S. (2007). Landscape heterogeneity shapes predation in a newly restored predator–prey system. *Ecology Letters*, **10**(8), 690–700. https://doi.org/10.1111/j.1461-0248.2007.01059.x

Leibold, M. A., & Chase, J. M. (2017). *Metacommunity Ecology* (Vol. 59). Princeton University Press. https://doi.org/10.2307/j.ctt1wf4d24

Leibold, M. A., Holyoak, M., Mouquet, N., Amarasekare, P., Chase, J. M., Hoopes, M. F., Holt, R. D., Shurin, J. B., Law, R., Tilman, D., Loreau, M., & Gonzalez, A. (2004). The metacommunity concept: A framework for multi-scale community ecology. *Ecology Letters*, **7**(7), 601–613. https://doi.org/10.1111/j.1461-0248.2004.00608.x

Lloyd, A. L., & May, R. M. (1999). Synchronicity, chaos and population cycles: Spatial coherence in an uncertain world. *Trends in Ecology & Evolution*, **14**(11), 417–418. https://doi.org/10.1016/S0169-5347(99)01717-6

Loreau, M., Mouquet, N., & Gonzalez, A. (2003). Biodiversity as spatial insurance in heterogeneous landscapes. *Proceedings of the National Academy of Sciences*, **100**(22), 12765–12770. https://doi.org/10.1073/pnas.2235465100

Loreau, M., & de Mazancourt, C. (2008). Species synchrony and its drivers: Neutral and nonneutral community dynamics in fluctuating environments. *The American Naturalist*, **172**(2), E48–E66. https://doi.org/10.1086/589746

Loreau, M., & de Mazancourt, C. (2013). Biodiversity and ecosystem stability: A synthesis of underlying mechanisms. *Ecology Letters*, **16**, 106–115. https://doi.org/10.1111/ele.12073

Lundberg, J., & Moberg, F. (2003). Mobile link organisms and ecosystem functioning: Implications for ecosystem resilience and management. *Ecosystems*, **6**(1),

0087-0098. https://doi.org/10.1007/s10021-002-0150-4

May, R. M. (1972). Will a large complex system be stable? *Nature*, **238**(5364), 413–414. https://doi.org/10.1038/238413a0

McCann, K. S., Hastings, A., & Huxel, G. R. (1998). Weak trophic interactions and the balance of nature. *Nature*, **395**(6704), 794–798. https://doi.org/10.1038/27427

Mills, L. S., & Doak, D. F. (1993). The keystone-species concept in ecology and conservation. *BioScience*, **43**(4), 219–224. https://doi.org/10.2307/1312122

Mouquet, N., Gravel, D., Massol, F., & Calcagno, V. (2013). Extending the concept of keystone species to communities and ecosystems (B. Blasius, Ed.). *Ecology Letters*, **16**(1), 1–8. https://doi.org/10.1111/ele.12014

Mouquet, N., & Loreau, M. (2002). Coexistence in metacommunities: The regional similarity hypothesis. *The American Naturalist*, **159**(4), 420–426. https://doi.org/10.1086/338996

Mouquet, N., & Loreau, M. (2003). Community patterns in source-sink metacommunities. *The American Naturalist*, **162**(5), 544–557. https://doi.org/10.1086/378857

Mouquet, N., Matthiessen, B., Miller, T., & Gonzalez, A. (2011). Extinction debt in source-sink metacommunities (T. Romanuk, Ed.). *PLoS ONE*, **6**(3), e17567. https://doi.org/10.1371/journal.pone.0017567

Neutel, A.-M., Heesterbeek, J. A. P., & De Ruiter, P. C. (2002). Stability in real food webs: Weak links in long loops. *Science*, **296**(5570), 1120–1123. https://doi.org/10.1126/science.1068326

Olff, H., Alonso, D., Berg, M. P., Eriksson, B. K., Loreau, M., Piersma, T., & Rooney, N. (2009). Parallel ecological networks in ecosystems. *Philosophical Transactions of the Royal Society B: Biological Sciences*, **364**(1524), 1755–1779. https://doi.org/10.1098/rstb.2008.0222

Patrick, C. J., Anderson, K. E., Brown, B. L., Hawkins, C. P., Metcalfe, A., Saffarinia, P., Siqueira, T., Swan, C. M., Tonkin, J. D., & Yuan, L. L. (2021). The application of metacommunity theory to the management of riverine ecosystems. *WIREs Water*, **8**(6). https://doi.org/10.1002/wat2.1557

Quévreux, P., Barbier, M., & Loreau, M. (2021a). Synchrony and perturbation transmission in trophic metacommunities. *The American Naturalist*, **714131**. https://doi.org/10.1086/714131

Quévreux, P., Pigeault, R., & Loreau, M. (2021b). Predator avoidance and foraging for food shape synchrony and response to perturbations in trophic metacommunities. *Journal of Theoretical Biology*, **528**, 110836. https://doi.org/10.1016/j.jtbi.2021.110836

Rebolo-Ifrán, N., Tella, J. L., & Carrete, M. (2017). Urban conservation hotspots: Predation release allows the grassland-specialist burrowing owl to perform better in the city. *Scientific Reports*, **7**(1), 3527. https://doi.org/10.1038/s41598-017-03853-z

Resetarits, E. J., Cathey, S. E., & Leibold, M. A. (2018). Testing the keystone community concept: Effects of landscape, patch removal, and environment on metacommunity structure. *Ecology*, **99**(1), 57–67. https://doi.org/10.1002/ecy.2041

Rooney, N., & McCann, K. S. (2012). Integrating food web diversity, structure and stability. *Trends in Ecology & Evolution*, **27**(1), 40–46. https://doi.org/10.1016/j.tree.2011.09.001

Rooney, N., McCann, K. S., Gellner, G., & Moore, J. C. (2006). Structural asymmetry and the stability of diverse food webs. *Nature*, **442**(7100), 265–269. https://doi.org/10.1038/nature04887

Ruokolainen, L., Abrams, P. A., McCann, K. S., & Shuter, B. J. (2011). The roles of spatial heterogeneity and
adaptive movement in stabilizing (or destabilizing) simple metacommunities. *Journal of Theoretical Biology*, **291**, 76-87. https://doi.org/10.1016/j.jtbi.2011.09.004

Schiesari, L., Matias, M. G., Prado, P. I., Leibold, M. A., Albert, C. H., Howeth, J. G., Leroux, S. J., Pardini, R., Siqueira, T., Brancalion, P. H., Cabeza, M., Coutinho, R. M., Diniz-Filho, J. A. F., Fournier, B., Lahr, D. J., Lewinsohn, T. M., Martins, A., Morsello, C., Peres-Neto, P. R., ... Vázquez, D. P. (2019). Towards an applied metaecology. *Perspectives in Ecology and Conservation*, **17**(4), 172-181.
https://doi.org/10.1016/j.pecon.2019.11.001

Schindler, D. E., & Scheuerell, M. D. (2002). Habitat coupling in lake ecosystems. *Oikos*, *98*(2), 177-189.
https://doi.org/10.1034/j.1600-0706.2002.980201.x

Schmitz, O. J. (2004). Perturbation and abrupt shift in trophic control of biodiversity and productivity: Perturbation and regime shift. *Ecology Letters*, **7**(5), 403-409.
https://doi.org/10.1111/j.1461-0248.2004.00592.x

Schmitz, O. J., Hawlena, D., & Trussell, G. C. (2010). Predator control of ecosystem nutrient dynamics. *Ecology Letters*, **13**(10), 1199-1209.
https://doi.org/10.1111/j.1461-0248.2010.01511.x

Shochat, E., Warren, P., Faeth, S., Mcintyre, N., & Hope, D. (2006). From patterns to emerging processes in mechanistic urban ecology. *Trends in Ecology & Evolution*, **21**(4), 186-191.
https://doi.org/10.1016/j.tree.2005.11.019

Shochat, E., Lerman, S. B., Anderies, J. M., Warren, P. S., Faeth, S. H., & Nilon, C. H. (2010). Invasion, competition, and biodiversity loss in urban ecosystems. *BioScience*, **60**(3), 199-208.
https://doi.org/10.1525/bio.2010.60.3.6

Soulé, M. E., Estes, J. A., Miller, B., & Honnold, D. L. (2005). Strongly interacting species: Conservation policy, management, and ethics. *BioScience*, **55**(2), 168.
https://doi.org/10.1641/0006-3568(2005)055[0168:SISCPM]2.0.CO;2

Steele, J. H. (1974). Spatial heterogeneity and population stability. *Nature*, **248**(5443), 83-83.
https://doi.org/10.1038/248083a0

Thaker, M., Vanak, A. T., Owen, C. R., Ogden, M. B., Niemann, S. M., & Slotow, R. (2011). Minimizing predation risk in a landscape of multiple predators: Effects on the spatial distribution of African ungulates. *Ecology*, **92**(2), 398-407.
https://doi.org/10.1890/10-0126.1

Vadeboncoeur, Y., McCann, K. S., Zanden, M. J. V., & Rasmussen, J. B. (2005). Effects of multi-chain omnivory on the strength of trophic control in lakes. *Ecosystems*, **8**(6), 682-693.
https://doi.org/10.1007/s10021-003-0149-5

Van Teeffelen, A. J., Vos, C. C., & Opdam, P. (2012). Species in a dynamic world: Consequences of habitat network dynamics on conservation planning. *Biological Conservation*, **153**, 239-253.
https://doi.org/10.1016/j.biocon.2012.05.001

Vasseur, D. A., & Fox, J. W. (2009). Phase-locking and environmental fluctuations generate synchrony in a predator-prey community. *Nature*, **460**(7258), 1007-1010.
https://doi.org/10.1038/nature08208

Wang, S., Haegeman, B., & Loreau, M. (2015). Dispersal and metapopulation stability. *PeerJ*, **3**, e1295.
https://doi.org/10.7717/peerj.1295

Wang, S., Lamy, T., Hallett, L. M., & Loreau, M. (2019). Stability and synchrony across ecological hierarchies in heterogeneous metacommunities: Linking theory to data. *Ecography*, **0**(0).
https://doi.org/10.1111/ecog.04290

Wang, S., & Loreau, M. (2014). Ecosystem stability in space: \textalpha, \textbeta and \textgamma variability. *Ecology Letters*, **17**(8), 891-901.
https://doi.org/10.1111/ele.12292

Wang, S., & Loreau, M. (2016). Biodiversity and ecosystem stability across scales in metacommunities (F. Jordan, Ed.). *Ecology Letters*, **19**(5), 510-518.
https://doi.org/10.1111/ele.12582

Wilcox, K. R., Tredennick, A. T., Koerner, S. E., Grman, E., Hallett, L. M., Avolio, M. L., Pierre, K. J. L., Houseman, G. R., Isbell, F., Johnson, D. S., Alatalo, J. M., Baldwin, A. H., Bork, E. W., Boughton, E. H., Bowman, W. D., Britton, A. J., Cahill, J. F., Collins, S. L., Du, G., ... Zhang, Y. (2017). Asynchrony among local communities stabilises ecosystem function of metacommunities. *Ecology Letters*, **20**(12), 1534-1545.
https://doi.org/10.1111/ele.12861
# S1 Complementary material and methods

## S1-1 Model description

The model has been originally developed by Barbier and Loreau (2019), who considered a food chain model with a simple metabolic parametrisation. Their model corresponds to the "intra-patch dynamics" part of equations (5a) and (5b) to which we graft a dispersal term to consider a metacommunity with two patches.

$$ \frac{dB_1^{(1)}}{dt} = B_1^{(1)}(\omega_1 g_1 - D_1 B_1^{(1)} - \gamma_1 \alpha_{2,1} B_2^{(1)}) + \delta_1(B_1^{(2)} - B_1^{(1)}) \quad (5a) $$

$$ \frac{dB_i^{(1)}}{dt} = \underbrace{B_i^{(1)}(-r_i - D_i B_i^{(1)} + \gamma_1 \epsilon \alpha_{i,i-1} B_{i-1}^{(1)} - \gamma_1 \alpha_{i+1,i} B_{i+1}^{(1)})}_{\text{intra-patch dynamics}} + \underbrace{\delta_i(B_i^{(2)} - B_i^{(1)})}_{\text{dispersal}} \quad (5b) $$

$B_i^{(1)}$ is the biomass of trophic level i in the patch #1, $\epsilon$ is the biomass conversion efficiency and $\alpha_{i,j}$ is the interaction strength between consumer i and prey j. Species i disperses between the two patches at rate $\delta_i$. The density independent net growth rate of primary producers $g_i$ in equations (5a), the mortality rate of consumers $r_i$ in equations (5b) and the density dependent mortality rate $D_i$ scale with species metabolic rates $m_i$ as biological rates are linked to energy expenditure.

$$ g_1 = m_1 g \quad r_i = m_i r \quad D_i = m_i D \quad (6) $$

In order to get a broad range of possible responses, we assume the predator-prey metabolic rate ratio $m$ and the interaction strength to self-regulation ratio $a$ to be constant. These ratios capture the relations between parameters and trophic levels. This enables us to consider contrasting situations while keeping the model as simple as possible.

$$ m = \frac{m_{i+1}}{m_i} \quad a = \frac{\alpha_{i,i-1}}{D_i} \quad d_i = \frac{\delta_i}{D_i} \quad (7) $$

<figure><img src="image_10.png" /></figure>

Varying $m$ leads to food chains where predators have faster or slower biomass dynamics than their prey and varying $a$ leads to food chains where interspecific interactions prevail or not compared with intraspecific interactions. As all biological rates are rescaled by $D_i$, we also define $d_i$, the dispersal rate relative to self-regulation (referred as scaled dispersal rate in the rest of the study), in order to keep the values of the dispersal rate relative to the other biological rates consistent across trophic levels. Finally, the time scale of the system is defined by setting the metabolic rate of the primary producer $m_1$ to unity. Thus, we can transform equations (5a) and (5b) into:

$$ \frac{1}{D} \frac{dB_1^{(1)}}{dt} = B_1^{(1)}(\omega \frac{g}{D} - B_1^{(1)} - \gamma m a B_2^{(1)}) + d_1(B_1^{(2)} - B_1^{(1)}) \quad (8a) $$

$$ \frac{1}{m^{i-1}D} \frac{dB_i^{(1)}}{dt} = \underbrace{B_i^{(1)}(-\frac{r}{D} - B_i^{(1)} + \gamma \epsilon a B_{i-1}^{(1)} - \gamma m a B_{i+1}^{(1)})}_{\text{intra-patch dynamics}} + \underbrace{d_i(B_i^{(2)} - B_i^{(1)})}_{\text{dispersal}} \quad (8b) $$

Thus, $\epsilon a$ and $m a$ defines the positive effect of the prey on its predator and the negative effect of the predator on its prey, respectively. These two synthetic parameters define the overall behaviour of the food chain and will be varied over the interval [0.1, 10] to consider a broad range of possible responses. Finally, the mortality rate is set to zero ($r = 0$) to remove the energetic limitations of the food chain and make interactions the dominant factors determining biomass distribution and stability patterns, as in Barbier and Loreau (2019).

## S1-2 Biomass at equilibrium when top predators populations are perfectly coupled

The system can be easily solved if we consider the total population of top predator instead of two populations connected by dispersal. Since the two populations are perfectly coupled by dispersal, top predator i biomass is constant across patches and we have $B_i^{(1)*} = B_i^{(2)*} = 0.5 B_i^{tot*}$. Thus we have the
The elements of $\vec{E}$ are defined by stochastic perturbations $E_i = \sigma_i dW_i$ with $\sigma_i$ their standard deviation and $dW_i$ a white noise term with mean 0 and variance 1. In our model, each species $i$ in each patch $k$ can receive demographic perturbations scaling with the square root of their biomass at equilibrium. Thus, $\vec{E}$ contains the white noise term $\sigma_i^{(k)} dW_i^{(k)}$ for each population of each species, $T$ is a diagonal matrix whose terms are $\sqrt{B_i^{* (k)}}$ and the matrix product $T \vec{E}$ results in the product of the white noise and the biomass scaling as in equation (2) in the main text.

## S1-5 Demonstration of the Lyapunov equation

The following demonstration of the Lyapunov equation has been taken from Oku and Aihara (2018). The continuous-time dynamics from equation (15) can be converted to a discrete-time dynamics by using Euler-Maruyama method:

$$ \vec{X}_{t+\Delta t} = \vec{X}_t + \Delta t J \vec{X}_t + \sqrt{\Delta t} T \vec{E}_t \quad (16) $$

$C^* = \mathbb{E}[\vec{X}_t \vec{X}_t^\top]$ (the expected value of the product $\vec{X}_t$ and its transpose $\vec{X}_t^\top$) is the stationary variance-covariance matrix of the system, therefore $dC^*/dt = 0$. We also have the following relation:

$$ \begin{aligned} \frac{dC^*}{dt} &= \lim_{\Delta t \to 0} \frac{\mathbb{E}[\vec{X}_{t+\Delta t} \vec{X}_{t+\Delta t}^\top] - \mathbb{E}[\vec{X}_t \vec{X}_t^\top]}{\Delta t} \\ &= \lim_{\Delta t \to 0} \frac{\mathbb{E}[(\vec{X}_t + \Delta t J \vec{X}_t + \sqrt{\Delta t} T \vec{E}_t)(\vec{X}_t + \Delta t J \vec{X}_t + \sqrt{\Delta t} T \vec{E}_t)^\top] - \mathbb{E}[\vec{X}_t \vec{X}_t^\top]}{\Delta t} \\ &= \lim_{\Delta t \to 0} \frac{\Delta t \mathbb{E}[\vec{X}_t \vec{X}_t^\top J^\top] + \Delta t \mathbb{E}[J \vec{X}_t \vec{X}_t^\top] + \Delta t^2 \mathbb{E}[J \vec{X}_t \vec{X}_t^\top J^\top] + \Delta t \mathbb{E}(T \vec{E}_t \vec{E}_t^\top T^\top)}{\Delta t} \\ &= \mathbb{E}[\vec{X}_t \vec{X}_t^\top J^\top] + \mathbb{E}[J \vec{X}_t \vec{X}_t^\top] + \mathbb{E}[T \vec{E}_t \vec{E}_t^\top T^\top] \\ &= C^* J^\top + JC^* + TV_E T^\top = 0 \end{aligned} \quad (17) $$

Because $\mathbb{E}[\vec{X}_t] = 0$, $\mathbb{E}[\vec{E}_t] = 0$, $\mathbb{E}[\vec{X}_t \vec{E}_t^\top] = 0$, $\mathbb{E}[\vec{E}_t \vec{X}_t^\top] = 0$ and $V_E = \mathbb{E}[\vec{E}_t \vec{E}_t^\top]$ the variance-covariance matrix of stochastic perturbations.

## S1-6 Resolution of the Lyapunov equation

In the vicinity of equilibrium, the Lyapunov equation links the variance-covariance matrix $V_E$ of the perturbation vector $\vec{E}$ to the variance-covariance matrix $C^*$ of species biomasses (see the appendix of Wang et al. (2015) for more details on the Lyapunov equation).

$$ JC^* + C^* J^\top + TV_E T^\top = 0 \quad (18) $$

The diagonal elements of $V_E$ are equal to $\sigma_i^2$ (variance of the white noises) and the non-diagonal elements are equal to zero because perturbations are independent. $\top$ is the transpose operator. $C^*$ can be calculated using a Kronecker product (Nip et al., 2013). The Kronecker product of an $m \times n$ matrix $A$ and a $p \times q$ matrix $B$ denoted $A \otimes B$ is the $mp \times nq$ block matrix given by:

$$ A \otimes B = \begin{pmatrix} a_{11}B & \cdots & a_{1n}B \\ \vdots & \ddots & \vdots \\ a_{m1}B & \cdots & a_{mn}B \end{pmatrix} $$

We define $C_s^*$ and $(TV_E T^\top)_s$ the vectors stacking the columns of $C^*$ and $TV_E T^\top$ respectively. Thus, equation (18) can be rewrite as:

$$ \begin{aligned} (J \otimes I + I \otimes J)C_s^* &= -(TV_E T^\top)_s \\ C_s^* &= -(J \otimes I + I \otimes J)^{-1}(TV_E T^\top)_s \end{aligned} \quad (19) $$

## S1-7 Coefficient of variation and correlation

Our different metrics of stability can be easily computed from the elements of the variance-covariance matrix $C^*$ defined by elements $w_{i(k)j(l)}$ that are the covariance between species $i$ in patch $k$ and species $l$.
j in patch $ℓ$.

$$C^* = \begin{pmatrix} w_{1(1)1(1)} & \cdots & w_{1(1)S(1)} & \cdots & w_{1(n)1(n)} & \cdots & w_{1(n)S(n)} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ w_{S(1)1(1)} & \cdots & w_{S(1)S(1)} & \cdots & w_{S(n)1(n)} & \cdots & w_{S(n)S(n)} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ w_{1(n)1(1)} & \cdots & w_{1(n)S(1)} & \cdots & w_{1(n)1(n)} & \cdots & w_{1(n)S(n)} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ w_{S(n)1(1)} & \cdots & w_{S(n)S(1)} & \cdots & w_{S(n)1(n)} & \cdots & w_{S(n)S(n)} \end{pmatrix} \quad (20)$$

The temporal variability of the metacommunity is assessed with the coefficient of variation (CV) of biomass at different scales: **population scale CV<sub>i</sub><sup>(k)</sup>**, which is the biomass CV of species *i* in patch *k*, **metapopulation scale CV<sub>i</sub>**, which is the biomass CV of the total biomass of species *i* across patches and **metacommunity scale CV<sub>MC</sub>**, which is the total biomass of the entire metacommunity (Wang and Loreau, 2014; Wang et al., 2019; Jarillo et al., 2022).

$$CV_i^{(k)} = \frac{\sqrt{w_{i(k)1(k)}}}{\mu_i^{(k)}} \qquad CV_i = \frac{\sqrt{\sum_{k\ell} w_{i(k)j(\ell)}}}{\sum_{k} B_i^{*(k)}} \qquad CV_{MC} = \frac{\sqrt{\sum_{ijk\ell} w_{i(k)j(\ell)}}}{\sum_{ik} B_i^{*(k)}} \quad (21)$$

<figure><img src="image_11.png" /><figcaption>Figure S1-1: Elements of the variance-covariance matrix $C^*$ used to compute the biomass CV at different scales defined in equation (21).</figcaption></figure>

The correlation matrix $R^*$ of the system, whose elements $\rho_{i(k)j(\ell)}$ are defined by:

$$\rho_{i(k)j(\ell)} = \frac{w_{i(k)j(\ell)}}{\sqrt{w_{i(k)1(k)}} \sqrt{w_{j(\ell)1(\ell)}}} \quad (22)$$

## S1-8 Asymptotic resilience

In addition to the response to stochastic perturbations, we consider asymptotic resilience to measure the long term return time of the metacommunity. Asymptotic resilience is measured by the opposite of the real part of the dominant eigenvalue $\lambda_{dom}$ of Jacobian matrix $J$ ($-\Re(\lambda_{dom})$). Since the dominant eigenvalue is the eigenvalue with the largest real part and we only consider ecosystems at equilibrium (i.e. all eigenvalues have negative real parts), the lower the real part of the dominant eigenvalue, the faster the long term return time.

Moreover, we can assess the influence of each species on asymptotic resilience by comparing the absolute value of the real part of each element $e_i$ of the dominant eigenvector ($E_{dom}$). Because $e_i$ is the contribution of species *i* to $E_{dom}$, $|e_i|/\sum_{j=1}^n |e_j|$ is the relative weight of species *i* in the dynamics of long term return to equilibrium (with *n* the number of populations in the metacommunity).

## References

* Arnoldi, J.-F., Loreau, M., & Haegeman, B. (2016). Resilience, reactivity and variability: A mathematical comparison of ecological stability measures. *Journal of Theoretical Biology*, 389, 47–59. https://doi.org/10.1016/j.jtbi.2015.10.012
Barbier, M., & Loreau, M. (2019). Pyramids and cascades: A synthesis of food chain functioning and stability. *Ecology Letters*, 22(2), 405-419. https://doi.org/10.1111/ele.13196

Jarillo, J., Cao-García, F. J., & De Laender, F. (2022). Spatial and ecological scaling of stability in spatial community networks. *Frontiers in Ecology and Evolution, 10*, 861537. https://doi.org/10.3389/fevo.2022.861537

Nip, M., Hespanha, J. P., & Khammash, M. (2013). Direct numerical solution of algebraic Lyapunov equations for large-scale systems using Quantized Tensor Trains. *52nd IEEE Conference on Decision and Control*, 1950–1957. https://doi.org/10.1109/CDC.2013.6760167

Oku, M., & Aihara, K. (2018). On the covariance matrix of the stationary distribution of a noisy dynamical system. *Nonlinear Theory and Its Applications, IEICE*, 9(2), 166-184. https://doi.org/10.1587/nolta.9.166

Wang, S., Haegeman, B., & Loreau, M. (2015). Dispersal and metapopulation stability. *PeerJ*, 3, e1295. https://doi.org/10.7717/peerj.1295

Wang, S., Lamy, T., Hallett, L. M., & Loreau, M. (2019). Stability and synchrony across ecological hierarchies in heterogeneous metacommunities: Linking theory to data. *Ecography*, 0(0). https://doi.org/10.1111/ecog.04290

Wang, S., & Loreau, M. (2014). Ecosystem stability in space: \texttt{alpha}, \texttt{beta} and \texttt{gamma} variability. *Ecology Letters*, 17(8), 891-901. https://doi.org/10.1111/ele.12292
# S2 Complementary results

## S2-1 General description of parameters

<figure><img src="image_12.png" /><figcaption>Figure S2-1: Distribution of parameters and their effects on an isolated food chain. A) Biomass distribution depending on the positive effect of prey on predator $\epsilon a$ and the negative effect of predator on prey $ma$. B) Value of the ratio of predator to prey metabolic rate $m = m_{i+1}/m_i$ for each combination of $\epsilon a$ and $ma$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure S2-2: Distribution of biomass in an isolated predator-prey system (without dispersal) depending on interaction strength $a$ relative to self-regulation. Increasing the asymmetry of the interaction strength $\gamma$ is equivalent to increasing $a$ ($m = 0.65$).</figcaption></figure>

Increasing the interaction strength $a$ relative to self-regulation decreases the biomass of prey because of the increased mortality due to predation (Figure S2-2). However, the biomass of predators follows a hump-shaped relationship: it first increases due to the increased resource consumption and then decreases because of prey overexploitation.

The effect of perturbations of populations within a community of $S$ species can be assessed by the ratio of the mean variance of species biomass $B_j$ to the variance of perturbations $\sigma_k$:

$$ \frac{\frac{1}{S} \sum_{j}^{S} \mathrm{Var}(B_j)}{\frac{1}{S} \sum_{k}^{S} \sigma_k^2} = \frac{\frac{1}{S} \sum_{j}^{S} \mathrm{Var}(B_j)}{\sigma_i^2} \quad \text{because we only consider one perturbation effecting species } i \quad (23) $$

As demonstrated by Arnoldi et al. (2019), exogenous perturbations affect more rare species, demographic perturbations evenly affect species regardless on the biomass distribution and environmental perturbations
affect more abundant populations. Therefore, we consider demographic perturbations to perturb the entire community with the same intensity regardless on the biomass variations caused by varying $\gamma$ (Figure S2-3).

<figure><img src="image_14.png" /><figcaption>Figure S2-3: Ratio of the mean variance of species biomass to the mean variance of environmental perturbations (see equation (23)) depending on the biomass of the perturbed species. Three types of perturbations with different scaling with the equilibrium biomass of the perturbed species $i$ ($B_i^{*z}$) are tested: exogenous perturbations ($z = 0$), demographic perturbations ($z = 0.5$) and environmental perturbations ($z = 1$).</figcaption></figure>

## S2-2 Dispersal of predators and perturbation of prey

### S2-2-1 Source-sink effect

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure S2-4: A) Distribution of the biomass of each species among patches depending on the asymmetry of interaction strength $\gamma$. B) Distribution of biomasses scaled by their value in a metacommunity without dispersal ($B_{scaled} = B_{d_2>0}/B_{d_2=0}$).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure S2-5: The asymmetry of interaction strength alters the biomass distribution between the two patches. Increased the interaction strength in patches #1 enhances the biomass production of predators that spill over patches from patch #1 to patch #2, therefore increasing prey biomass in patch #1 and decreasing it in patch #2 ($\gamma = 3$). Plain rectangles represent the biomass of each population, while dashed rectangles represent the same population in a metacommunity without dispersal (no spillover effect)</figcaption></figure>

Varying the asymmetry of interaction strength $\gamma$ is equivalent to varying the interaction strength $a$ in patch #1 and has the same effects on species biomass: first increasing $\gamma$ increases predator biomass by increasing prey consumption, then it decreases predator biomass because of resource overexploitation (Figure S2-2). This leads to different biomass distributions in patches #1 and #2 (Figure S2-4A). Predator biomass increases with $\gamma$ and is the same in both patches because their high dispersal rate balances any difference. Prey biomass is higher in patch #1 than in patch #2, and both decrease with $\gamma$, except in patch #1, where we first observe a small increase for $\gamma < 2$. This response is due to source-sink effects: the increase in prey consumption in patch #1 increases predator biomass (source) that spills over patch #2 (sink) due to dispersal (Figure S2-4B and Figure S2-5). Therefore, predator biomass is lower in patch #1 and higher in patch #2 compared to what we expect in the same food chains in isolation (i.e., without dispersal). This also prevents predators from overexploiting prey in patch #1 by spreading the increased predator biomass across the metacommunity, which explains why we do not observe a decrease in predator biomass for high values of $\gamma$, as shown in Figure S2-2. Conversely, the distribution of prey biomass across the two patches is opposite (higher in patch #1 and lower in patch #2).

### S2-2-2 Conditions of coexistence

Asymmetry and dispersal lead to competition, apparent competition and source-sink dynamics that can rescue or drive local populations to extinction. Therefore, we consider limit cases in which dispersal is infinite (well mixed populations across the metacommunity) to analytically calculate biomasses at equilibrium and determine the range of values of $\omega$ and $\gamma$ enabling the coexistence of all populations of each species.

We consider the total biomass of predators $B_2^{tot} = B_2^{(1)} + B_2^{(2)}$ and because the very high dispersal of predators equally distributes its biomass among the two patches, we have $B_2^{tot} = 2 \times B_2^{(1)}$. Then, we can define the system:

$$ \frac{dB_1^{(1)}}{dt} = DB_1^{(1)} \left( \frac{\omega g}{D} - B_1^{(1)} - \gamma ma \frac{B_2^{tot}}{2} \right) \qquad (24a) $$

$$ \frac{dB_1^{(2)}}{dt} = DB_1^{(2)} \left( \frac{g}{D} - B_1^{(2)} - ma \frac{B_2^{tot}}{2} \right) \qquad (24b) $$

$$ \frac{dB_2^{tot}}{dt} = mD \frac{B_2^{tot}}{2} \left( -\frac{r}{D} - B_2^{tot} + \varepsilon a(\gamma B_1^{(1)} + B_1^{(2)}) \right) \qquad (24c) $$

Since $r = 0$, we remove it from the equations for the sake of simplicity. We define $\lambda = \varepsilon ma^2$, which is
the intensity of top-down control defined by Barbier and Loreau (2019). At equilibrium, we obtain:

$$B_1^{(1)*} = \frac{g}{D} \left( \frac{2\omega + \omega\lambda - \gamma\lambda}{2 + \lambda(\gamma^2 + 1)} \right) \qquad (25a)$$

$$B_1^{(2)*} = \frac{g}{D} \left( \frac{\lambda\gamma^2 - \omega\lambda\gamma + 2}{2 + \lambda(\gamma^2 + 1)} \right) \qquad (25b)$$

$$B_2^{tot*} = \frac{2\varepsilon ag(1 + \omega\gamma)}{D(2 + \varepsilon a^2 m(\gamma^2 + 1))} \qquad (25c)$$

Prey biomass in patch #1 $B_1^{(1)}$ is positive only if:

$$\gamma < \frac{\omega(2 + \lambda)}{\lambda} \xrightarrow{\lambda \to \infty} \omega \qquad (26)$$

Prey biomass in patch #2 $B_1^{(2)}$ is positive if $f(\gamma) = \lambda\gamma^2 - \omega\lambda\gamma + 2 > 0$. $f$ opens upwards: thus, if $\omega < \sqrt{8/\lambda}$, $f$ has no roots and is always positive. Otherwise, $B_1^{(2)}$ is positive if:

$$\gamma > \frac{\lambda\omega + \sqrt{\lambda(\lambda\omega^2 - 8)}}{2\lambda} \xrightarrow{\lambda \to \infty} \omega \quad \text{if} \quad \omega > \sqrt{\frac{8}{\lambda}} \quad \text{or} \qquad (27a)$$

$$\gamma < \frac{\lambda\omega - \sqrt{\lambda(\lambda\omega^2 - 8)}}{2\lambda} \xrightarrow{\lambda \to \infty} 0 \quad \text{if} \quad \omega > \sqrt{\frac{8}{\lambda}} \qquad (27b)$$

Predators $B_2$ thrive in each patch for all values of $\omega$ and $\gamma$ (Figure S2-6C). Hence, coexistence is ensured for all values of top-down control $\lambda$, asymmetry of resource supply $\omega$ and asymmetry of interaction strength $\gamma$ only if $\gamma = \omega$. In the main text, we always consider $\gamma = \omega$ (Figure S2-6D), but their independent effects are detailed in the following.
<figure><img src="image_18.png" /><figcaption>Figure S2-6: Distribution of parameters, asymmetry of resource supply $\omega$ and asymmetry of interaction strength $\gamma$, leading to the coexistence of predator and prey in each patch. Only predators are able to disperse at an infinite rate (well-mixed predator populations). This distribution is assessed for different values of the positive effect of prey on predator $\epsilon a$ and negative effect of predator on prey $ma$. A) Biomass of prey in patch #1 $B_1^{*(1)}$ and B) in patch #2 $B_1^{*(2)}$. C) Biomass of predator in patches #1 and #2 ($B_2^{*(1)} = B_2^{*(2)}$ because predator populations are well mixed). D) Coexistence of predator and prey in each patch.</figcaption></figure>

### S2-2-3 Nontransitivity of correlation

To explain the correlation between prey populations, we can track the transmission of perturbations in the metacommunity. Increasing the asymmetry of interaction strength $\gamma$ tends to decorrelate predator and prey dynamics within each patch (Figure S2-7A). When prey are perturbed in patch #1, the dynamics of predator and prey biomass are correlated in patch #1 and anticorrelated in patch #2 due to the bottom-up and top-down transmissions of perturbations, respectively. Although we would expect the two populations of prey to be anti-correlated according to the mechanism described by Quévreux et al. (2021) (see Figure S2-28 in the following), we actually observe a weak correlation of these two populations (Figure S2-7B). In the same way, the intermediate correlation and anti-correlation of predator and prey when prey are perturbed in patch #2 do not explain the strong anti-correlation of prey populations (Figure S2-7C). Therefore, other mechanisms are acting in our system.
<figure><img src="image_19.png" /><figcaption>Figure S2-7: Correlation between predator and prey with each patch ($εa = 1$, $ma = 1$, $d_2 = 10^6$, $ω = γ$). A) Correlation in each patch depending on asymmetry of interaction strength $γ$ when prey in patch #1 (left panel) or in patch #2 (right panel) are perturbed. Labels and vertical dashed lines represent the correlation values used in panels B and C. B) Schematic representation of correlations when $γ = 7$. Coloured double arrows and their associated number (see also the colour scale) represent the correlations between populations. Prey are perturbed in patch #1. C) Prey are perturbed in patch #2.</figcaption></figure>

### S2-2-4 Complete effects of $εa$ and $ma$

The stability patterns observed in Figures 2, 3 and 5 in the main text are also observed for a wide range of ecological and physiological parameters aggregated into the positive effect of prey on predators $εa$ and the negative effect on predators on prey $ma$. Therefore, our results are robust and the identified mechanisms are specific to a particular combination of parameters.

The response of the asymptotic resilience to the asymmetry of interaction strength (Figure S2-11A) is not similar to the results of Rooney et al. (2006). Indeed, we do not observe minimum of resilience for $γ = 1$ for all combinations of $εa$ and $ma$. The variations in asymptotic resilience depend on the relative contribution of each population of each species (Figure S2-11B), which is governed by the biomass distribution of each species among patches (Figure S2-8A) and the ratio of predator to prey metabolic rate ratio $m$ (Figure S2-1B). As demonstrated by Haegeman et al. (2016) and Arnoldi et al. (2018), rare species control the long term response to perturbations of the metacommunity (i.e., the asymptotic resilience), as well as species with a slow pace of life (i.e., a slow metabolism).

<figure><img src="image_20.png" /><figcaption>Figure S2-8: Biomass distribution of each species in each patch depending on asymmetry of interaction strength $γ$, positive effect of prey on predator $εa$ and negative effect of predator on prey $ma$. A) Biomass distribution. B) Biomass scaled by the biomass in the metacommunity without dispersal ($d_2 = 0$).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure S2-9: Correlation between populations depending on asymmetry in interaction strength $\gamma$ when predators disperse and prey are perturbed in patch #1 or #2. Predators have a high scaled dispersal rate ($d_2 = 10^6$), which strongly couples their two populations ($\gamma = \omega$).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure S2-10: Biomass CV at different scales depending on asymmetry of interaction strength $\gamma$, positive effect of prey on predator $\varepsilon a$ and negative effect of predator on prey $ma$ ($d_2 = 10^6$ and $\omega = \gamma$). A) Biomass CV of the population of each species in each patch. B) CV of the total biomass of each species.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure S2-11: Linear stability depending on asymmetry of interaction strength $\gamma$, positive effect of prey on predator $\varepsilon a$ and negative effect of predator on prey $ma$. **A)** Asymptotic resilience (real part of the dominant eigenvalue of the Jacobian matrix) **B)** Contribution of the populations of each species to the dominant eigenvector.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure S2-12: Time series of biomasses rescaled by their value at equilibrium after an increase in prey biomass by 20% in patch #1 (left panel) or patch #2 (right panel) for two values of interaction strength asymmetry ($\gamma = 1$ or $\gamma = 10$, $\varepsilon a = 1$, $ma = 1$, $d_2 = 10^6$ and $\omega = \gamma$).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure S2-13: Time series of biomasses rescaled by their value at equilibrium after an increase in prey biomass by 20% in patch #1 (left panel) or patch #2 (right panel) depending on asymmetry of interaction strength $\gamma$, positive effect of prey on predator $\epsilon a$ and negative effect of predator on prey $ma$ ($\gamma = 3$, $d_2 = 10^6$ and $\omega = \gamma$).</figcaption></figure>

### S2-2-5 Effect of asymmetry of resource supply $\omega$

According to Figure S2-6, we set $\omega = \gamma$ to ensure the coexistence of prey and predators in each patch for all combinations of $\epsilon a$ and $ma$. Varying the asymmetry of resource supply $\omega$ does not qualitatively alter the response of biomass (Figure S2-14), correlation (Figure S2-15A) and biomass CV (Figure S2-15B) to the variations in the asymmetry of interaction strength $\gamma$.

<figure><img src="image_26.png" /><figcaption>Figure S2-14: Biomass distribution of each species in each patch depending on asymmetry of interaction strength $\gamma$ and biomass production $\omega$ ($\epsilon a = 0.1$, $ma = 0.1$ and $d_2 = 10^6$). A) Biomass distribution. B) Biomass scaled by the biomass in the metacommunity without dispersal ($d_2 = 0$).</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure S2-15: Stability depending on asymmetry of interaction strength $\gamma$ and biomass production $\omega$ ($\epsilon a = 0.1$, $ma = 0.1$ and $d_2 = 10^6$). **A)** Correlation between populations. **B)** Biomass CV of the population of each species in each patch.</figcaption></figure>

### S2-2-6 Effect of perturbation of predators

The perturbation of predators leads to the same response regardless of the perturbed patch because the very high dispersal of predators perfectly synchronises their population dynamics. The asymmetry of interaction strength leads to different dynamics in each patch that decreases the correlation of prey dynamics (Figure S2-16) and stabilises predator dynamics by decreasing their biomass CV (Figure S2-17A).

<figure><img src="image_28.png" /><figcaption>Figure S2-16: Correlation between populations depending on asymmetry of interaction strength $\gamma$, positive effect of prey on predator $\epsilon a$ and negative effect of predator on prey $ma$. Predators disperse and are perturbed in patch #1 or patch #2 ($d_2 = 10^6$ and $\omega = \gamma$).</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Figure S2-17: Biomass CV at different scales depending on asymmetry of interaction strength $\gamma$, positive effect of prey on predator $\varepsilon a$ and negative effect of predator on prey $ma$. Predators disperse and are perturbed in patch #1 or patch #2 ($d_2 = 10^6$ and $\omega = \gamma$). A) CV of the total biomass of each species. B) CV of the total biomass of the metacommunity.</figcaption></figure>

### S2-2-7 Effect of food chain length

Here, we consider three trophic levels to extend our results to metacommunities with longer food chain lengths. In this setup, only top predators (species 3) are able to disperse, and basal species (species 1) receive stochastic perturbations. $\gamma$ also has the same value across trophic levels. We observe a similar response to the case with two trophic levels for the correlations of the dynamics of the biomass of species 2 and 3 (Figure S2-18) as well as for biomass CV (Figure S2-19). However, the response of species 1 in completely different. Therefore, the mechanisms described in the main text are only acting for the dispersing species and the species directly interacting with it.

<figure><img src="image_30.png" /><figcaption>Figure S2-18: Correlation between populations in a three trophic level food chain depending on the asymmetry of interaction strength $\gamma$, positive effect of prey on predator $\varepsilon a$ and negative effect of predator on prey $ma$ ($d_3 = 10^6$ and $\omega = 1$). Top predators disperse, and the basal species is perturbed in patch #1 or #2.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Figure S2-19: Biomass CV in a three trophic level food chain at different scales depending on asymmetry of interaction strength $\gamma$, positive effect of prey on predator $\varepsilon a$ and negative effect of predator on prey $ma$ ($d_3 = 10^6$ and $\omega = 1$). **A)** Biomass CV of the population of each species in each patch. **B)** CV of the total biomass of each species.</figcaption></figure>

## S2-3 Dispersal of prey and perturbation of predators

In this section, we consider a setup mirroring the metacommunity model described in the main text. Here, only prey are able to disperse at a very high rate ($d_1 = 10^6$), and predators receive stochastic perturbations. We also set $\gamma = \omega$ to be consistent with the results in the main text. In the following, we find the same responses to the asymmetry of interaction strength $\gamma$, which demonstrates that the mechanisms described in the main text are not conditioned by the trophic position of the dispersing species.

### S2-3-1 Conditions of coexistence

When only prey are able to disperse, all populations of each species have positive biomasses for all values of $\omega$, $\gamma$, $\varepsilon a$ and $ma$ (Figure S2-20).
<figure><img src="image_32.png" /><figcaption>Figure S2-20: Distribution of parameters, asymmetry of resource supply $\omega$ and asymmetry of interaction strength $\gamma$, leading to the coexistence of predator and prey in each patch. Only prey are able to disperse at an infinite rate (well mixed prey populations). This distribution is assessed for different values of the positive effect of prey on predator $\varepsilon a$ and the negative effect of predator on prey $ma$. The product of $\varepsilon a$ and $ma$ is the strength of top-down control $\lambda$ ($\lambda = \varepsilon a^2 m$, see Barbier and Loreau (2019)). **A)** Biomass of prey in patches #1 and #2 ($B_1^{*(1)} = B_1^{*(2)}$ because prey populations are well mixed). **B)** Biomass of predator in patch #1 ($B_2^{*(1)}$ **C**) in patch #2 $B_2^{*(2)}$.</figcaption></figure>
S2-3-2 Complete effects of εa and ma

<figure><img src="image_33.png" /><figcaption>Figure S2-21: Biomass distribution of each species in each patch depending on asymmetry of interaction strength $\gamma$, positive effect of prey on predator $\varepsilon a$ and negative effect of predator on prey $ma$. **A)** Biomass distribution. **B)** Biomass scaled by the biomass in the metacommunity without dispersal ($d_1 = 0$). Prey and predator curves perfectly overlap</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure S2-22: Correlation of population when prey disperse and predators are perturbed in patch #1 or in patch #2 depending on asymmetry of interaction strength $\gamma$, positive effect of prey on predator $\varepsilon a$ and negative effect of predator on prey $ma$.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Figure S2-23: Biomass CV at different scales depending on asymmetry of interaction strength $\gamma$, positive effect of prey on predator $\varepsilon a$ and negative effect of predator on prey $ma$ ($d_1 = 10^6$ and $\omega = \gamma$). A) Biomass CV of the population of each species in each patch. B) CV of the total biomass of each species.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure S2-24: Time series after pulse perturbation of predators in patch #1 or in patch #2. Biomasses are scaled by their value at equilibrium, and dispersal is high ($\gamma = 3$ and $d = 10^6$).</figcaption></figure>
S2-4 Correlated environmental perturbations

<figure><img src="image_37.png" /><figcaption>Figure S2-25: Ratio of the mean variance of species biomass to the mean variance of environmental perturbations (see equation (23)) in each patch depending on asymmetry of interaction strength $\gamma$. Each prey population receives spatially correlated environmental perturbation (colour gradient scale) scaling with equilibrium biomass $B_i^*(z = 1)$.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Figure S2-26: Stability at different scales depending on asymmetry of interaction strength $\gamma$ when predators disperse and prey are perturbed in patch #1 or #2 with environmental perturbations ($\epsilon a = 1$, $ma = 1$, $\omega = \gamma$). **A)** Spatial correlation between the populations of each species. **B)** Biomass CV at the population scale. **C)** Biomass CV at the metapopulation scale (CV of the total biomass of each species). **D)** Biomass CV at the metacommunity scale (CV of the total biomass of the metacommunity).</figcaption></figure>
<figure><img src="image_39.png" /><figcaption>Figure S2-27: Stability at different scales depending on asymmetry of interaction strength $\gamma$ when predators disperse and prey are perturbed by a spatially correlated environmental perturbations (colour gradient scale) ($\epsilon a = 1$, $ma = 1$, $\omega = \gamma$). A) Spatial correlation between the populations of each species. B) Biomass CV at the population scale. C) Biomass CV at the metapopulation scale (CV of the total biomass of each species). D) Biomass CV at the metacommunity scale (CV of the total biomass of the metacommunity.</figcaption></figure>

Here, we consider the same metacommunity as in the main text (see Figure 1 in the main text), but prey receive spatially correlated environmental perturbations. Environmental perturbations correspond to the synchronous response of all individuals of the same population to an environmental factor (e.g., drought), and they scale with equilibrium biomass $B_i^*$ (see the supporting information of Quévreux et al. (2021) for the demonstration). In our metacommunity, we also consider that environmental perturbations are spatially correlated since it is reasonable to assume that different populations of the same species will respond in a similar way to environmental perturbations.

The effect of a perturbation on a population within a community can be assessed by the ratio of the mean variance of species biomass $j$ to the variance of the perturbation $i$ by equation (23). As demonstrated by Arnoldi et al. (2019), environmental perturbations affect abundant populations the most, which is the prey population in the fast patch in our case (Figures S2-3 and S2-25). Therefore, we can approximate the effect of environmental perturbations by the effect of the perturbation of prey in the fast patch. Perturbing a single population with demographic or environmental perturbations leads to exactly the same qualitative results (Figure S2-26 and Figure 2 and Figure 3 in the main text), and only the CV values change because of the different biomass scaling.

Increasing the correlation of perturbations increases the correlation of the dynamics of prey populations (Figure S2-27A) because of the Moran effect (Moran, 1953). The increase in synchrony explains the increase in the biomass CV observed at each scale for all species (Figure S2-27B-D), except for prey in the slow patch (Figure S2-27B). The Moran effect is particularly strong at low asymmetry ($\gamma < 4$), but once asymmetry is high enough, two mechanisms disrupt the Moran effect. First, when asymmetry is high, the dynamics in each patch become so different that correlated perturbations are not able to generate similar responses. Second, because of the discrepancy in the distribution of prey biomass among the two patches, environmental perturbations mostly affected prey in the fast patch (Figure S2-25). Therefore, with increasing asymmetry of interaction strength $\gamma$, the response of the metacommunity to correlated environmental perturbations converges towards the response of a metacommunity in which only prey in the fast patch are perturbed.
## S2-5 Reminder of the symmetric case

A

<figure><img src="image_40.png" /></figure>

B

<figure><img src="image_41.png" /></figure>

C

<figure><img src="image_42.png" /></figure>

Figure S2-28: Summary of the main results from Quévreux et al. (2021), who considered a two patch predator-prey metacommunity with passive dispersal. In the setup presented in A), prey are perturbed in patch #1 and only predators are able to disperse. Thus, perturbations have a bottom-up transmission in patch #1 (i.e. transmission from lower to upper trophic levels). This leads to the temporal correlation of the biomass dynamics of predators and prey in patch #1 showed in B)(1) because if a perturbation increases the biomass of prey, it also increases the biomass of predators due to the vertical transfer of biomass. The passive dispersal of predators transmits the perturbations and spatially correlate their populations as shown in B)(2). Then, perturbations have a top-down transmission in patch #2 (i.e. transmission from upper to lower trophic levels). This leads to the temporal anticorrelation (negative coefficient of correlation) of the biomass dynamics of predators and prey in patch #2 showed in B)(3) because if a perturbation increases the biomass of predators, it decreases the biomass of prey due to the negative effect of predators on prey. Eventually, prey populations are spatially anticorrelated, as shown in B)(4). Hence, by knowing which species is perturbed, which species disperses and how perturbations propagate within a food chain, Quévreux et al. (2021) were able to explain the spatial synchrony of the various populations of a metacommunity, summarized by the correlation matrix in C).

## References

* Arnoldi, J.-F., Bideault, A., Loreau, M., & Haegeman, B. (2018). How ecosystems recover from pulse perturbations: A theory of short- to long-term responses. *Journal of Theoretical Biology*, *436*, 79-92. https://doi.org/10.1016/j.jtbi.2017.10.003

* Arnoldi, J.-F., Loreau, M., & Haegeman, B. (2019). The inherent multidimensionality of temporal variability: How common and rare species shape stability patterns (J. Chase, Ed.). *Ecology Letters*, *22*(10), 1557-1567. https://doi.org/10.1111/ele.13345

* Barbier, M., & Loreau, M. (2019). Pyramids and cascades: A synthesis of food chain functioning and stability. *Ecology Letters*, *22*(2), 405-419. https://doi.org/10.1111/ele.13196

* Haegeman, B., Arnoldi, J.-F., Wang, S., de Mazancourt, C., Montoya, J. M., & Loreau, M. (2016). Resilience, invariability, and ecological stability across levels of organization. *bioRxiv*. https://doi.org/10.1101/085852

* Moran, P. (1953). The statistical analysis of the Canadian Lynx cycle. *Australian Journal of Zoology*, *1*(3), 291. https://doi.org/10.1071/ZO9530291

* Quévreux, P., Barbier, M., & Loreau, M. (2021). Synchrony and perturbation transmission in trophic metacommunities. *The American Naturalist*, *714*131. https://doi.org/10.1086/714131

* Rooney, N., McCann, K. S., Gellner, G., & Moore, J. C. (2006). Structural asymmetry and the stability of diverse food webs. *Nature*, *442*(7100), 265-269. https://doi.org/10.1038/nature04887