As mentioned above, stability can be measured by the temporal variability of biomass generated by stochastic perturbations (Arnoldi et al., 2016, 2019). Stochastic perturbations are applied to one or several populations to measure the variance of the biomass of each population following the transmission of perturbations in the entire metacommunity. Mathematically, stochastic perturbations are modeled as follows:

f_i(B₁, ..., B_S) is the deterministic part of the dynamics of the biomass of the species i, which depends on the biomass of the S species present in the metacommunity (see Equations 8a, 8b bellow). Stochastic perturbations are defined by their standard deviation σ_i and a white noise term dW_i with mean 0 and variance 1. In addition, perturbations scale with each species' biomass with an exponent z depending on the type of perturbation considered (Haegeman and Loreau, 2011; Arnoldi et al., 2019): exogenous stochasticity (from harvesting for instance) corresponds to z = 0, demographic stochasticity (from birth-death processes) to z = 0.5, and environmental factors to z = 1 (see demonstration in Lande et al., 2003 and in the Supplementary Material of Quévreux et al., 2021a). Arnoldi et al. (2019) showed that when a species is perturbed, the ratio of its biomass variance to perturbation variance increases with the species' biomass in the case of environmental perturbations, while it is independent of its biomass in the case of demographic perturbations. Therefore, Quévreux et al. (2021a) and Quévreux et al. (2021c) chose demographic perturbations in their analysis as these enabled them to perturb different species with the same relative intensity regardless of their abundance.

Synchrony in the metacommunity can be evaluated from the covariance between the temporal variations of the various species and patches, which are encoded in the variance-covariance matrix C^*. If the metacommunity is at equilibrium and perturbations are small enough, we can linearise the system in the vicinity of the equilibrium to get Equation (2) where Xi=Bi-Bi* is the deviation from equilibrium [see the Supplementary Material of Quévreux et al. (2021a) for more mathematical details].

J is the Jacobian matrix and T defines how the perturbations E_j = σ_jdW_j affect the dynamics. For instance, T_ij tells us how perturbation j affects species i. In the case of independent demographic perturbations, T is a diagonal matrix whose elements are Tii=Bi*0.5.

Then, we get the variance-covariance matrix C^* of species biomasses (variance-covariance matrix of X⃗) from the variance-covariance matrix of perturbations V_E (variance-covariance matrix of E⃗) by solving the Lyapunov (Equation 3) (Arnold, 1974; Wang et al., 2015; Arnoldi et al., 2016; Shanafelt and Loreau, 2018, and see the Supplementary Material of Quévreux et al., 2021a for a detailed description of the solution of Equation 3).

From the variance-covariance matrix C^*, whose elements are w_ij, we can compute the correlation matrix R^* of the system whose elements ρ_ij are defined by:

Correlation coefficients represent the pairwise synchrony of species i and j: if ρ_ij > 0 (correlation), species i and j are synchronous and if ρ_ij < 0 (anti-correlation), species i and j are asynchronous. Note that this metric is different from the community-wide synchrony defined by Loreau and de Mazancourt (2008), which measures the overall synchrony of many species.

Most metacommunity models consider passive dispersal for the sake of simplicity. However, many species emigrate to find food or avoid predators (Fronhofer et al., 2015; Jacob et al., 2015). Quévreux et al. (2021c) included the density-dependent dispersal function of Hauzy et al. (2010) in the metacommunity model developed by Quévreux et al. (2021a) (see Equation 10 in the Supplementary Material).

Thus, the dispersal of species i can be modulated by its own density, prey density and/or predator density. In the following, we only consider the prey density-dependent dispersal of predators (f_prey,i detailed in Equation 9a), but the two other dependencies are thoroughly described in Quévreux et al. (2021c).

S_i,j is the sensitivity of the dispersal of species i to species j biomass density, with S_0,i a constant and Bj* the biomass of species j at equilibrium. f_prey,i(B_i−1) is a decreasing function of B_i−1, which means that a higher biomass density of prey in patch #1 decreases the emigration of predators to patch #2 (negative density-dependent dispersal). In addition, varying S_0,i enables us to tune the sensitivity of dispersal to biomass density. For instance, for prey density-dependent dispersal f_prey,i, a small value of S_0,i leads to a strong response to prey biomass while a high value leads to a weak response. In the latter case, dispersal is similar to passive dispersal.

