A Generalized Spatial Measure for Resilience of Microbial Systems.

The emergent property of resilience is the ability of a system to return to an original state after a disturbance. Resilience may be used as an early warning system for significant or irreversible community transition; that is, a community with diminishing or low resilience may be close to catastrophic shift in function or an irreversible collapse. Typically, resilience is quantified using recovery time, which may be difficult or impossible to directly measure in microbial systems. A recent study in the literature showed that under certain conditions, a set of spatial-based metrics termed recovery length, can be correlated to recovery time, and thus may be a reasonable alternative measure of resilience. However, this spatial metric of resilience is limited to use for step-change perturbations. Building upon the concept of recovery length, we propose a more general form of the spatial metric of resilience that can be applied to any shape of perturbation profiles (for example, either sharp or smooth gradients). We termed this new spatial measure "perturbation-adjusted spatial metric of resilience" (PASMORE). We demonstrate the applicability of the proposed metric using a mathematical model of a microbial mat.

Introduction

Complex networks of interacting components produce outcomes that cannot be easily predicted, even when the state of the network components and the inputs to the network are known. These difficult-to-determine outcomes have come to be known as emergent phenomena or higher-order properties, which emerge from the functioning of the whole network, rather than as a simple sum of the individual states of the parts (De La Fuente et al., 2008; Schubert, 2014). Emergent properties arise in complex networks that impact our lives, such as in the communities of microbes and higher organisms that compose ecosystems (Tilman et al., 2001), social networks (Lusseau, 2003; Mitrovic and Tadic, 2010), military-political structures (Porter et al., 2005), commercial systems (Carey and Carville, 2003; Cimellaro et al., 2010), and climate and weather systems (Higgins et al., 2002; Easterling and Kok, 2003). Consequently, developing the ability to understand and predict emergent properties from complex networks is of great importance.

Microorganisms are commonly found physically associated with one another in spatially structured communities such as biofilms or microbial mats (Curtis and Sloan, 2004; Wagner et al., 2006; Fuhrman, 2009; Robinson et al., 2010; Renslow et al., 2011). These communities may range from monocultures to highly diverse assemblages of species; however, in either case, individual members operate and interact in ways governed by their individual functional responses and local microenvironmental conditions (Schramm et al., 1996; Hibiya et al., 2003; Bernstein et al., 2013, 2014). Formation of a biofilm matrix composed of extracellular polymeric substances confers significant fitness advantages on the microbes it shelters, including physical protection, such as from predation or shearing forces, reduction of environmental stresses, such as from rapid changes in environmental conditions or exposure to antibiotics, facilitation of beneficial interspecies relationships, and rapid exchange of genetic material (Flemming and Wingender, 2001; Laspidou and Rittmann, 2002; Czaczyk and Myszka, 2007; Cao et al., 2011). It is likely that emergent properties arise from the spatial organization of microbes, forming microenvironments and promoting interconnectedness of a multi-species metabolic network through resource exchange and intercellular communication (Xavier and Foster, 2007; Konopka, 2009; Wintermute and Silver, 2010).

Resilience is a higher-order property in microbial communities, characterized by the ability to recover from a perturbation or disturbance (Allison and Martiny, 2008; Shade et al., 2012; Griffiths and Philippot, 2013; Hawkes and Keitt, 2015). While resilience is not yet conclusively defined in microbial communities, we use this term to imply the rate of recovery of a given function in a community (or, more generally, their functional relationship with environmental variables) after perturbation. The concept of functional resilience in both engineered and ecological microbial systems, and the relationship between state and functional properties is examined in detail in our companion paper (Song et al., 2015). Attempts at quantifying resilience have primarily been done by monitoring functional recovery over time, with resilience being negatively correlated to recovery time (i.e., the time required for function to recover a defined percent of original function), or being correlated to recovery speed (i.e., the rate of recovery of function with units of slope). Faster recovery of function is therefore an indicator of higher resilience. These ideas have been explored theoretically and experimentally (Wissel, 1984). For example, recovery rate was monitored for cyanobacterial cultures exposed to a dilution perturbation by flushing out 10% of the population volume; the recovery rate decreased as the population lost resilience and approached a tipping point where function became irrecoverable (Veraart et al., 2012). In this experiment where the recovery over time was readily observable over a relevant timescale, assessing the system’s resilience is straightforward. Such quantification of recovery time is important because it is known that systems nearing the verge of collapse or those trending toward unstable dynamics frequently exhibit reduced rates of recovery, or decreasing resilience (van Nes and Scheffer, 2007). This is known as “critical slowing down,” meaning that, as a system nears a tipping point, recovery of function after a perturbation slackens (Scheffer et al., 2009, 2012). Thus quantitative measures of resilience are essential for early warning of impending system collapse, from which function cannot be regained.

