The addition of spatial structure to ecological concepts and theories has spurred integration between sub-disciplines within ecology, including community and ecosystem ecology. However, the complexity of spatial models limits their implementation to idealized, regular landscapes. We present a model meta-ecosystem with finite and irregular spatial structure consisting of local nutrient–autotrophs–herbivores ecosystems connected through spatial flows of materials and organisms. We study the effect of spatial flows on stability and ecosystem functions, and provide simple metrics of connectivity that can predict these effects. Our results show that high rates of nutrient and herbivore movement can destabilize local ecosystem dynamics, leading to spatially heterogeneous equilibria or oscillations across the meta-ecosystem, with generally increased meta-ecosystem primary and secondary production. However, the onset and the spatial scale of these emergent dynamics depend heavily on the spatial structure of the meta-ecosystem and on the relative movement rate of the autotrophs. We show how this strong dependence on finite spatial structure eludes commonly used metrics of connectivity, but can be predicted by the eigenvalues and eigenvectors of the connectivity matrix that describe the spatial structure and scale. Our study indicates the need to consider finite-size ecosystems in meta-ecosystem theory.
The concepts of the population, the community  and the ecosystem  are fundamental to ecological understanding. In order to operationalize these concepts into usable components of theory, ecologists have added temporal, genetic and spatial structure to the concepts. The successful incorporation of space has led to metapopulation , metacommunity  and meta-ecosystem  theories, which have in turn spurred new experiments, observations and models and have renewed hope for the integration of community and ecosystem ecology through spatial structure [6–8].
Modifying population, community or ecosystem models to include space has been done through limiting cases such as two connected systems [9,10], implicit space [11,12] or infinite continuous [13,14] or discrete spatial domains . These limiting cases can help us analyse spatial processes by simplifying the spatial structure. For example, analysing the effects of the movement of nutrients on ecosystem dynamics and functioning in meta-ecosystems was simplified by using idealized spatial structures [9,10,16].
However, the spatial structure of ecological systems are finite and irregular [17–20] and these finite, irregular features of the physical landscape can affect how organisms and materials are distributed across space by affecting both interactions and movement [6,17,20]. Population persistence and dynamics are affected by realistic landscapes in metapopulation [21–25], epidemic  and predator–prey  models, and metrics capturing these effects of the landscape [21–25] on population persistence have been derived and related to landscape connectivity [21,28]. What is lacking are equivalent metrics capturing the effects of spatial structure on ecosystem dynamics and functioning, which are likely to be affected by the different movement rates of organisms and materials.
The goal of our study is to fill this gap by expanding meta-ecosystem theory to include finite landscapes and to examine how the movement of organisms and materials interacts with landscape connectivity to impact the stability, dynamics and functioning of ecosystems. We do so by creating a meta-ecosystem model that consists of nutrient–autotroph–herbivore ecosystems that exchange materials and organisms and has a spatial structure determined by a finite spatial network that mimics aspects of real landscapes.
Our results show that high nutrient and high herbivore movement rates can destabilize the meta-ecosystem and can lead to spatially heterogeneous dynamics, but the destabilization and its associated dynamics are dependent on spatial structure and the autotroph movement. The effect of spatial structure is revealed by the ‘scales of spatial interactions’ (non-zero eigenvalues of the connectivity matrix) that emerge from the differences in connectivity between ecosystems, and which scale is associated with the destabilization can predict the dynamics seen in the meta-ecosystem. For example, the dynamics associated with small scales of spatial interaction have large oscillations in highly connected ecosystems, and smaller oscillations in less connected ecosystems. Furthermore, our analysis reveals how the scales of spatial interactions cannot be easily explained through other network connectivity metrics. In addition, the spatial structure of a meta-ecosystem can affect its primary and secondary productivity, indicating complex effects of spatial structure on ecosystem function. Our results provide new ways of integrating finite spatial structure into meta-ecosystem theories and of interpreting its impact on ecosystem function.
2. Material and methods: the meta-ecosystem model
(a) Regional and local processes in meta-ecosystems
We modify a meta-ecosystem model  to highlight the effects of spatial structure on meta-ecosystem dynamics and functioning. The model can be broken up into regional processes that connect ecosystems and local processes, which describe the internal dynamics of the ecosystems (figure 1). The regional processes are the movements of materials and organisms between the ecosystems, whereas the local processes are trophic and non-trophic interactions (nutrient recycling; figure 1). As in previous work , we limit our local ecosystems to one limiting nutrient R, and we track the stocks of that nutrient in autotrophs A and herbivores H (figure 1).
The movements of materials and organisms between ecosystems are determined by the connectivity matrix (C) and the movement matrix (D), which allows us to separate out the effect of spatial structure from the effects of movement rates . The connectivity matrix is an n × n matrix, where n is the number of ecosystems in the meta-ecosystem, whose off-diagonal entries (i.e. cij, i ≠ j) indicate the links between the different ecosystems. The diagonal entries (i.e. cii) represent the total number of connections that ecosystem i has, normalized by the total possible connections it could have. Because of certain beneficial properties for analysis (see the electronic supplementary material, appendix A), we will consider only symmetric connectivity matrices (i.e. cij = cji) that allow for no loss of materials and organisms during movement between ecosystems (i.e. ).
The movement matrix is a k × k matrix, where k is the number of ecosystem compartments in each ecosystem. The diagonal entries of the movement matrix are the movement rates of each compartment, whereas the off-diagonal entries would indicate cross-movement, which occurs when the movement of one ecosystem compartment (e.g. autotrophs) is dependent on another compartment (e.g. herbivores). However, we will not consider cross-movement in this study, making the movement matrix a 3 × 3 diagonal matrix with entries dR, dA and dH, which are the movement rates of the limiting nutrient, the autotrophs and the herbivores, respectively.
At the local level, we have the available limiting nutrient R at the base of the ecosystem. It is supplied at a constant rate I from rock weathering and other abiotic sources and is lost at rate E proportional to its concentration in the medium (e.g. soil, water). Part of the available limiting nutrient R is assimilated into the autotrophs based on their uptake function U(R, A), but this is balanced by nutrient recycling from autotroph losses herbivore losses χL(H) and assimilation inefficiencies from herbivory γW(A, H). Autotroph nutrient stocks increase through the uptake function U(R, A), but decrease through their intrinsic losses M(A) and herbivory W(A, H). Herbivore nutrient stocks increase through herbivory (1 − γ)W(A, H) and decrease through intrinsic losses L(H).
Combining regional and local processes gives us the following system of ordinary differential equations: 2.1a 2.1band 2.1cFor analytic simplicity, we assume that the parameters within the functions are the same across the meta-ecosystem (but see the electronic supplementary material, appendix). While much of our mathematical analysis can be done with the functions of equation (2.1) (see the electronic supplementary material, appendix B), our numerical simulations require specified functions. We assume type II/Michaelis–Menten functional responses based on their widespread prevalence in plants  and herbivores  for uptake, and density-independent losses for both the autotrophs and herbivores, which transform equation (2.1) into 2.2a 2.2band 2.2cwhere αA and αH are maximum uptake rates, βA and βH are the half-saturation constants and m and l are density-independent loss rates.
(b) Effects of spatial structure on meta-ecosystem stability
The meta-ecosystem model presented above can exhibit a wide range of dynamical behaviour, even when no movement occurs . However, our focus is on how meta-ecosystem stability is modulated by its spatial network structure. As in Marleau et al. , we use movement parameters to perturb the meta-ecosystem, allowing us to highlight spatial processes instead of local processes.
For simplicity, we restrict ourselves to parameter ranges that allow a unique stable equilibrium when there are no regional processes (but see the electronic supplementary material, appendix). In other words, each ecosystem in the meta-ecosystem will be in the same state if no nutrients and no organisms are moving between them. If we add regional processes, each ecosystem will return to this state after perturbations as long as the Jacobian matrix describing the dynamics of the whole meta-ecosystem at the spatially homogeneous equilibrium state has only eigenvalues with negative real parts. This Jacobian matrix, which would normally be a 3n × 3n matrix, can be broken up into n matrices of the following form : 2.3where J is the Jacobian matrix of a local ecosystem without regional processes and λi is the ith eigenvalue of the connectivity matrix C (see the electronic supplementary material, appendix B). If each V(i) matrix has eigenvalues, ϕik, with negative real parts, then each ecosystem will return to its original equilibrium state after perturbation.
The stability of the ecosystem equilibrium is lost when the real part of at least one eigenvalue ϕik of one of the V(i) matrices becomes positive (bifurcates) as the movement rates of the nutrients and organisms change. Such bifurcations can lead to: (i) spatially heterogeneous equilibria (see the electronic supplementary material, appendix B), with individual ecosystems at different equilibrium values, which can include some local ecosystems having no autotrophs, or (ii) spatially heterogeneous oscillations across individual ecosystems . The minimum, positive movement rates necessary to cause the bifurcations are defined as the minimum critical movement rates, where X = R, H, as autotroph movement cannot cause a bifurcation. The critical minimum herbivore movement rate, is associated with spatially heterogeneous equilibria, while the critical nutrient movement rate, is associated with spatially heterogeneous oscillations. Our analysis will focus on the spatially heterogeneous oscillations to compare with previous work .
(c) Scales of spatial interaction in meta-ecosystems
The non-zero eigenvalues of the connectivity matrix play an important role in creating spatial heterogeneity (equation (2.3)), similar to how the dominant eigenvalue (λ1 = 0) of any connectivity matrix represents the dynamics without spatial structure . These eigenvalues represent how the spatial structure of the meta-ecosystem influences the response of the ecosystems to perturbations. In equilibrium contexts, they are the equivalent of the wavenumber or spatial frequency in reaction–diffusion models, which indicates the spatial scale of perturbations to the equilibrium . This spatial scale of interaction between the ecosystems is thus key in predicting the ability of a perturbation to propagate across a meta-ecosystem through movement of individuals and matter .
Therefore, we define a scale of spatial interaction to be a non-zero eigenvalue λi of the connectivity matrix C, with more negative eigenvalues representing smaller spatial scales. For convenience, we order the non-zero eigenvalues from the largest (i.e. less negative) to the smallest (i.e. most negative) such that such that we go from large scales of spatial interaction to the smallest. Each meta-ecosystem could have multiple, unique (up to n − 1, excluding λ1) scales of spatial interaction, each corresponding to an eigenvalue that lies between 0 and −2 (see the electronic supplementary material, appendix A).
Each scale of spatial interaction has an associated eigenvector that can provide information on the amplitudes and frequencies of the emergent spatial perturbation. For finite meta-ecosystems observed here, the eigenvectors would be associated with the oscillations within the individual ecosystems. The signs of the elements of the eigenvector indicate the synchrony between ecosystems, and their relative magnitudes represent the amplitudes of fluctuations in an ecosystem (see the electronic supplementary material, appendix C). Combined, these tools can help us examine the temporal and spatial stability of finite-size meta-ecosystems, though such information is of limited value for scales of spatial interaction that are repeated as the associated eigenvectors can offer different predictions.
Previous work on metapopulations showed the effect of the dominant eigenvalue of the connectivity matrix, which should represent the shortest scale of spatial interaction (λn: [23–25]). Here, we examine the general relationship between the scales of spatial interaction and our critical movement rates to discern which eigenvalues of the connectivity matrix play a role in determining meta-ecosystem stability.
(d) Meta-ecosystem dynamics and functioning
We use our general model (equation (2.1)) to quantify the critical minimum movement rates that capture the shift from spatial homogeneity to spatial heterogeneity. We also use numerical simulations of our specified model equations (2.2) to detail the dynamics of individual ecosystems within the meta-ecosystem across critical movement rates for stability. We illustrate our results with 5 × 5 connectivity matrices with equal link density, but differing spatial network structure. In addition, we examined how the spatial scale of interaction associated with the spatial heterogeneity can determine meta-ecosystem dynamics by altering the movement rate of the autotrophs.
The implication of dynamical responses to movement for meta-ecosystem function is assessed by measuring the average primary and secondary production at the meta-ecosystem level for increasing rates of nutrient movement (dR) in meta-ecosystems with differing network structures. Average primary and secondary production was measured by evaluating U and (1 − γ)W over 5000 time steps, respectively. In addition, we compare the rankings of network structures for primary and secondary production across scales of spatial interactions causing meta-ecosystem destabilization.
(e) Relating network connectivity properties to scales of spatial interaction
Landscape ecology and network theories have produced a number of metrics to characterize connectivity. We use our model to determine whether the scales of spatial interaction are related to two common metrics of connectivity associated with network stability: link density and maximum node degree. There are other metrics available , but we focus on these metrics in order to capture meta-ecosystem-level connectivity with a single number for use in prediction.
Link density is the number of links divided by the number of nodes in the network, which in our case is the number of ecosystems in the meta-ecosystem. Maximum node degree is the number of links found at the most connected node in the network. For our purposes, we derived a relative scale of maximum node degree that goes from 0 (minimum) to 1 (maximum; electronic supplementary material, appendix A).
We used 694 071 randomly generated 30 × 30 connectivity matrices to discern whether any or all of these connectivity measures can predict the scales of spatial interaction and can provide a link between meta-ecosystem properties and connectivity (see the electronic supplementary material, appendix A).
(a) Meta-ecosystem stability: the interaction between movement and spatial structure
We first analyse the stability of the spatially homogeneous meta-ecosystem following its perturbation by the movement of nutrients and other organisms. We derive functions of critical movement rates corresponding to a transition from a spatially homogeneous ecosystem state throughout the meta-ecosystem to one with significant spatial heterogeneity (see the electronic supplementary material, appendix B) 3.1aand
where jlk is the row l and column k element of Jacobian matrix J, det(J) is the determinant of the Jacobian matrix, and B and C are complex functions of movement rates and of eigenvalues of the connectivity matrix (see the electronic supplementary material, appendix B). From these functions, we derive the minimum movement rates required to create spatial heterogeneity for given movement rates and a given connectivity matrix 3.2aand 3.2bwhere n is number of patches in the meta-ecosystem, which means that there is only a finite number (n − 1) of non-zero eigenvalues for a specific connectivity matrix. Furthermore, the number of unique non-zero eigenvalues can range from 1 to n − 1, which means that few scales of spatial interaction (i.e. few unique λi) can be present even in large n meta-ecosystems.
The relationship between the scales of spatial interaction (λi) and the minimum critical movement rates ( and ) depends strongly on the movement rate of autotrophs (dA; figure 2). When dA is low, and decrease with decreasing λi (figure 2a,b). Therefore, for a given metaecosystem, the value of λn (the shortest scale of spatial interaction) determines the minimum critical movement rates at low dA (figure 2a,b). As meta-ecosystems with high maximum node degree and high link density have smaller λn, meta-ecosystems with greater connectivity, according to those metrics, are more easily destabilized by nutrient and herbivore movement at low dA (figure 3a,b).
At higher dA values, and show a parabolic relationship with λi, such that minima are reached at intermediate values of λi (figure 2c,d). For a given meta-ecosystem, this can lead to either a longer scale of spatial interaction (e.g. λn−1) determining the minimum critical movement rates or it can result in no destabilization being possible as the λi's all lead to negative and values (figure 2c,d). In other words, meta-ecosystems lacking longer scales of spatial interaction will not be destabilized at high dA. Therefore, all the scales of spatial interaction would need to be evaluated to determine the stability, not just the shortest (figure 2c,d). However, network connectivity metrics provide little guidance in predicting what scales of spatial interaction to expect at given connectivity levels and hence provide little help in determining meta-ecosystem stability (figure 3c,d).
(b) Meta-ecosystem dynamics: dependence on scale of spatial interaction and spatial structure
The realized dynamics of local ecosystems after meta-ecosystem destabilization depend on their spatial structure and the scales of spatial interaction (figure 4; see the electronic supplementary material, appendix D for parameter values). If the destabilization is associated with the shortest scale of spatial interaction, λn, the ecosystems with higher node degree (and with neighbours with higher node degree) within the meta-ecosystem have greater amplitude oscillations than ecosystems with lower node degree (figure 4a,b). As the node degree of each ecosystem depends directly on network structure, meta-ecosystems with different network properties display different synchrony patterns between ecosystems (figure 4a,b). In particular, ecosystems with the same connectivity properties (e.g. node degree and node degree of immediate neighbours) exhibit synchronized and identical oscillations (figure 4a,b). These dynamics can be seen even with spatially heterogeneous nutrient supplies or high nutrient supplies that lead to local oscillations without movement (see the electronic supplementary material, appendix E).
When the scales of spatial interaction between ecosystems are longer (e.g. λn−1), the resulting dynamics are more complex (figure 4c,d). For example, it is possible for central ecosystems within the network to be barely affected by the instability, whereas the outer ecosystems show large and anti-phase oscillations (figure 4c). Or there can be little discernible pattern in the spatial distribution of oscillations across the meta-ecosystem (figure 4d). Furthermore, ecosystems with similar connectivity properties no longer demonstrate synchronized dynamics nor do they necessarily oscillate at the same amplitude as they did at shorter scales of spatial interaction (figure 4c,d).
The dynamics shown above can, in certain cases, be predicted through the eigenvectors associated with the scales of spatial interaction (see the electronic supplementary material, appendix D). For example, meta-ecosystems with the same scale of spatial interaction and the same associated eigenvector will have the same spatial and temporal dynamics close to (see the electronic supplementary material, figure C.1). However, the predictive ability of the eigenvectors weaken as the movement rates increase beyond . The reason for this is that the other scales of spatial interaction could destabilize the meta-ecosystem at the new higher rates of movement independently of the original destabilizing scale and their contributions to the dynamics become significant (see the electronic supplementary material, appendix C).
(c) Meta-ecosystem production
The differences in meta-ecosystem dynamics owing to spatial structure lead to differences in meta-ecosystem functioning (figure 5). When small-scale interactions cause the spatially homogeneous equilibrium to lose stability, both primary and secondary production generally increase with increasing nutrient movement (figure 5a,b). The specific network structures show differences in terms of their production, with the network with the highest maximum node degree consistently having the highest production at all levels of nutrient movement rate (figure 5a,b).
Similar to the results involving meta-ecosystem dynamics, the destabilization associated with large-scale interactions results in patterns in meta-ecosystem production that differ from those associated with small-scale interactions (figure 5c,d). Both primary and secondary production show small, non-monotonic increases with increasing nutrient movement relative to the equilibrium case (figure 5c,d). Furthermore, the network with the highest maximum node degree consistently has the lowest primary production, which is in opposition to the small-scale spatial interaction case, though it does not always hold for secondary production (figure 5c,d).
Our analysis reveals that non-zero eigenvalues of the connectivity matrix describing the finite spatial structure of the meta-ecosystem, and hence multiple scales of spatial interaction, can determine meta-ecosystem stability and productivity in response to increasing movement of matter and organisms. We also show how the scales of spatial interactions driving the loss of meta-ecosystem stability can predict the distribution of dynamical regimes among local ecosystems. The study of finite-size landscapes escapes predictions relating ecosystem dynamics to the dominant scales of spatial interaction and made under the assumption of infinite or well-mixed space. Our results demonstrate that the importance of multiple scales of interactions for resolving the complex response of stability and productivity to fluxes of matter and individuals across meta-ecosystems of finite size.
(a) Scales of spatial interaction: the importance of finite space in ecological models
Our study uses the non-zero eigenvalues and eigenvectors of the connectivity matrix to help characterize the spatial structure of meta-ecosystems. These eigenvalues represent the scales of spatial interaction between local ecosystems, indicate how the ecosystems will respond to perturbations and each of them can lead to the destabilization of meta-ecosystesms. In addition, these eigenvalues are not well predicted by other measures of spatial structure commonly used in ecology, which makes them a novel tool for research [17,20]. Furthermore, when modelling multi-level movement models, the need to examine all scales of spatial interaction leads to the inadequacy of two-patch  and infinite domain  models to capture these critical elements of spatial structure.
Two-patch models can be represented by a connectivity matrix with a single eigenvalue or scale of spatial interaction . This makes the two-patch model similar to a fully connected network as every ecosystem is connected to every possible neighbour and its results do not scale up when networks are not fully connected. By contrast, models of infinite domain can be formulated to contain all possible scales of spatial interaction such that movement could always destabilize them if destabilization is possible . The issue here is that metapopulations, metacommunities and meta-ecosystems are not infinite in size, but finite [4,8,21,22,27]. Such finite size leads to a limited number of associated scales of spatial interaction where perturbations can manifest and destabilize meta-ecosystems at equilibrium. Therefore, infinite domain models may predict that a perturbation will destabilize an ecological system, while a more realistic finite network model will predict its stability.
(b) Unequal movement rates of materials and of organisms across landscapes
Spatial instabilities are produced and their impacts on dynamics and functioning are modulated through the unequal movement of nutrients and of organisms [29,32,33]. Our study shows how multiple scales of spatial interaction resulting from such movement are determined by the spatial structure of the meta-ecosystem. Our model predicts that if the shortest scale of spatial interaction drives destabilization, increasing connectivity (i.e. the link density or the maximum node degree of the meta-ecosystem) can further promote instabilities. However, increasing the movement rate of the autotrophs allows for longer scales of spatial interaction to have a dominant destabilizing role, which then results in the relationship between instabilities and connectivity to be highly irregular.
Interactions between the structure and rates of movement have profound implications for ecosystem functioning. The loss of regional stability driven by spatial structure and movement rates leads to the emergence of source-sink ecosystems for autotrophs [10,34] and to increased productivity at the regional level. However, the production gains are offset by increases in variability at local and regional scales, which may lead to local loss of autotrophs and herbivores . Our results support the importance of integrating movement properties of ecosystem compartments to the spatial structure of landscapes.
Our analysis emphasizes the importance of unequal movement across ecosystem compartments. Other efforts to discern the importance of movement on community and ecosystem processes have postulated that spatial structure could be subsumed through the coupling of fast moving, high trophic-level organisms, which would have the strongest support in marine ecosystems [12,35,36]. However, coupling can occur at all trophic levels, and we predict that even limited movement can have large impacts on dynamics, functioning and stability. The study of differential movement in finite-size ecosystems should allow for greater integration between food web and landscape ecology [5,6,8].
Our model has several limitations with regards to the unequal movement rates of materials and of organisms that should be addressed in future studies. First, the impacts of cross-diffusion were ignored, even though herbivores can serve as vectors for autotroph movement . Second, adding another autotroph or herbivore with a different movement rate should be explored in order to determine how this can impact dynamics, functioning and stability. Lastly, we recognize that ecosystems are not immobile spatial patches, but are spatially distributed and formed through the interactions between their biotic and abiotic components [8,38]. Given these features of natural ecosystems, our study reinforces the notion that more research is needed to understand the impacts of connectivity on communities and ecosystems if we are to develop better conservation strategies [39,40]. Our study points towards the multiple scales over which ecosystems interact across landscapes as a tool to predict and understand their complex dynamics.
J.N.M. was supported by a NSERC PGS-D2 and FQRNT B2. F.G. was supported by the Natural Science and Engineering Research Council (NSERC) of Canada. M.L. was supported by the TULIP Laboratory of Excellence (ANR-10-LABX-41).
- Received August 12, 2013.
- Accepted December 4, 2013.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.