Quévreux et al. (2021c) showed that taking density-dependent dispersal into account can lead to synchrony patterns in total contradiction to those predicted by models considering passive dispersal. As an example, we consider a metacommunity with two trophic levels where only predators are able to disperse depending on prey density (see Equation 10 in the Supplementary Material), and prey in patch #1 receive stochastic perturbations. If predator dispersal is not sensitive to prey density (S0,i=103), the dynamics of prey populations are anti-correlated (Figure 4), according to the results of Quévreux et al. (2021a) presented in Figure 3. If predator dispersal is sensitive to prey density (S0,i=10-3), the dynamics of prey populations become correlated and predator populations tend to be anti-correlated (Figure 4).

To explain these correlation patterns, Quévreux et al. (2021c) proposed two mechanisms, which we detail in the case of a perturbation that increases the biomass of prey in patch #1. First, when dispersal depends on the density of several species, the species with the highest biomass CV drives dispersal, i.e., prey in the case we consider (Supplementary Figure S1-1A). Note this mechanisms is acting for prey and predator density-dependent dispersal (Figure 5A). This makes predators immigrate in patch #1 because prey are abundant there, and leads the opposite variations of predator biomass between the two patches (Supplementary Figure S1-1B). Thus, perturbations are directly transmitted from prey in patch #1 (Figure 5B ①) to patch #2 and bypass the classic vertical transmission through trophic interactions (Figure 5B ②). Finally, perturbations are transmitted vertically in patch #2 (Figure 5B ③), and the release of predator pressure in patch #2 increases prey biomass in patch #2.

Representing spatial heterogeneity is one of the major purposes of the metacommunity framework. Quévreux et al. (in prep.) represented spatial heterogeneity in the same way as Rooney et al. (2006) by increasing interaction strength by a factor γ in patch #1 (Figure 6). They also increased biomass production by primary producers by a factor ω, which is set equal to γ in the following. They showed that varying ω does not change synchrony qualitatively but strongly affects species persistence.

Quévreux et al. (in prep.) mainly showed that, in heterogeneous metacommunities, synchrony depends on which patch is perturbed, i.e., perturbing the population of a species in patch #1 or #2 does not lead to the same correlation pattern between patches. Here, we consider a metacommunity where interaction strength is higher in patch #1 and only predators are able to disperse. Increasing asymmetry of interaction strength γ increases the correlation of prey populations when prey are perturbed in patch #1 while it decreases it when prey are perturbed in patch #2 (Figure 7). This discrepancy is due to the source-sink dynamics of predators (predator biomass produced in excess in patch #1 spills over patch #2), which increases apparent competition for prey and leads to a lower prey biomass in patch #2 than in patch #1 (Figure 8 ① and Supplementary Figures S1-2A,B). This alters predator-prey dynamics: when prey are perturbed in patch #1, high prey biomass and high interaction strength generate strong responses of prey and predator biomasses (Figure 8A ② and Supplementary Figures S1-2C) that are transmitted to prey in patch #2, which are completely under predator control (Figure 8A ③). This increases synchrony because predators drive the dynamics in the two patches. When prey are perturbed in patch #2, predator-prey dynamics are weak because of low prey biomass and weak interaction strength (Figure 8B ②). Thus, perturbations are weakly transmitted to patch #2 (Figure 8B ③) and prey populations slowly return to equilibrium, which leads to asynchrony (Supplementary Figure S1-2C).

Dispersal is a complex process that has received much attention in the last two decades in order to understand the dynamics of metacommunities. Many mechanisms have been identified (Amarasekare, 2008; Gross et al., 2020), but most of them have been studied with different models (metacommunity or coupled food webs) displaying different dynamical regimes (equilibrium or limit cycles) and with different measures of stability (asymptotic resilience, occurrence of Hopf bifurcation, or temporal variability). Here, we propose a general framework based on a food web model that covers a broad range of ecological contexts, which enables us to derive robust conclusions. In addition, the flexibility of the model, which can be easily modified, and the use stochastic perturbations, which can target specific trophic levels in a particular patch, enable us to have a deep understanding of the mechanisms acting in the metacommity. Most previous studies described the factors governing synchrony and stability without explaining the mechanisms underlying the observed effects. For instance, Hauzy et al. (2010) and Rooney et al. (2006) did not explain how density-dependent dispersal and asymmetry affect synchrony.

Our model fills this gap: it is easy to handle because it is linearised and most of the mechanisms have additive effects. For instance, the effects of density-dependent dispersal and spatial heterogeneity can be disentangled properly. It appears that spatial heterogeneity has a weaker effect on predator (species 3) synchrony than density-dependent dispersal since the overall correlation of their populations is qualitatively similar to that obtained in a model with density-dependent dispersal only. Thus, in Rooney et al. (2006) model, prey preference is the major driver of synchrony and structural asymmetry, which is how they defined spatial heterogeneity, has minor effects. The comparison with Rooney et al.'s (2006) model stresses the continuum between coupled food chains and metacommunities. Many studies consider a unique population of mobile top predators feeding on two distinct energetic channels (Post et al., 2000; McCann et al., 2005; Rooney et al., 2006; Vasseur and Fox, 2007; Attayde and Ripa, 2008), but these models are actually a limiting case of a metacommunity in which top predators have a density-dependent dispersal and a high dispersal rate.