Although recovery rate provides a direct measurement of resilience, in some systems it is not possible or practical to measure recovery time on time scales relevant to the perturbation. Additionally, some systems do not allow for measurement of resilience because the experiments needed for quantify recovery time are impossible, unethical (e.g., when imposing a perturbation and monitoring function would endanger ecosystems or humans), or otherwise detrimental to the function of neighboring systems. In such cases not amenable to temporal analysis, resilience can alternatively be quantified in terms of spatial recovery as a proxy for temporal recovery (Dai et al., 2013). This means that resilience may still be predicted by observing spatial, rather than temporal, features of the system in cases where temporal and spatial recovery are reasonably correlated. The ability to monitor microbial function over relevant spatial scales has become increasingly possible due to advances in imaging and sensor capabilities, enabling a link between function, structure, and microbial identity (Behrens et al., 2008; Li et al., 2008; Pett-Ridge and Weber, 2012; Babauta et al., 2014; Vanwonterghem et al., 2014).

Recently, Dai et al. (2013) put forward recovery length as a measure for resilience. Recovery length is characterized by the distance required for function to recover from a spatial perturbation, and it was initially based on the observation that recovery length, which was correlated to recovery time, increased when a yeast system was on the verge of collapse. In their experiment, populations of Saccharomyces cerevisiae were connected spatially along a one-dimensional array through discrete dispersal events. Population stability was measured as the dilution factor was increased, imposing a perturbation. Recovery length was quantified by measuring the population density across distance, which provided a warning signal of imminent population collapse as the dilution factor approached unsustainable levels. Furthermore, these indicators increased with gradual increases in dilution factor, revealing deterioration of resilience in the system in real time. The recovery length concept proposed by Dai and coworkers faithfully quantified resilience to perturbations with sharp boundaries (such as step-change perturbations). As the authors pointed out, however, the same may not be effective for more realistic forms of perturbations that follow gradients or which are blurred across the spatial dimension. In reality, perturbation may come as gradients across both time and space. For example, spatial perturbations that form gradients include the intrusion of substrate or toxins into soil and its subsequent removal (Cao et al., 2012; Moran et al., 2014) or sunlight penetrating into a microbial mat (Bhaya et al., 2001; Häder and Lebert, 2001; Lindemann et al., 2013).

In this study we propose a perturbation-independent resilience metric, termed perturbation-adjusted spatial metric of resilience (PASMORE), to be applicable to perturbations that may be either step-change or gradient-based, and thus extending the recovery length work by Dai et al. (2013). To test this new metric, we developed a simple model of a microbial mat containing species with differing resilience. We discuss cases where the previously proposed metric fails to faithfully measure resilience under gradient perturbations, and then demonstrate how PASMORE provides an appropriate assessment of resilience for the same and other types of perturbation profiles. We also explored the impact of limited knowledge about the shape of perturbation profiles on the effectiveness of PASMORE. Finally, we investigated a test application of PASMORE to systems of two interacting species, and observed how resilience is affected by interspecies competition.

Model System

We considered a simple mathematical model, which simulates a microbial mat with populations of microbial species: one with faster motility, and another with slower motility (see Figure for illustration). Each species responds to a perturbation (e.g., light) with movement proportional and opposite to the perturbation; during perturbation, both species move away from the perturbation toward a region where the perturbation is not present. The key aspect of this simulation is that the motility allows the microbial species to recover their initial population density after the perturbation occurs, with the more motile species recovering more quickly, and thus displaying higher resilience. While resilience is bestowed by motility rate in our model system, there are many other properties of a microbial species that can confer resilience beyond motility, and these features may present themselves within spatial constructions besides microbial mats. Furthermore, our proposed resilience measure has been formulated to be agnostic to the type of function being monitored. As discussed in Song et al. (2015), it is the researcher’s job to identify the function of interest that is measured in a spatially defined community. Function may be defined as any trait or variable of interest, and resilience of the chosen function must be relative to a specific perturbation (Felix and Wagner, 2008). This is the cardinal “What to What?” question discussed by Carpenter et al. (2001) and Lesne (2008). For an in-depth discussion on the biological features that contribute to microbial resilience see the review by Shade et al. (2012). For the sake of evaluating our proposed resilience measure, motility suffices to provide a good working model, with population density as the example function of interest (or as a proxy for community functions that are closely linked to it).

