Illuminating the intrinsic and extrinsic drivers of ecological stability across scales

Department of Zoology, School of Natural Sciences, Trinity College Dublin, Dublin, Ireland Biodiversity and Biocomplexity Unit, Okinawa Institute of Science and Technology Graduate University, Okinawa, Japan Graduate School of Life Sciences, Tohoku University, Sendai, Japan Integrated Bioresource Information Division, RIKEN BioResource Research Center, Ibaraki, Japan Biological Complexity Unit, Okinawa Institute of Science and Technology Graduate University, Okinawa, Japan Centre for Biological Diversity, University of St Andrews, St Andrews, UK

K E Y W O R D S community dynamics, empirical dynamic modeling, energy landscape analysis, metacommunity, time series 1 | ECOLOGICAL STABILITY ACROSS SCALES Ecological stability (Table 1) considers how ecosystems resist or recover from environmental change Donohue et al., 2013;Hillebrand et al., 2018;Holling, 1973). Understanding the drivers of ecological stability is then key to predicting how communities and ecosystems respond to disturbance (Holling, 1973;Pimm, 1984). Yet, the concept of stability is complex and multidimensional, requiring clear definition to be of value as a quantitative tool Kéfi et al., 2019). There exists a long history of measuring ecological stability in the temporal domain (MacArthur, 1955;Pimm, 1984), either using time series data to quantify the temporal variability of some measure of interest (e.g., Tilman et al., 2006), or by measuring an ecological response to a specific perturbation (e.g., Solé & Montoya, 2001; Figure 1a). However, perturbations have both a spatial and temporal signature which influences ecological responses (Pinek et al., 2020;Ryo et al., 2019;Waldock T A B L E 1 Glossary of key terms Community network structure The composition of pairwise interspecific interactions in a community; i.e., the graph representing structure of interactions between species in a community. Community state The composition of species in a community at a given point in time. Alternative stable states represent stable coexistence of multiple community states under the same environmental conditions (Fukami, 2015). Attractive basins represent domains in state space in which all initial conditions will tend towards an equilibrium state (Walker et al., 2004).
Compositional variability Variability in community composition through time (i.e., community state) and/or space.

Demographic stochasticity
The probabilistic dynamics (in time and/or space) of vital rates (e.g., births, deaths) of individuals underpinning population growth rates (Neubert & Caswell, 2000).
Ecological stability A multidimensional concept typified either as a measure of variability in time and/or space (e.g., Tilman et al., 2006) or as the response to a perturbation, quantified as resistance, resilience, recovery etc. . To operationalize the stability concept, it is important to clarify the stability dimension(s) of interest (Donohue et al., 2013).
Empirical dynamic modeling A data-driven framework to model nonlinear dynamic systems from time series data. It is equation-free and captures system dynamics and relationships from the attractor manifold reconstructed by taking time-lag coordinates of the observed variables (see Deyle et al., 2016).
Energy landscape analysis A method to determine the likelihood of a given (community) state. The method assigns an energy value to each state, where states with low energy are taken with a large probability during stepwise state transitions, e.g., community assembly (Watanabe et al., 2014).
Environmental variability Variability in the environment through time or space. Includes the dynamics of this variability (environmental stochasticity/noise) and the range of environmental conditions experienced over a given spatial and/or temporal scale (environmental heterogeneity).
Extrinsic driver Drivers of ecological dynamics external to the biological component of the system (e.g., environmental stochasticity).
Interaction variability Variability in pairwise or multi-species interspecific interactions through time and/or space. This variability could include variation in the identity, strength, density-dependence, asymmetry etc. of interactions determining community network structure.
Intrinsic driver Drivers of ecological dynamics within the biological component of the system (e.g., demographic stochasticity, community network structure, interaction density-dependence).

Metacommunity
A set of local communities linked by dispersal of multiple potentially interacting species (Leibold et al., 2004).
Spatial/temporal (sampling) extent The total area of observation encompassing all observation units in space or time (Viana & Chase, 2019). Zelnik et al., 2018), and communities may respond heterogeneously in space. The increasing frequency and intensity of such disturbances in the Anthropocene justify continued efforts by ecologists to quantify temporal stability across a range of scales (Jacquet & Altermatt, 2020, Newman, 2019, Wolkovich et al., 2014Figure 1b). Here, we discuss the consequences of scale-the spatial or temporal window size of a sampling unit or sampling extent (Table 1)-in studies of ecological stability and dynamics, and highlight nascent methods with the potential to better integrate scale into stability research. Existing models in community ecology have integrated spatial context, including dispersal and environmental heterogeneity, instigating the development of metacommunity theory (Hubbell, 2001, Leibold et al., 2004, Vellend, 2016 Table 1). This framework is still being refined, with recent models aiming to encompass multiple concepts and processes across different spatial and temporal scales . Such models spurred an extensive discourse contrasting selection and niche versus neutrality in ecological communities (e.g., Legendre et al., 2009;Mori, 2018). Yet, the mechanisms by which such metacommunity processes influence community structure in systems with high spatial complexity or intricate temporal dynamics are still poorly understood. Recent technical advances such as remote sensing (e.g., Bush et al., 2017) and ecological data synthesis (e.g., Dornelas et al., 2018), in combination with metacommunity theory Vellend, 2016), should allow us to untie interlinked processes that shape metacommunity dynamics. Interactions between processes, including between selection and dispersal (Mouquet & Loreau, 2002, and dispersal and drift (Ron et al., 2018), are increasingly better studied, as is the link between biodiversity and stability in the metacommunity context Wang & Loreau, 2016). Furthermore, the interaction between various ecological processes and dispersal hinges on the specific environmental conditions and spatial scale considered, thereby shaping biodiversity patterns and ecological responses idiosyncratically (Mori, 2018;Ratcliffe et al., 2016). This means environmental context-one source of extrinsic variation-is F I G U R E 1 A scale-dependent understanding of ecological stability. (a) Ecological stability can be measured with multiple metrics using a range of ecological variables (including measures of biodiversity, ecosystem functioning, or ecosystem service supply). Variability (1) is often measured as the coefficient of variation of a variable in time or space (Tilman et al., 2006); resistance (2) is the degree to which a variable of interest is altered in response to a perturbation (dashed line) (Pimm, 1984); recovery time (3) and recovery rate (engineering resilience) (4) are the time (or space) taken for a value to reach some threshold and the rate of return to this threshold, respectively (e.g., Barnthouse, 2004). There are many more dimensions of stability including persistence and asymptotic stability each with specific definitions and formulations Donohue et al., 2013). (b) When considering the environment in which ecological communities are embedded, the spatial (x-axis) and temporal (y-axis) scale of measurements influences the form of the environment-time (or space) relationship, which in turn results in scale-dependent extrinsic (environmental) variation (inset plots).
(c) Environment is thus a function of both time and space (x-axis) and the complexity of this environment (including our scale-dependent consideration of extrinsic variation) can hinder our ability to understand how such environments influence ecological stability (different lines) [Color figure can be viewed at wileyonlinelibrary.com] unavoidably coupled with the processes shaping community stability through their interaction with spatial scale (Chase, 2003, Gonzalez et al., 2020Figure 1b,c). Studies of community dynamics and the processes that shape them must therefore seek to reveal dynamics at the spatial scale of interest. Identification of the dominant processes driving ecological dynamics is key for modeling and prediction. The dominant processes shaping population or community dynamics depend explicitly on the spatial scale considered and on environmental context, including the temporal signal of environmental conditions (Pinek et al., 2020;Villa Martín et al., 2019;Zelnik et al., 2018). Patterns of biodiversity change in ecological systems are also scale dependent . Thus, spatial and temporal dynamics and their drivers interact to produce a complex stability landscape in space and time. However, work todate has mostly focused on a single or few spatiotemporal scales Waldock et al., 2018), demanding increased focus on and knowledge about how scale interacts with the processes shaping ecological communities (Gonzalez et al., 2020, McGill et al., 2015Figure 1b).
Here, we discuss metrics and emerging tools in studies of ecological stability from theoretical and empirical perspectives, with a particular focus on temporal and spatial context. Our aim is not to provide a novel framework for considering ecological dynamics in light of spatiotemporal scale, but rather we highlight recent methodological advances with the potential to better integrate scale into studies of ecological stability and dynamics. Particularly, we describe conditions under which spatial or temporal scale may influence our understanding of various intrinsic (interaction variability, compositional variability) and extrinsic (environmental variability) drivers (Table 1) of ecological dynamics, and some methodological considerations using recent case-studies at different levels of ecological organization. We describe methods by focusing on the main source(s) of variability being imposed on ecological communities, which often depend on spatiotemporal scale. We consider first how variation in interspecific interactions (an intrinsic driver of community dynamics) are a function of time and provide some suggestions for applying empirical dynamic modeling to ecological time series. Next we move from pairwise interspecific interactions to consider the entire suite of species in a metacommunity, focusing on how spatial scale influences compositional variability as an intrinsic driver of assembly dynamics. We then move to the determinants of species' ranges at a global scale by discussing how population expansions, adaption and species invasions may be driven by the fluctuating nature of (extrinsic) environmental conditions. Finally, we touch on how these sources of intrinsic and extrinsic variation interact to collectively shape population, community and ecosystem dynamics, and how scaledependence might propagate across levels of ecological organization to influence our understanding of how global change influences ecological dynamics and stability.