Our results based on a two-patch metacommunity may well explains what would happen in a larger landscape with many patches since the overall response to perturbations would be the additive effect of pairwise connections between patches. Synchrony in multi-patch metacommunities has been addressed by models with limit cycles (Marleau et al., 2014; Hayes and Anderson, 2018), but the response of such large metacommunities to perturbations remains to be explored.

Synchrony in metacommunities is tightly linked to trophic interactions and biomass distribution. Barbier and Loreau (2019) and Quévreux et al. (2021a) showed that perturbations tend to have a bottom-up (top-down) transmission when biomass pyramid is bottom (top) heavy. This is mostly due to the distribution of metabolic rates among species: when m is high, high trophic levels have a faster metabolic rate than lower trophic levels, which accelerates their biomass dynamics (Equation 8b). Thus, these fast top predators will quickly dampen perturbations while perturbed lower trophic levels with a slower pace of life will have a stronger impact on the rest of the food chain. In addition, some types of perturbations can affect species depending on their abundance. Arnoldi et al. (2019) showed that immigration perturbations strongly affect rare species while environmental perturbations strongly affect abundant species. Demographic perturbations are generally used because they affect species with the same intensity regardless of their biomass but environmental perturbations are more relevant from an ecological point of view. Thus, synchrony in metacommunities under environmental stressors affecting every species may be more easily understood by tracking the transmission of perturbations affecting the most abundant species.

The effect of dispersal is conditioned indirectly by trophic dynamics. Quévreux et al. (2021a) showed that perturbation transmission depends on the relative importance of dispersal compared with the other local demographic processes. Thus, in a metacommunity, dispersal will have the biggest impact not on the species with the highest dispersal rate but on the species for whom dispersal has the highest relative contribution to biomass dynamics. This distinction is critical because many studies reported that the dispersal of particular trophic levels is stabilizing (Koelle and Vandermeer, 2005; Pedersen et al., 2016).

The trophic context of dispersal is even more important for density-dependent dispersal. Quévreux et al. (2021c) showed that vertical transmission through trophic interactions is bypassed only if the species interacting with the dispersing species have a higher biomass CV than the dispersing species (Figure 5A). Species biomass CV depends on perturbation intensity but also on its transmission across the food chain. Shanafelt and Loreau (2018) demonstrated that in food chains displaying a trophic cascade, species experiencing a strong top-town control have a high biomass CV because they have a low biomass. Thus, trophic cascades also translate into stability cascades with species being more or less stable depending on their trophic position. In the case of density-dependent dispersal, these rare species should have a major impact on dispersal in many situations.

Finally, we were able to draw general conclusions from the wide range of responses displayed by the model. A key factor of Barbier and Loreau (2019) model is self-regulation, which encapsulates various negative intra-specific interactions such as territoriality, cannibalism and disease transmission that lead to a density-dependent mortality rate. Self-regulation remains neglected in many food web models, although it has been identified in natural marine phytoplankton communities (Picoche and Barraquand, 2020), where it promotes species coexistence (Picoche and Barraquand, 2019) and food web stability (Barabás et al., 2017) by dampening interspecific interactions relative to intraspecific interactions. Thus, self-regulation definitively interferes with dispersal and must be taken into account in future metacommunity models.

As we explained above, our model is easy to handle because it is linearised and thus most of its responses to perturbations have additive effects. However, this may not hold with systems that are not at equilibrium or subject to perturbations pushing them to far from equilibrium. For instance, Post et al. (2000), McCann et al. (2005), Hauzy et al. (2010), and Marleau et al. (2014) did not consider self-regulation and their food web models were characterized by strong inter-specific interactions leading to limit cycles.

In systems with limit cycles, phase-locking, which is a constant phase difference between patches through time, can occur and lead to emergent metacommunity properties (Jansen, 1999; Lloyd and May, 1999). For instance, Goldwyn and Hastings (2009) showed that even a 0.1% variation in species parameters among patches (small spatial heterogeneity) is enough to generate substantial asynchrony. Two components of the structure of the metacommunity are critical for its stability. First, the dispersal of particular trophic levels determines the occurrence of limit cycles and the synchrony of their phase (Koelle and Vandermeer, 2005; Hauzy et al., 2010; Pedersen et al., 2016). Second, in spatially explicit metacommunities, connectivity between patches (Guichard et al., 2019) and their spatial organization (Marleau et al., 2014; Hayes and Anderson, 2018) drive phase synchrony. For instance, in riverine metacommunities, patch linked to multiple upstream and downstream patches have communities with less temporal variability (Anderson and Hayes, 2018).