Schematic of model system. A microbial mat has a population of species with differing resilience. In this system, resilience is conferred by movement, with rapid movement (green) having higher resilience than slower movement (red). Both cells respond to a perturbation, in this case light penetration into the mat, by moving out of the zone affected by the perturbation. After the perturbation concludes, the cells recover their initial cell density by moving back into the previously affected region. The faster moving cell has a higher resilience, i.e., the ability to recover population density more quickly, as it is capable of returning more quickly to the previously vacated space.

The simulation was built using Comsol Multiphysics (v. 5.0.1.276) finite element analysis software with the chemical reaction engineering module. The microbial mat model was implemented using the diffusion application mode, with a one-dimensional geometry, similar to models previously described by our group (Renslow et al., 2013a,b,c).

The Concept of the Recovery Length and Proposed Extension

Dai et al. (2013) proposed the half-point recovery length definition of spatial recovery, which increases in proportion to the loss of resilience. The half-point recovery length is defined as L in Equation 1:

where f is the function of interest, x is the spatial distance at the boundary of the perturbation, f(x is the functional value at the corresponding equilibrium (e.g., unperturbed) condition (a diagram of these terms is provided in Dai et al. (2013) supplementary information). For our simulation, we define f(x as the functional value when it has reached 99% of its true equilibrium value, similar to how we calculate temporal recovery, described below. The equilibrium condition may also be considered to be f(x where the function is fully perturbed, especially in the case where recovery length is measured at a boundary or where a full recovery profile (i.e., one that includes both perturbed and unperturbed equilibrium regions) may not be available (see Dai et al., Supplementary Figures 6 and 8 for a more detailed discussion on this special case). In cases where the perturbation boundary is not well defined or unknown, such as during a gradient perturbation, we assign x to be the spatial distance at the start of the perturbation. This does not alter the recovery length–resilience correlation for step-change perturbations in our model, but simply allows for quantifying L in gradient perturbation cases. Also note that, as described in Dai et al., function profiles are normalized by f(x, and all results are shown on normalized scales for clarity and ease of comparison.

A new measure proposed in this work, PASMORE takes into account the shape of the perturbation, whether it follows a sharp step-change or a gradient. PASMORE is defined as a weighted integral, where the perturbation profile is the weighting function:

where f is the function of interest, p is the perturbation profile, and PASMORE is the integral of pf over the defined system space, S. In our simulation, the defined system space is the microbial mat where the perturbation has effect.

Comparison of Recovery Length and Pasmore in Microbial Mats with Non-Interacting Species

In Figure , we demonstrate how a microbial population responds to the application and removal of a perturbation. It is possible to see the relationship between temporal recovery and spatial recovery. Temporal recovery is defined as the time required for the recovery of function, in this case population density, to 99% of its original value. In cases where spatial recovery is correlated to temporal recovery, such as in our simulation (Pearson product-moment correlation coefficient, i.e., Pearson’s r, of 1.0), spatial recovery metrics may be used as a proxy for quantifying resilience.

Recovery over time is associated with recovery through the simulated microbial mat. The top profile shows the bulk population response of a species over time: before (A), during (B), and after (C,D) a step-change perturbation is applied. Bottom figures show the species spatial distribution at each time point. Before the perturbation (A) the species has a homogenous population density. During the perturbation (B), with perturbation profile shown as dashed gray line), the population density is lowest in the perturbed region. After the perturbation has been removed, the species begins to recover (C) until its population density approaches the pre-perturbed profile (D).

It is known that perturbations, whether occurring over time or over spatial domains, will occur frequently across gradients of varying sharpness. For example, using the example of light as the source of perturbation, attenuation follows an exponential decay profile. Even when the perturbation occurs across a gradient, rapidly motile species exhibit high resilience, meaning in this case that they recover population density more quickly than species with slow motility. However, as discussed by Dai et al. (2013), half-point recovery length may fail to accurately quantify resilience under conditions of a blurred or gradient perturbation. We tested step-change, linear, polynomial, and exponential perturbation gradients (data shown only for the exponential case) (Figure ). The ability for spatial recovery to correlate to temporal recovery decreases as the perturbation profiles approaches an exponential function, where all correlation is lost (Pearson’s r of 0.0). Figure demonstrates that PASMORE is able to accurately quantify spatial recovery and maintain the correlation to temporal recovery (Pearson’s r of 1.0) regardless of the perturbation profile. This was true for all perturbation gradients that we tested. Resilience arises due to the motility of the population, thus the metric should maintain its relationship to motility regardless of the perturbation to which it is exposed.

Spatial population density profiles during a perturbation and recovery time – motility relationship for two different perturbation profiles: step-change The faster species (green) recovers more quickly after a perturbation than the slow species (red). Regardless of the shape of the perturbation profile, the recovery time decreases as motility increases, verifying that resilience is higher for the faster species.

Recovery time – recovery length and recovery time – PASMORE relationships for two different perturbation profiles: step-change For a step-change perturbation, recovery length correctly correlates with recovery time (A), however, it is unable to correctly correlate when the perturbation profile is exponential (C). PASMORE correctly correlates to recovery time regardless of the perturbation profile (B,D).

Quantifying Resilience with Limited Perturbation Information

Properly quantifying resilience using spatial metrics requires information about the perturbation’s shape. This is true whether PASMORE or other metrics like half-point recovery length are used. However, in experimental systems, it may not be possible to obtain the entire profile of a perturbation gradient. Here we wanted to quantify the effect of limited information about the gradient. Therefore, we examined PASMORE for use when the exact shape of the gradient imposed by a perturbation was uncertain. As shown in Figure , we compared two cases against the full-profile PASMORE: (1) a two-point dataset to generate a linear approximation between the start and end of the observed perturbation and (2) a three-point dataset to generate two linear approximations between the start and middle and between the middle and end of the observed perturbation. Note that only the perturbation profile was approximated, not the population density, since it is assumed that the function of interest, for which resilience is to be evaluated, is quantifiable. Using sensitive measurement instruments in real systems, it is likely that more than two or three data points of the perturbation profile could be obtained, for example by using microelectrodes (Nguyen et al., 2012) or NMR imaging (Renslow et al., 2010). However, as will be demonstrated, even using a three-point approximation, the calculation of PASMORE quickly approaches the full-profile PASMORE. Figure shows that PASMORE is able to maintain the correlation to temporal recovery with a two-point approximation (Pearson’s r of 0.96), and this correlation improved with a three-point dataset (Pearson’s r of 0.99).

In some cases, a full perturbation profile, i.e., However, it may be possible to make a few measurements and interpolate the perturbation profile, for example at the start and end of the perturbation profile. (A) Three demonstration cases for a measured perturbation profile: Start and end (s.e., linear approximation), start-middle-end (s.m.e., two linear lines approximation), and full profile (full, exponential). (B) Recovery Time – PASMORE correlation to the three measured perturbation profiles. Even in the cases where only limited information about the true perturbation profile is known, PASMORE may still correlate well with recovery time.

Analysis of a Microbial Mat with Interacting Species

We investigated how resilience as measured by PASMORE would change if two species in the mat interacted with each other. Therefore, we modified the model to simultaneously include two motile species that exhibited a difference in their rate of motility, but which was dependent on the consumption of a substrate. The substrate diffused from the top of the mat and consumption followed Monod-type kinetics, based on the diffusion and reaction parameters and setup of a previous model (Renslow et al., 2013b). Such a configuration broadly approximates the structure of a cyanobacterial mat, where the carbon fixation upon which heterotrophs depend is maximal near the mat surface (Lindemann et al., 2013). The case where both species competed for the same substrate was compared to the case where the species consumed two independent substrates noncompetitively. Figure shows the percent change of PASMORE comparing the non-interacting to competitive interaction cases. Compared to the non-interacting case, the interacting species displayed an increase in PASMORE, and thus were found to have lower resilience. Furthermore, the species with the faster motility exhibited a larger loss of resilience compared to the slower species. This is due to the fact that, during perturbation, the species with higher resilience occupy a lower nutrient zone, while the species with lower resilience are closer to the nutrient source and are thus less affected. As such, the model illustrates a trade-off between resource availability and stress avoidance in these motile species; we anticipate that similar trade-offs will exist in many spatially organized communities where local interactions lead to gains in fitness. Finally, this simple simulation demonstrates how PASMORE can help quantify resilience in a case where previous metrics would not have been meaningful.

Schematic of the model system, which has been modified to force competition between two species (one with higher resilience and one with lower resilience) by making motility dependent on consumption of a sole nutrient, which is supplied in the liquid above the microbial mat. During a perturbation, the nutrient concentration far from the perturbation becomes limited. The right plot shows the percent change of PASMORE when a slower species competes for the nutrient with a faster species, compared to a non-interacting case (i.e., when the species rely on differing nutrients). The competitive interaction led an increase in PASMORE, i.e., lower resilience, for both the slow and fast species due to the competition between species. The faster species exhibited a greater decrease in resilience compared to the slower species.

Pasmore in the Broader Context of Spatial Resilience Measures

The capability to monitor spatial patterns of microbial function can help elucidate the resilience of a community in cases where temporal measurements may not have previously been practical. Spatial measures remove the need to monitor a given function over time, and thus a snapshot of the present community stability can be gauged. This type of analysis may have implications beyond the microbial world, as there are many forms of complex networks and communities that display properties of emergent phenomena with spatially-relevant functions. For example, the spatial effects of social dilemmas have been investigated for several decades (Nowak and May, 1992; Hauert, 2006) and the frequency of human cooperation and collaboration or selfishness and exploitation across spatial dimensions impacts the resilience of groups of peoples and their ability to maintain high productivity (Alvard, 2004; Jimenez et al., 2008). Furthermore, spatial patterns of human settlements and population densities may be related to community stability in the face of perturbation events such as loss of water, tillable soil, hunting ground, food, energy and other resources, as well as the related incidents of overpopulation and politico-military conflicts (Tir and Diehl, 1998; Vandam et al., 2013). In this context, densities of refugee settlements across countries and subsequent patterns of repopulation of cities after wars may reveal human community and ethno-regional group stability properties in the face of significant devastation. Other possible application areas for spatial measures of resilience could include aquifer-groundwater recharge rates (Guglielmi and Mudry, 1996; Katic and Grafton, 2011), soil and vegetation health (Seybold et al., 1999; van de Koppel and Rietkerk, 2004; Alongi, 2008; Jiang et al., 2012), permafrost vulnerability (Jorgenson et al., 2010), coral reef habitats (Nystrom and Folke, 2001; Cheal et al., 2013), and fishery management (Carpenter and Brock, 2004; Kerr et al., 2010). Indeed, understanding the resilience of biological and non-biological systems will be critical for determining the impact of climate change and the necessary policy changes to alleviate negative consequences (Bloetscher et al., 2010; Cumming, 2011). Future experiments and observations will need to be done to determine the extent to which spatial measures of resilience, such as PASMORE, may be used to provide insight into stability of complex systems beyond microbial communities. Furthermore, the ability of PASMORE to measure resilience in highly heterogeneous, discontinuous, or physically-restricted systems will need to be tested. It is unclear how barriers that constrain functional recovery, whether social, cultural, financial, political, or physical (e.g., soil and minerals at the microbial scale), will impact the capacity of PASMORE to measure resilience.

Conclusion

‘Perturbation-adjusted spatial metric of resilience’ (PASMORE) is capable of accurately quantifying spatial recovery and maintaining the correlation to temporal recovery regardless of the perturbation profile. In many environmental and ecological systems, the gradient shape of the perturbation may not be quantifiable and even the identity or type of perturbation may not be known. However, even with limited information, PASMORE can be calculated with approximation of the perturbation profile to reveal a close estimate of the system’s resilience. We envision that establishing the concept of PASMORE as a practically useful metric of resilience requires further studies, including investigation of the limits of using spatial recovery in lieu of temporal recovery to quantify resilience in real natural and engineered microbial communities and relating these measurements to predictions of when systems are on the verge of collapse or nearing an irreversible transition across a tipping point. The present work provides initial foundations along this direction.

Author Contributions

All authors listed have made substantial, direct, and intellectual contribution to the work, and approved it for publication.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.