| COMMUNITY DYNAMICS: INTERSPECIFIC INTERACTION VARIABILITY
Studies of ecological stability typically focus on aggregated variables at the community level . At this level of organization, interspecific interactions are a principle driver of community dynamics and stability. Pairwise interactions (e.g., predator-prey) have long been a focus of population ecologists (e.g., Berryman, 1992). However, in nature, species do not simply interact in pairs. Rather, species' population dynamics are influenced by densities of various other species, including their prey, predators, competitors and mutualistic partners, as their vital rates are affected by those interacting species. An ecological community can therefore be viewed as a community network of interspecific interactions, in which a change in one population leads to cascading changes in other populations. Community network structure (Table 1) can be quantified as the number of other species a given species interacts with (i.e., degree distribution, Newman, 2003), or based on motifs or modules, each of which is a subset of a network with a small number of interactions and species arranged in a specific way, to reveal indirect interspecific interactions (Amarasekare, 2008;Simmons et al., 2019). Community network structure defines the nature of interactions within the community by determining the paths through which interspecific effects are transmitted. Network structure, therefore, has the potential to influence dynamics and stability at the community level (Higashi & Nakajima, 1995;Yeakel et al., 2020).
Community network structure has been a major research theme in recent decades, particularly when considering the relationship between network complexity and ecological stability. May (1972) predicted that a community with more species or more interactions, that is, a more complex community, has less chance to be stable (in terms of reproducible periodic oscillations of the populations therein). How complex community networks arise and persist has since been a major theme of theoretical ecology, and many factors are candidates for stabilizing complex community networks, including weak interactions (McCann et al., 1998), flexibility in interactions arising from adaptation (Kondoh, 2003), nestedness (Thébault & Fontaine, 2010), mixtures of antagonistic and mutualistic interactions (Mougi & Kondoh, 2012), meta-community structure (Mougi & Kondoh, 2016) and redundancy among species effects on environmental conditions (ecosystem engineering, Yeakel et al., 2020). Moreover, these network properties and their relationship with community dynamics are likely influenced by spatial and temporal scale. For example, the presence, directionality or strength of interactions between species vary through timeparticularly across species' life-stages (Yang & Rudolf, 2010) or during community assembly (Yeakel et al., 2020)-which leads to time dependence in network structure (Deyle et al., 2016;Moreno-Mateos et al., 2020).
It remains a challenge to empirically test the hypothesis that community network structure itself can be a determinant of community dynamics and stability. There are several fundamental difficulties in performing such empirical tests. The first is in identifying the interspecific interactions necessary to build a community network. Presence of an interspecific interaction between Species X and Species Y is defined by the fact that the density of Species X affects the population dynamics of Species Y. However, these cases are hard to examine in nature. Even with the continued rise in field experiments manipulating community composition, the fully resolved community network in which species are embedded is rarely available (Pringle & Hutchinson, 2020). That is, it is still difficult to identify all pairwise combinations of interspecific interactions in a natural community. Empirically tracking the dynamics of many-species systems involves considerable effort. Even if monitoring the population dynamics of many species, it is not straightforward to isolate community stability-which intrinsically arises from the full suite of interspecific interactions between all species in a community-from variability caused by environmental noise (Table 1). A further difficulty is in identifying causation in the link between structure and dynamics. Even with success in identifying interaction network structure and tracking community dynamics, how best can the two be related?
The minimalist approach (Ushio et al., 2018) is a framework recently proposed to overcome such difficulties, enabling the empirical study of community network structure and dynamics via empirical dynamic modeling (Table 1). Empirical dynamic modeling can be used to predict temporal change in community network structure using time series of multispecies abundance data. The method assumes that system dynamics are deterministic while randomness arises from chaotic behavior of nonlinear dynamics. What enables the minimalist approach is a model that allows interactions to be inferred from available time series data (Chang et al., 2017). Time series have rich information about the processes that drive community dynamics. For example, several authors have applied Takens' embedding theorem and related state space reconstruction theories (Deyle & Sugihara, 2011;Takens, 1981) in various ways, such as in determining causal relationships between variables (Sugihara et al., 2012), forecasting (Sugihara & May, 1990), quantifying temporally varying interaction strengths (Deyle et al., 2016) and for evaluating the stability of community trajectories (Ushio et al., 2018). By combining time series data and appropriate modeling methods, one can infer interspecific interactions, their signs and strengths, and even evaluate the dynamical consequences of community structure using only observational data from the field. Empirical dynamic modeling also has the potential to reveal higher-order interactions, where an interaction between two species is moderated by the presence of a third (Bairey et al., 2016), because information encoded in this unobserved driver is embedded in the observed time series (Takens, 1981, but see Blanchet et al., 2020). Ushio et al. (2018) recently used the minimalist approach to determine the structure of large community networks based on 12 years of fortnightly coastal fish data. In this study, the analysis of a multi-species observational time series allowed detection of 14 interspecific interactions, of which some were positive and others negative. Quantification of interaction strengths revealed the temporally varying nature of interspecific interactions; some interactions fluctuated between positive and negative values over time. The authors not only found that temporally varying interspecific interactions produced temporally fluctuating local Lyapunov stability (therein termed dynamic stability, Ushio et al., 2018), but also provided empirical evidence that community network complexity is important for community stability. More specifically, higher species diversity led to a dominance of weak interactions, which in turn stabilized community dynamics (Ushio et al., 2018). Critically, Ushio et al. (2018) provide evidence of the potential for empirical dynamic modeling to successfully link complexity and stability in a causal framework using only field data.
Interspecific interactions, community network structure and stability all depend on the spatial or temporal scale of observation. Population dynamics which are each unstable over short temporal scales (e.g., monthly) can give rise to stable dynamics at another temporal scale (yearly). Similarly, unstable dynamics at local spatial scales can be stable when considering a regional scale, as suggested by classical meta-population theory and extensions of the theory across levels of organization (e.g., Hanski, 1998;Wang & Loreau, 2014. Given the scale-dependent nature of ecological dynamics, the minimalist approach has several desirable characteristics for the study of community structure and stability; interactions and their long-term changes that actually drive dynamics at the spatiotemporal scale of choice can be identified because time series of interspecific interactions are inferred from observed short-term (transient) interactions (multi-species abundance time series). That being said, when using empirical dynamic modeling to understand community network complexity, it is imperative that species' time series be sampled at a temporal scale that matches that of the underlying mechanisms of interest (e.g., dispersal vs. birth rate); it is not enough simply to select an appropriate temporal scale for a target species' life history or a system's disturbance regime, as the mechanisms behind community dynamics and both the properties of and responses to disturbance vary with temporal scale (Wolkovich et al., 2014). Short time series have little power to reveal meaningful dynamics, because this method requires time series to be several times longer than the system return time, that is, the time taken to return to a local attractor (Munch et al., 2019). This led to the suggestion that time series longer than 30 data points are most appropriate for use in empirical dynamic modeling, as this length reduces prediction error significantly (Chang et al., 2017;Munch et al., 2019;Sugihara et al., 2012), although this suggestion remains naïve to specific disturbance regimes. Conversely, long time series may violate the model's assumption of stationarity (Munch et al., 2019). System dynamics themselves are also inherently scale-dependent, meaning that analyses at different temporal scales can produce different results (Chang et al., 2017), and that larger sampling intervals have the potential to miss relevant dynamics leading to apparent indeterminism (Munch et al., 2019). Finer resolution sampling intervals have the advantage of accurate detection of transient ecological interactions, while also increasing cross-map skill (coherence) as time series length increases (Sugihara et al., 2012). A key open question remains how empirical dynamic modeling can be applied and integrated across scales to capture the dynamics of species with vastly different life histories.
It is clear, then, that temporal scale plays a critical role in the variation in pairwise interspecific interactions, and the consequences of this variation for community network dynamics and stability. Thinking next about the composition of ecological communities rather than interspecific interactions therein, spatial scale is a key driver of compositional variability and thus the intrinsic dynamics and stability of community assembly.

| COMMUNITY ASSEMBLY DYNAMICS: COMPOSITIONAL VARIABILITY
The next level up is composed of networks of interacting communities or metacommunities. At this level of organization, a key driver of stability is the dynamic interaction between the composition of species in a community (hereafter a community state, Table 1) and habitat types in different patches. Metacommunity models often focus on patch dynamics which arise from shifts in the balance of local extinction and (re)colonization and result in a stable regional species pool (Leibold et al., 2004). Because of patch dynamics, community composition may pass through transient states as community assembly progresses . This results in sites harboring communities of different developmental stages at a given point in time. A stable state may no longer be perceived as a permanent endpoint, but instead may represent a dynamic equilibrium recognized in cycling community states.
Typical approaches to understanding species cooccurrences, such as joint species distribution models (Tikhonov et al., 2017), do not directly describe dynamic relationships between different community states (Baselga & Araújo, 2009;Elith & Leathwick, 2009;Norberg et al., 2019). These approaches cannot then be applied as models for community assembly dynamics. Dynamical models, such as differential equations, are able to describe shifts in community state based on the changes in species abundance (Gravel et al., 2011). However, it is often difficult to develop fully mechanistic models for multi-species communities from abundance time series. To overcome these difficulties, it is possible to integrate a pairwise maximum entropy model (i.e., a Markov network, Araújo et al., 2011, Harris, 2016 with energy landscape analysis (Wales et al., 1998, Watanabe et al., 2014; Table 1, Figure 2).
Briefly, consider the probabilities of observing different community states, which describe how often a specific combination of species is observed and are adjacent to other states, that is, differing in the presence-absence status of a single species. Community assembly is then assumed to proceed in a stepwise manner through this network of possible states. With this network of possible community states and their frequencies-plus abiotic variables if available-an energy landscape can be calculated where energy represents the probability of observing a given community state in nature; states with lower energy have a greater probability of being realized (Watanabe et al., 2014). The form of the energy landscape can then provide systematic understanding of the developmental pathways that constrain natural community assembly. Applied to ecology, the method starts with community data including observations on the occurrence of species in a set of temporal and/or spatial samples (Figure 2a), converted to a presence-absence community matrix (Figure 2b). Using this matrix, a pairwise maximum entropy model is then produced, which estimates the probability with which a given community state is realized (Figure 2c), and energy landscape analysis can be applied on the product of this model (Figure 2d). This process identifies possible stable states and tipping points and can be used to form a disconnectivity graph summarizing hierarchical relationships between the stable states and possible tipping points, defining an attractive basin of stable states (Table 1). It also allows emulation of community dynamics constrained on the energy landscape. This method considers only the consequence of species extinctions or migrations in local communities, that is, changes in presence-absence of species. Hence it ignores population dynamics and approximates the continuous state space as a network of all possible community states. By discarding detailed descriptions of population dynamics, this approach offers a practical way to study community assembly dynamics from observational data. As this method uses species occurrence data to infer an energy landscape instead of directly estimating ecological interactions-though the imprint of ecological interactions is included in the pairwise maximum entropy modelresearchers can avoid the issue that species co-occurrence may not inform ecological interactions (Blanchet et al., 2020), since the pairwise maximum entropy model is treated as an approximate model of governing dynamics (e.g., population dynamics) rather than a predictive function as in the case of machine learning. F I G U R E 2 Energy landscape analysis of community assembly pathways. (a) Appropriate datasets for energy landscape analysis include occurrence of species in local communities sampled from multiple sites and/or timepoints, possibly accompanied by values representing local abiotic environment (explicit abiotic factors). Colors and sizes of ellipses represent differences in local environmental conditions. (b) The dataset is converted to matrices of presence/absence (black/white) status, and explicit abiotic factors (if available) (colors show values of abiotic factors). (c) These matrices are used to fit parameters in a pairwise maximum entropy model. (d) The fitted model specifies an energy landscape network with nodes representing community states and links representing transitions between community states. The energy landscape analysis acknowledges: (i) the stable states (purple/yellow nodes) and tipping points (green node), (ii) a disconnectivity graph as the summary of hierarchical relationships between the stable states and tipping points, (iii) attractive basin(s) of stable states (purple/yellow nodes indicate attractive basins of the two stable states). (iv) We are also able to emulate assembly dynamics constrained on the energy landscape. Given how spatial scale interacts with the drivers of community assembly (e.g., Chase, 2014), the results of energy landscape analysis should be interpreted in light of the spatial scale of both the sampling unit and total sampling extent (Viana & Chase, 2019). If sampling units are at an appropriate spatial or temporal scale for a driving process of interest to operate, single or multiple attractive basins may be observed (Figure 3). However, if the sampling units are too small or too large, these basins may be less distinct, even if the interaction between species is strong. Moreover, explicit consideration of spatial extent in studies of energy landscape analysis can disentangle the drivers of community assembly across scales. If the spatial extent of observations is large and the effect of environmental heterogeneity significant, by choosing a relevant environmental parameter (or appropriately reducing multiple parameters into a single dimension), the shift of stable community states can be unfolded with respect to the parameter. One could then explore whether the model's explanation is consistent with our processbased understanding of community assembly dynamics (Figure 3). The results of energy landscape analysis across spatial scales should reveal that under a small spatial extent, compositional variability (Table 1) stems from both local biotic interactions and stochastic processes, but at a larger spatial extent, compositional variability depends more on environmental conditions including habitat connectivity. Consequently, energy landscape analysis potentially offers a practical way to understand the drivers of community dynamics across spatial scales. In principle, this approach could also be applied across temporal scales to reveal temporal contingencies in community assembly drivers . Such an approach would complement existing methods proposed to assess the relative importance of the ecological processes that drive community assembly (D'Amen et al., 2018;Legendre et al., 2009;Mertes & Jetz, 2018;Meynard et al., 2013).
While the outcome of methods described in this section can reflect environmental variation, the focus of such methods is on understanding how community state changes dynamically. In the next section, we shift our focus to methods for revealing the processes through which environmental variability contributes to community assembly through single species population dynamics.

| POPULATION DYNAMICS: ENVIRONMENTAL VARIABILITY
Community assembly is driven by a collective of multiple species' population establishments and dynamics, but the interaction between spatiotemporal scale and dynamics is just as relevant when considering each species in isolation. Just as community assembly can be driven by intrinsic characteristics of the community, populations are governed by their intrinsic demographic rates (e.g., births, deaths, dispersal). We have so far overlooked a key driver of dynamics and stability across multiple levels of biological organization: environmental variability (Table 1). At the global scale, this extrinsic variation has the ability to shape species' distributions since virtually all mobile species can move towards suitable exogenous conditions, thereby expanding their range. Range expansions into new environments are often a consequence of environmental change. Global environmental change is rearranging species distributions through population expansions, contractions and poleward shifts (Chen et al., 2011;Parmesan et al., 1999). Such population changes can occur either as a direct result of individual-level actions (e.g., migrating to more suitable environmental conditions) or as a consequence of human-mediated dispersal which can introduce species to otherwise unreachable habitats.
Predicting the stability of ecosystems undergoing range expansions and continuous population rearrangements is a challenge, not only because of the difficulty in inferring causation, but also because the study of range expansions F I G U R E 3 Relationship between the spatial extent of observation and the stable state diagram resulting from energy landscape analysis. (Top) When a larger spatial extent is considered, a larger potential range of abiotic conditions is captured (colors). In turn, a larger range of community states can be observed on the energy landscape. (Bottom) This corresponds to a stable state diagram representing a more comprehensive range of attractive basins under different abiotic conditions. The black lines in the figure represent the position of the stable state in the energy landscape for the abiotic conditions represented by the horizontal axis [Color figure can be viewed at wileyonlinelibrary.com] must generally consider large spatiotemporal scales. Until recently, large-scale spatial and temporal data on population diversity and abundance were scarce and often inaccurate, but the continued improvement of global collaborative networks (e.g., BioTIME, Dornelas et al., 2018) facilitates a deeper understanding of the causes and consequences of population rearrangements.
When considering the factors underpinning population dynamics, scale-dependencies are pervasive. For example, spatiotemporal scale affects the demographic processes influencing dynamics at the population level. Likewise, extrinsic factors such as environmental suitability vary as a function of spatiotemporal scale to further moderate population dynamics. The larger the sampling scale considered, the greater the likelihood of detecting environmental variation in space or time (Wiens, 1989). Consequently, environmental change is characterized not only by its mean and variance, but by the dynamics of its variance in space and time (its stochasticity, Shoemaker et al., 2019; Table 1), its rate of change (Pinek et al., 2020) and its velocity, which comprises both the spatial and temporal dimensions of environmental variation (Burrows et al., 2014). These heterogeneities affect both deterministic and stochastic processes influencing ecological populations (Shoemaker et al., 2019;Waldock et al., 2018) and thus, greater effort should be made to incorporate such heterogeneities into ecological theory and experiments.
Another factor contributing to range expansions is the diversity of the expanding population. Although expansion mechanisms of phenotypically diverse populations are well explored (Barto n et al., 2012;Fu et al., 2018;Keller & Segel, 1971;Neubert & Caswell, 2000), there has been little emphasis on the potential impact of spatial or temporal environmental heterogeneities on such mechanisms. Recently, Villa Martín et al. (2019) considered how environmental heterogeneity can modify the success of population range expansions. In particular, they considered a population front for which expansion velocity depends on environmental conditions that vary in time or space at different frequencies. Under infrequent environmental variation, phenotypically diverse populations expanded faster than homogeneous populations for a wide range of environmental conditions. However, under frequent environmental variation, homogeneous populations expanded faster than did populations with higher phenotypic diversity. Moreover, they found that diversity had a stronger effect on expansion velocity under spatial rather than temporal environmental heterogeneity (Villa Martín et al., 2019).
Environmental heterogeneity, and its structure in space and time, emerges as a critical determinant of population expansions. Temporal and spatial environmental variability in combination cause nontrivial effects on the success of population expansions. For example, slow expansions through unfavorable environments will have either stronger or weaker effects on expansion velocity depending on the temporal fluctuations of the environment therein (Villa Martín et al., 2019). To better relate environmental conditions to population dynamics and subsequently to ecological stability, it is necessary to recast the multiple dimensions of environmental change as a single disturbance framework (Jacquet & Altermatt, 2020;Newman, 2019;Pinek et al., 2020;Shoemaker et al., 2019;Zelnik et al., 2018). In considering more realistic environmental conditions and their interaction with population dynamics in this way, predictions are likely to be of greater utility (Shoemaker et al., 2019). For example, the temporal and spatial distribution of species at a given location can be revealed through integration of environmental fluctuations based on observed or projected climate data. While this is straightforward in principle, the majority of experiments manipulating climate variables misrepresent local climate projections (Korell et al., 2019). Despite this mismatch, the potential of locally informed environmental variables in experiments and simulations is clear (e.g., Fry et al., 2013).
Finer resolution population dynamics can also enable predictions of greater potential accuracy. Specifically, individual based models-in which each individual is characterized by a different phenotype-can predict whether particular species or phenotypes will succeed in establishing in new habitats. In this case, both the initial location of individuals in space, and initial environmental conditions will dramatically affect establishment success. Individuals or populations at their range margins have a greater probability of invading new habitats (Hallatschek & Nelson, 2008;Waters et al., 2013). However, the successful establishment of individuals at the invasion front hinges on both initial environmental conditions and the temporal dynamics of the environment, particularly when fitness depends on environmental conditions (leaving aside phenotypic evolution, see Kubisch et al., 2014, Norberg et al., 2012. Both spatial and temporal scale can be clearly integrated into studies of range dynamics to better understand the consequences of spatiotemporal scale for range shifts and contractions in light of environmental suitability (e.g., Keith et al., 2008;Sarmento Cabral et al., 2013), and for expansion velocity in the expanding population. The likelihood of successful population expansion differs across spatiotemporal scales for a given population (Villa Martín et al., 2019). At small scales, for which environmental change is less observable, homogeneous populations are predicted to be more successful. However, when considering a larger spatial or temporal extent, phenotypic heterogeneity will improve the performance of an invading population, since diverse phenotypes have a greater range of responses to environmental conditions and so buffer populations against environmental change (Elmqvist et al., 2003). Outstanding questions then are to what degree enlarging the spatial (or temporal) extent of the study increases the advantageousness of population heterogeneity, and whether there is a clear threshold in scale beyond which phenotypic heterogeneity clearly outperforms homogeneous populations as a strategy for invasion success. Clearly, disentangling the causes and consequences of species invasions and population expansions or contractions involves considering spatiotemporal scale as well as the intrinsic and extrinsic mechanisms governing population dynamics.

| INTERACTING SCALE-DEPENDENCIES ACROSS LEVELS OF ORGANIZATION
Thus far, we have discussed the sources of variation in ecological dynamics (both intrinsic and extrinsic) in isolation, and some pathways through which spatial or temporal scale can influence understanding of these dynamics (Figure 4). Yet these sources of variation do not act in isolation to generate the dynamics inherent to natural systems. In reality, intrinsic and extrinsic drivers of variation act in concert across levels of ecological organization (Leibold et al., 2004). For example, studies of deterministic community assembly recognize that the non-neutral processes structuring ecological communities include both biotic and abiotic factors (Gravel et al., 2011;Vellend, 2016). The spatial scale of consideration then feeds back to influence the primary mechanism structuring communities. For example, environmental constraints on functional and phylogenetic diversity of island birds may be relaxed on larger islands (Ross et al., 2019), and biotic homogenization may be more detectable at larger spatial scales when considering biodiversity change . Sources of variation at one organizational level also can propagate across other levels. For example, the intrinsic demographic processes acting on a metacommunity hub or connector species can modulate metacommunity structure (Toju et al., 2017). However, the paucity of studies investigating ecological dynamics or stability at more than a single scale F I G U R E 4 Overview of our discussion. Figure illustrating the methods and sources of variation we discussed with respect to ecological stability and dynamics. Different colored circles represent different species. At small spatial scales, interaction variability is an intrinsic source of variation affecting community network structure. Considering larger spatial scales, compositional variability is an intrinsic source of variation influencing community assembly. At macroecological scales, species' ranges are influenced by demographic stochasticity (intrinsic variation) and its interaction with environmental variability (extrinsic variation) in space and/or time. We discussed how temporal scale influences community network structure inferred through empirical dynamic modeling, how spatial scale influences community assembly dynamics measured via energy landscape analysis and how both time and space influence species' ranges through environmental heterogeneity. In the final section, we discuss how these sources of variation may interact with each other and with spatiotemporal scale to add further complexity to our understanding of ecological stability and dynamics [Color figure can be viewed at wileyonlinelibrary.com] or organizational level makes it difficult to resolve the interdependencies of ecological dynamics on sources of variation across levels of organization .
Better integration of extrinsic variables into models of intrinsic sources of variation may allow us to tease apart the effects of environmental from intrinsic factors governing population or community dynamics (Bocedi et al., 2014;Harfoot et al., 2014;Pagel & Schurr, 2012). The approaches described above are each amenable to the inclusion of multiple sources of variation. For example, when reconstructing the state space of interspecific interactions using empirical dynamic modeling, the dynamics of these interactions are affected not only by intrinsic changes to the abundance of each species, but also by the environmental conditions in which these interactions are embedded. Munch et al. (2019) discuss how the state space reconstruction may be influenced by both direct and indirect effects on the species, and provide some guidance on interpreting results in light of this. Considering community assembly pathways, energy landscape analysis can explicitly model environmental conditions when defining the site-by-species community matrix, assigning a value representing local environmental conditions at each site ( Figure 2). The analysis then produces a separate matrix of explicit abiotic factors and models the relationship between the environment and a given community state (Figure 3). Perhaps the approach with the clearest potential for considering multiple sources of variation is to explicitly define model parameters including demographic rates and environmental suitability in population models as in Villa Martín et al. (2019). Such models can elegantly consider the multiple sources of intrinsic and extrinsic variation responsible for shaping population dynamics at species' range margins (Bocedi et al., 2014;Pagel & Schurr, 2012;Singer et al., 2016).
Disturbances are inherently scale dependent (Pinek et al., 2020;Ryo et al., 2019;Zelnik et al., 2018). Yet, responses to disturbance are also a function of the spatial or temporal scale of consideration Jacquet & Altermatt, 2020;Moreno-Mateos et al., 2020;Waldock et al., 2018). Accordingly, scale underpins our understanding of the responses of ecological communities to disturbance, including apprehension of dynamics and stability across levels of ecological organization. Increasing recognition of the scale-dependent nature of most, if not all, studies of ecological dynamics and stability will be critical (Gonzalez et al., 2020), particularly since there is presently little we can do to negate the impact of scale in empirical studies. When considering scaling in ecology, ecosystems are mostly assumed to be in equilibrium, but this may be the exception rather than the rule (Newman, 2019). Nonequilibrium dynamics add another layer of complexity to the consideration of ecological systems, as ongoing alterations to disturbance regimes are systematically driving biodiversity change across the globe (McGill et al., 2015;Newman, 2019;Wolkovich et al., 2014). Given that our understanding of ecological stability is built on assumptions of scale, it is then imperative that we better recognize and begin to account for scale-dependence in the observation of ecological dynamics and stability under global environmental change. Doing so will facilitate informed decisions about conservation policy and the successful management of nature's contributions to people.

ACKNOWLEDGMENTS
This contribution arose from discussions at the Ecological Research Symposium "Ecological Stability: spatial and temporal dynamics" organized by Samuel R. P.-J. Ross and Yuka Suzuki at the 66th meeting of the Ecological Society of Japan in Kobe, Japan. We thank Masahiro Nakaoka, Katsuhiko Yoshida and the Ecological Society of Japan's ER symposium committee for facilitating discussions through our symposium. Funding to attend the conference was provided for Maria Dornelas by Ecological Research, and for Samuel R. P.-J. Ross by a British Ecological Society Training and Travel Grant (TT19/1029). Samuel R. P.-J. Ross was supported by Trinity College Dublin and an Irish Research Council Postgraduate Scholarship (GOIPG/2018/3023). Yuka Suzuki was supported by JSPS KAKENHI (JP20J10699) with additional support to Yuka Suzuki and Paula Villa Martín from the Okinawa Institute of Science and Technology Graduate University. Michio Kondoh was supported by CREST (JPMJCR13A2), Japan Science and Technology Agency and JSPS KAKENHI (19H05641 and 16H04846). Kenta Suzuki was supported by the Management Expenses Grant for RIKEN BioResource Research Center, MEXT, and JSPS KAKENHI (20K06820 and 20H03010). We also thank three anonymous reviewers for their positive and constructive comments on earlier versions of this manuscript.