These responses of oscillating systems are not disconnected from responses to stochastic perturbations. Vasseur and Fox (2007) demonstrated that small synchronous stochastic perturbations can continuously deviate the system from its regular oscillations and dampen predator dynamics. Thus, stochastic perturbations and limit cycles can strongly interact to drive food web stability. Finally, in systems at equilibrium but with strong predator-prey interactions (high values of εa and ma), perturbations lead to dampened oscillations (see Rooney et al., 2006 and Supplementary Figure S1-2). Large perturbations affecting this kind of system could then lead to responses at the interface between stochastic perturbations and limit cycles. Thus, the theoretical framework used to understand synchrony in systems with limit cycles may be used in future studies to understand the response to perturbations of systems close to the Hopf bifurcation.

Theoretical models help us to understand the response of ecosystems to perturbations, but empirical tests are required to confirm the effects of the described mechanisms in real systems. Fortunately, empiricists have already developed methods and equipments for manipulating metacommunities. Microcosms offer immense possibilities because of their high level of control, as extensively described by Altermatt et al. (2015). Bacteria and protists represent huge populations with predator-prey dynamics that can be quantified efficiently by flow cytometry and image analysis to record fine scale times series (Vasseur and Fox, 2009; Fox et al., 2013). For instance, Vasseur and Fox (2009) manually drove dispersal by pipetting medium from one microcosm to another and applied stochastic perturbations by randomly changing temperature. Harvey et al. (2016) and Harvey et al. (2020) coupled microcosms by exchanging subsidies (microwaved aliquots of medium) and oriented dispersal by asymmetric transfer of medium between microcosms.

Larger facilities are required to perform experiments with populations of large mobile organisms such as flying insects and vertebrates. The terrestrial metatron consists of 48 interconnected caged patches (10 × 10 × 2 m) where temperature, humidity and illuminance can be controlled (Legrand et al., 2012). There, Fronhofer et al. (2018) measured density-dependent dispersal of butterflies and damselflies in the presence of predators such as toads and frogs. Its aquatic equivalent consists of 144 basins of 2 m³ each that can be connected by aquatic or aerial corridors and in which temperature can be controlled. The PLANAQUA experimental facility includes an aquatic macrocosm platform consisting of 16 artificial lakes with an individual volume of 700 m³, connected to each other by dispersal corridors, and equipped with an automated instrumentation (Mougin et al., 2015). With such a platform, long-term experiments with large fish populations can be performed in a controlled environment, which is close to a natural shallow pond integrating littoral, benthic, and pelagic zones. However, a strong sampling effort is required to get the time series of predator-prey dynamics and assess their temporal stability (Rakowski et al., 2019).

We summarized a few mechanisms that govern synchrony in metacommunities but these mechanisms are population centered and do not account for the rest of the ecosystem. Ecosystems are linked by flows of matter and energy (Gruber and Galloway, 2008; Gounand et al., 2018b), turning metacommunities into metaecosystems (Loreau et al., 2003; Massol et al., 2011; Guichard, 2017; Gounand et al., 2018a). Nutrient cycling strongly affects ecosystem dynamics (DeAngelis, 1992; Loreau, 1994, 2010; Quévreux et al., 2021b) and the response of local ecosystems to perturbations (Theis et al., 2022). Nutrient flows between ecosystems create subsidies that alter local dynamics and can be destabilizing by increasing the amplitude of biomass oscillations (Marleau et al., 2010), while the dispersal of particular trophic levels in such metaecosystems can dampen these negative effect (Leroux and Loreau, 2012; Gounand et al., 2014). However, the general response of metaecosystems to stochastic perturbations is yet to be studied.

In the context of global changes, in a world where habitat is fragmented and subject to anthropogenic perturbations, understanding the drivers of species dispersal and how spatial structure shapes the dynamics at landscape scale are crucial for conservation policies. The metacommunity framework is a powerful tool to achieve this goal because it enables us to disentangle the various mechanisms governing the synchrony of species dynamics at local and regional scales. This paper has summarized some aspects of metacomminities such as density-dependent dispersal and spatial heterogeneity, which are key to predict the response of metacommunities to perturbations. These mechanisms can be studied in isolation and their individual effects can be added up to have an insight in the response of realistic metacommunity to perturbations. Ultimately, our results help us to identify the species and the patches for which perturbations must be mitigated to improve stability at local and regional scales.