Stochastic disturbance regimes alter patterns of ecosystem variability and recovery.

Altered ecosystem variability is an important ecological response to disturbance yet understanding of how various attributes of disturbance regimes affect ecosystem variability is limited. To improve the framework for understanding the disturbance regime attributes that affect ecosystem variability, we examine how the introduction of stochasticity to disturbance parameters (frequency, severity and extent) alters simulated recovery when compared to deterministic outcomes from a spatially explicit simulation model. We also examine the agreement between results from empirical studies and deterministic and stochastic configurations of the model. We find that stochasticity in disturbance frequency and spatial extent leads to the greatest increase in the variance of simulated dynamics, although stochastic severity also contributes to departures from the deterministic case. The incorporation of stochasticity in disturbance attributes improves agreement between empirical and simulated responses, with 71% of empirical responses correctly classified by stochastic configurations of the model as compared to 47% using the purely deterministic model. By comparison, only 2% of empirical responses were correctly classified by the deterministic model and misclassified by stochastic configurations of the model. These results indicate that stochasticity in the attributes of a disturbance regime alters the patterns and classification of ecosystem variability, suggesting altered recovery dynamics. Incorporating stochastic disturbance processes into models may thus be critical for anticipating the ecological resilience of ecosystems.

Introduction

Understanding and anticipating ecosystem responses to disturbance is a fundamental goal in ecology and conservation that has become even more pressing in the face of global environmental changes that alter disturbance regimes. Although most theoretical and empirical studies focus on mean responses, there is increasing recognition that patterns of variability also contain valuable information about disturbance effects [1, 2]. For example, ecologists have shown that changes in variability can indicate differences in the predictability of responses [3, 4], recovery time [5, 6], and the internal behavior of an ecosystem [7]. Shifts in the variability of an ecological response can also signal changes in ecosystem recovery dynamics and impending ecosystem state transitions [8-10], thereby providing insights into the ecological resilience of an ecosystem, i.e., the magnitude of disturbance that a system can absorb before it changes stable states [11]. Monitoring such ecological resilience indicators may enable actions to avert these outcomes [12, 13]. The ability to predict how natural and anthropogenic disturbances affect variability is therefore of great interest, particularly from a policy perspective, where patterns of variability are increasingly being used to guide decision-making [14-16] and to promote the ecological resilience of ecosystems in a management context [17].

Disturbances can alter patterns of ecosystem variability for several reasons. Within and among ecosystems, spatial and temporal variation in the physical environment and biota modulate disturbance effects and recovery [18, 19]. For example, species’ responses to and recovery from disturbance events depend on their life history, dispersal capacity, and competitive ability [20-23]. Additionally, disturbance regimes are typically stochastic and differ with respect to frequency, severity, and extent, as well as other properties [24]. Variation along any of these axes can subsequently result in different recovery dynamics [25, 26].

Given the inherent complexity of changes in even one of these attributes, it is common to turn to qualitative conceptual models for insight, especially when attempting to generalize across systems of differing properties. Such models can help refine and prioritize questions for further empirical work, as well as inform additional modeling efforts. These conceptual models can also be implemented quantitatively to provide a more explicit understanding of how alternative formulations of key processes influence responses to altered drivers, thus generating testable predictions for future work e.g., [27-29]. Of particular importance to the work on the relationship between disturbance and ecosystem variability is the model developed by Turner et al. [30], which sought to identify the spatial and temporal attributes of a disturbance regime that could lead to qualitatively different landscape dynamics. Turner et al. [30] characterized landscape dynamics by computing the constancy of seral stages (i.e., temporal variability in the amount of cover) under varying combinations of disturbance spatial extent and frequency at a range of scales. The model provides a solid and tractable baseline for discussion of the effects of disturbance regime properties on ecosystem dynamics. However, because the model is deterministic, it does not address how stochasticity in disturbance regime attributes affect ecosystem variability and resultant dynamics. Disturbances are only partially deterministic, and a constrained model may obscure relationships between the properties of disturbance regimes and ecosystem variability. Understanding the effects of stochastic disturbance regimes on ecosystem variability is even more important under a changing environment where altered disturbance regimes will be the norm rather than the exception [24].

Previously, we used the fully deterministic framework of Turner et al. [30] to predict how disturbance might alter the variability of terrestrial and aquatic ecosystems [1]. We compared model predictions with the results of published studies that documented changes in ecosystem variability following disturbance. Although we found good overall agreement between predictions and observations, many outcomes were not as expected. We highlighted the need to incorporate more realistic estimates of species recovery as well as the difficulties associated with properly characterizing chronic disturbance events in the framework as possible improvements, but did not elaborate on the role that stochasticity might play in improving the agreement between the model and empirical studies.

To address this gap, we extended the deterministic model of Turner et al. [30] by modifying the attributes of disturbance regimes to include stochasticity in three principle disturbance parameters: frequency, severity and spatial extent, and investigated which stochastic disturbance parameters might be most influential in generating changes in ecosystem variability and recovery dynamics. We also explored whether stochastic parameterizations improve the agreement between simulated dynamics based on the deterministic and extended stochastic model and results from published empirical studies. Our work demonstrates that accounting for stochasticity in the attributes of a disturbance regime significantly affects the patterns and classification of ecosystem variability, suggesting altered ecosystem dynamics.

Methods

Base model

We duplicated the deterministic model of Turner et al. based on their description of the parameters and processes [30]. The model environment consisted of a set of 100 × 100 gridded cells. Eight vegetation classes representing successional stages were included in the model. Each cell within the grid had a variable that tracked the successional stage and the modelled environment began with all cells at the most mature stage (successional stage 8). The model used a discrete time representation; at each time step a disturbance could be introduced into the environment and any previously disturbed grid cell could recover. Disturbances could occur in all successional stages, were of a fixed size with four equal sides (s) and were imposed at a fixed temporal frequency (f). In this base model all disturbed regions recovered deterministically by passing through the eight successional stages (one per time step) until the mature stage was reached. Eight is therefore also the value we use for ecosystem recovery time (see “S and T parameter space” below). Although Turner et al. [30] created the model to mimic successional dynamics in a plant assemblage, they also encouraged its extension to other ecosystems.

Disturbance spatial extent and frequency have been shown to strongly affect recovery dynamics [7, 20, 30]. In the model, the spatial extent (s2) determines the number of cells affected when a disturbance event occurs. The location of the center of the disturbance extent is determined by selecting coordinates randomly from a uniform distribution. Frequency (f) determines the number of disturbance events per time step during a simulation, and is inversely related to disturbance return interval (1/f), which is the number of time steps between successive disturbance events. Disturbances are evenly spaced, every 1/f time steps throughout the duration of a simulation, and in the case of multiple events per time step (when f > 1), each is executed during that time step (see “Stochastic Frequency” below for further explanation and examples). In the deterministic base model, all disturbances within a simulation use the same spatial-extent and frequency parameters (i.e., the values are fixed). The spatial boundaries of the environment are considered toroidal such that any disturbances that extend beyond the horizontal or vertical limits “wrap” around to the opposite side, eliminating boundary effects. Each model time step proceeded as follows: (1) increase any previously disturbed grid locations by one successional stage (i.e., recover); (2) if a new disturbance event is scheduled, determine center of disturbance randomly; (3) set all newly disturbed grid locations to lowest successional stage. This latter step represents the modelled severity (i.e., impact on the ecosystem) and, in the deterministic base model, a disturbance always resets a grid cell back to the lowest successional stage.

We extended this deterministic model through the addition of stochasticity into all three key disturbance properties: frequency, severity, and spatial extent, as described below. To maintain generality, all parameter values are drawn from a pseudorandom uniform distribution. While some observable phenomena may be more accurately replicated using a more specific distribution (normal, exponential, Weibull, etc.), a uniform random selection provided a reasonable starting point for extending the base model and analyzing its sensitivity.

Stochastic frequency

It is often difficult to predict specifically when a disturbance event will occur because initiating events are typically random (e.g., lightning strike), although over time the average number of events can often be predicted with reasonable accuracy. To ensure comparability with the deterministic model, we maintained the same average number of events over the duration of a simulation, as determined by the frequency parameter f. Before the simulation was run, a timeline was established and disturbance events were scheduled on randomly selected time steps. Consider a simple example with a deterministic frequency of one event every 6 time steps for simulation with 100 time steps. In comparison with the deterministic base model scenario, 16 distinct, evenly spaced events will occur. In our extended stochastic model, the 16 scheduled events will occur but are unlikely to be evenly spaced given the random selection process. Instead they may be clustered or potentially occur within the same time step, as events are drawn with replacement.

Stochastic severity

Like frequency, the severity of a disturbance event was also modelled stochastically. The deterministic model assumes that all disturbance events “reset” a location on the landscape to the lowest successional stage (successional stage 1). While this is reasonable under some conditions, it is not always realistic. For example, fire severity spans a gradient from stand-replacing crown fires to a lightly burned or scorched event [31], and a single disturbance can result in the manifestation of all of these burn severities. Here we implemented stochastic disturbance severity by differentially adjusting the successional stage of a cell. We randomly selected a disturbance severity in each time step and a disturbance was imposed such that the disturbed locations were set back to a successional stage lower than the current stage but not necessarily successional stage 1 (and never lower than successional stage 1). The same level of severity was applied across the entire disturbance extent if more than a single event occurred at the same time. As a result, cells differed in their susceptibility to disturbance. For example, a cell at successional stage 2 would not be “reset” if disturbance severity was ≥ 2. This is analogous to a recovering cell being “skipped” by disturbance. Further, because individual disturbance events could have low or high severity, simulations executed with stochastic severity had lower mean severity than simulations with deterministic severity.

Stochastic spatial extent

Similarly, the spatial extent of a disturbance event is often stochastic. A deterministic disturbance extent is likely for only a limited number of cases (harvesting, controlled burn, etc). We simulated stochastic spatial extent by randomly selecting the width and length (within the 2D environment) of the disturbance from a uniform distribution. The width and length of each event are calculated independently and can range from 1 cell to 2 times the maximum dimension specified in the deterministic model. For example, a parameterization of the deterministic model for disturbance events with spatial extent 20 × 20 cells produces multiple, identically sized disturbances of 400 cells. In our extended model, each disturbance event would have dimensions between 1 cell and 40 cells for each side, producing disturbances with a spatial extent of 420 cells on average. Thus, each event has a unique rectangular shape and size, but the average size of the events over the duration of the simulation are comparable to the corresponding deterministic parameterization.

Fully stochastic

Disturbance regimes are likely to occur on a continuum ranging from largely deterministic disturbances (e.g., spruce budworm infestations, prescribed burning and other forest management practices) to those that are random in frequency, severity and spatial extent (e.g., extreme storm events, floods). Further, because disturbances events can co-occur and interact, we also confronted the model with stochasticity in all parameters by applying the above protocols to all variables simultaneously.

Simulation setup and computations

Model validation

To validate the model and permit visual inspection of recovery dynamics under different model configurations, each simulation was initially run for 100 time steps with deterministic or stochastic parameters (S1 Table). Under the deterministic configuration, disturbance frequency was set to 0.125, resulting in 8 time steps between events, and disturbance extent was 35 × 35 units, while severity was parameterized such that each disturbance reset the successional stage of the disturbed area to 1. We re-parameterized the model for stochastic execution with frequency, spatial extent, and severity of the disturbances simulated separately and in combination as described above. Results were expressed as the mean and variance of the proportion (Vp) of the landscape occupied by the mature successional stage (successional stage 8) over the duration of the simulation. This provides a measure of stability, defined as constancy in structural attributes [32, 33] and forms the basis of our comparisons of ecosystem variability and dynamics in deterministic and stochastic scenarios. The results of model validation are reported in S1 Appendix.

S and T parameter space

To provide general insights about how disturbance regime attributes affect patterns of ecosystem variability and compare across ecosystems and disturbance regimes which can vary greatly in spatial and temporal dynamics, the spatial extent and return interval of disturbance events were expressed relative to the size of the environment and its temporal recovery interval respectively. Specifically, the spatial parameter (S) is disturbance size (s2) / ecosystem size (where ecosystem size is constant at 10,000 (100 x 100 units) and the temporal parameter (T) is disturbance return interval (1/f) / ecosystem recovery time, where recovery time is the number of time steps required for a disturbed cell to achieve the “mature” stage (8 times steps). Using these two key parameters, we defined the state space characterizing the response of ecosystem variability and recovery dynamics to various disturbance regimes.

To compare the effects of stochasticity in this S-T space, we configured and ran the deterministic and stochastic models to generate a wide range of S and T values. Using a factorial approach, we quantified the mean and variance resulting from all combinations of S and T for values of S ranging from 1–100 in steps of 1, and values of T ranging from 0.01–1000 in steps of 10i, where 3 ≥ i ≥ -2 and varies by -0.1, for a total of 50 steps. The small incremental values of S and T enabled us to quantify the mean and variance across the state space at a high resolution (S2 Table). For each unique combination of S and T, we performed a simulation for 100,000 time steps to compute a statistically robust mean and variance. Generally, the model stabilized in < 100 time steps, so the influence of transient dynamics on the calculated metrics (i.e., mean and variance) was minimal.

The results of these simulations were used to identify regions within the S and T parameter space that suggested qualitatively different dynamics based on quantitative patterns in the mean and variance of the proportion of the landscape occupied by the mature successional stage [sensu 30]. These regions represent the range of observed or theoretically predicted dynamics, including steady state dynamics [34, 35], stable dynamics with repeating recovery sequences characterized by low or high response fluctuations [36, 37], and unstable dynamics [38]. The definition of each region is as follows: region A, where the landscape is relatively undisturbed, the mature stage is dominant (> 50% of the landscape) and variance is low (< 5), indicating equilibrium dynamics; region B, where disturbance is more frequent, the mature stage is dominant and variance is low (5–10), indicating stable dynamics with low variance; region D, where disturbance is still more frequent, the mature stage is < 50% of the landscape and variance is low (5–10), indicating stable dynamics with low variance; and regions C and E, where coverage of the mature stage fluctuates substantially producing high (10–20) or very high (> 20) variance, respectively, indicating stable dynamics with high variance; and region F, where disturbances are both frequent and large, producing low variance (< 5) and unstable dynamics, as the mature stage is no longer dominant over the duration of the simulation. Turner et al. [30] suggested that disturbance regimes that place ecosystems in region F might result in state changes.

To evaluate changes in the patterns of variability from stochastic disturbance regimes, we quantified differences in the values of Vp produced by deterministic and stochastic simulations to highlight the parameter space where deviations were largest. We also evaluated how the classification of recovery dynamics changed by comparing the number of pixels falling within the ranges of S and T that defined each qualitatively unique region A-F among the state-space diagrams. For all simulations, the environment was defined as a 100 × 100 grid and the recovery interval was characterized by 8 successional stages. To explore the sensitivity of the model to the number of successional stages and differences in mean severity, we repeated the simulations above using a model with 4 and 16 successional stages.

Comparison of published empirical results and simulated results

We tested whether incorporating stochasticity improved the agreement between empirical results reported in the literature (S2 Table) and results simulated by the deterministic model and our extended model. Studies were located through electronic searches of the Web of Science (Clarivate Analytics) citation index using the keywords: disturbance, heterogeneity, stability, variance and variability, and by examining the references in these studies. Although these search methods may have missed studies, those included here represent a wide range of disturbances, ecosystems and scales, and thus allow for broad inference and testing of our model. We restricted the studies to those that focused on structural (cover, abundance) or compositional responses to disturbance and recorded how the mean and variance changed in response to disturbance to align with the types of responses simulated by the model. We found that the mean did not change much between our models; thus, we focused on variance as the response of greatest interest because it signals where dynamics may shift qualitatively. In sum, we identified 45 suitable empirical cases that focused on either forest (6), grassland (9), freshwater (12), or marine (18) ecosystems (S2 Table).

To classify the response of ecosystem variability and dynamics to the disturbance regime described in each empirical study, we executed deterministic and stochastic simulations using the base and extended models. Models were parameterized using the reported values for disturbance return interval and extent, and five simulations were performed for each study. The first simulation was fully deterministic, with all disturbance parameters having fixed values. In the next three simulations, one disturbance parameter varied stochastically (frequency, severity, or extent) while all other parameters were fixed. For simulations with stochastic extent or frequency, we incorporated stochasticity as described above, by drawing values from a uniform distribution with restrictions to ensure that the average parameter value over all times steps was comparable to the value reported in the study being considered (see “Stochastic Frequency” and “Stochastic Spatial Extent”). For simulations with stochastic severity, severity was randomly selected as described above (see “Stochastic Severity”). In the final simulation, all three disturbance parameters could vary stochastically. We ran each simulation for 100,000 time steps and calculated the mean and variance (Vp) of the proportion of the landscape occupied by the mature successional stage over the duration of the simulation.

The mean and variance computed from each simulation were used to characterize qualitative patterns of ecosystem dynamics for each study and evaluate data-model agreement. As in our previous work [1], we focused on qualitative rather than quantitative patterns to allow comparison across ecosystems and disturbances which can vary greatly in spatial and temporal scales, and to provide insights into conditions under which disturbance produces state shifts in ecosystem dynamics. We used the mean and Vp to classify the dynamics produced by a given simulation as either: equilibrium (A); stable, low variance and dominated by the mature stage (B); stable, low variance and not dominated by the mature stage (D); stable, high variance (C/E combined); or unstable, low variance (F), as described above (see “S and T parameter space”). We likewise classified the dynamics observed in each of the empirical studies based on the reported changes in the mean and variance of the response following disturbance. For example, if we found a minor change in the mean and a major change in the variance, then we classified the study as ‘stable/high variance (C/E).’ To evaluate data-model agreement, we compared the number of studies correctly classified by a deterministic model configuration, a stochastic model configuration, both a deterministic and stochastic model configuration, or incorrectly classified by all models using a chi-square goodness of fit test, where the expected frequency of each category was identical (i.e., 0.25) to indicate an equal chance of being in one of the four groups. Analyses were performed using R statistical software [39].

To facilitate comparison among the different studies, we also computed S and T parameters for each study via the estimation of ecosystem extent and recovery interval using the information supplied in the study. For example, in Coleman et al. [40] variability of marine algal cover increased substantially following release from grazing but algae remained the dominant cover type. The perturbation occurred over half of the experimental area at a rate of once per month, with an estimated recovery of three months. We thus coded the dynamics as “stable, high variability (C/E)” with S = 0.5 and T = 0.3; see S3 Table. When more than one response was examined in a study, we evaluated them separately, apart from responses that had identical values for S and T and showed qualitatively equivalent changes in variance.

Results

S and T parameter space

Under deterministic conditions, the mean coverage of the mature successional stage varied with S and T. When the interval between disturbances was long relative to recovery time (i.e., high T) and/or disturbance spatial extent was small relative to landscape extent (low S), the mature stage dominated (> 50% coverage; Fig 1). Overall, the mature stage dominated across 56% of the state space defined by S and T. Other stages dominated when disturbance spatial extent was large relative to landscape extent (higher S) and relatively frequent (lower T).

Fig 1

Mean coverage of the mature successional stage produced by deterministic and stochastic simulations.

Scenarios include a) deterministic, b) stochasticity in disturbance frequency, c) stochasticity in disturbance spatial extent, d) stochasticity in disturbance severity, and e) stochasticity in all disturbance parameters: frequency, spatial extent and severity. Colors indicate the values of mean coverage of the mature successional stage, which were calculated under varying disturbance regime characteristics summarized by S (expressed as % of the landscape) and T. In the region above the solid white line, the mature successional stage dominates (> 50% of the landscape); in the region below the solid line, the other successional stages dominate.

Implementing stochastic parameters had modest effects on the dominance of the mature stage (Fig 1b–1d). Similar to deterministic simulations, stochastic simulations resulted in the mature stage dominating across 57, 57, 61 and 62% of the state space for stochastic frequency, stochastic extent, stochastic severity, and fully stochastic simulations, respectively. Comparing the mean values derived from deterministic and stochastic simulations, we found that the coverage of the mature stage increased when disturbance extent was large and of intermediate frequency (Fig 2). Smaller and more frequent disturbances also increased the coverage of the mature successional stage, but only under stochastic severity and fully stochastic scenarios (Fig 2).

Fig 2

Quantitative changes in the mean coverage of the mature successional stage from stochastic model configurations compared to deterministic model configurations.

Changes are characterized as the difference in mean coverage of the mature successional stage between stochastic and deterministic simulations, under varying disturbance regime characteristics. Colors represent the magnitude of differences in mean coverage. The spatio-temporal attributes of the disturbance regimes producing the differences are summarized by S (expressed as % of the landscape) and T.

In contrast, stochasticity in disturbance parameters had large effects on variance patterns. Comparing the values of Vp derived from deterministic and stochastic simulations, we found that the variance of simulated dynamics increased most when disturbance regimes had intermediate frequencies and large spatial extents (Fig 3). The smallest changes occurred when disturbance regimes had a relatively small extent, and either low or high frequency. This is likely because when disturbances are small or occur very infrequently, there is less potential for interactions between disturbances, or when they are very frequent, the system is almost always at the lowest successional stage. Stochastic severity decreased the mean severity of the disturbance regime and generally caused more modest changes in variance compared to other stochastic attributes (Fig 3).

Fig 3

Quantitative changes in the variance of the mature successional stage from stochastic model configurations compared to deterministic model configurations.

Changes are characterized as the difference in Vp between stochastic and deterministic simulations, where Vp is the variance in occupancy of the mature successional stage under varying disturbance regime characteristics. Colors represent the magnitude of differences in Vp. The spatio-temporal attributes of the disturbance regimes producing the differences are summarized by S (expressed as % of the landscape) and T.

Comparing a cross-section from the simulations (see Fig 3) at a fixed spatial extent (35 x 35) across the full range of f values represented by disturbance return interval (1/f), we found that variance in the occupancy of the mature successional stage diverged the most for values of disturbance return interval between 0.1 and ~7 (Fig 4). Between these values, stochastic disturbance frequency and spatial extent resulted in higher variance than stochastic disturbance severity, and all these configurations produced higher variance than deterministic parameterizations. Increasing f, which equates to decreasing disturbance return interval, caused a decline in variance because the landscape was increasingly dominated by early successional stages. Decreasing f, which equates to increasing disturbance return interval, caused a rapid decline in variance for all parameter configurations because the landscape was increasingly dominated by the late successional stage. Under these conditions, however, stochastic disturbance extent produced higher variance because of the potential for disturbances to vary in size and to overlap with previous disturbance events, thus causing a greater delay in recovery to the mature successional stage. As we increase or decrease disturbance return time to their extremes, variance decreases as the mature successional stage is either dominant due to lack of disturbance or perpetually in a state of recovery due to multiple overlapping disturbances respectively. However, all models with stochastic elements, show a maximized variance around when disturbance return interval ⁓ 1; i.e., when a single disturbance occurs every time step. At this value, the mature successional stage vacillates between the greatest range of recovery and disturbance conditions, thus producing the highest variance.

Fig 4

Effects of changing disturbance frequency on variance in occupancy of the mature successional stage for different sets of model parametrizations.

For all simulations, disturbance frequency (f) is plotted on the x-axis. In the deterministic parameterization, all parameters are fixed, with disturbance extent = 35 × 35 cells, and disturbance severity is maximal (reset to 1 when disturbed). In the remaining parameterizations, the manipulated parameter is stochastic and all other parameters are fixed with values used in the deterministic parameterization.

One other notable result was the behavior of the fully deterministic version, which showed a sharp decrease in variance as disturbance return interval approached 8 and then a spike in variance as it exceeded 8 (Fig 4). This behavior emerged as an artifact of the deterministic model. Since the number of successional stages and thus ecosystem recovery was set to 8, when 1/f ⁓ 8, variance went to zero because a new disturbance occurred just as the landscape neared full recovery. As a result, the proportion of landscape at successional stage 8 was nearly constant over time. This behavior occurs only in the deterministic model and only when the disturbance return interval equals ecosystem recovery time.

Stochasticity also altered the domains of the S-T parameter space where qualitatively different dynamics emerged (Fig 5). Compared to the regions derived from deterministic simulations (Fig 5a), the regions of the parameter space with stable but high variance (region C) or very high variance (region E) expanded when stochastic parameters were applied (Fig 5b–5e). The size of region C increased by 58, 62, 29, and 55% and the size of region E increased by 130, 111, 28, and 142% for stochastic spatial extent, stochastic frequency, stochastic severity, and fully stochastic simulations, respectively. Thus, the size of regions associated with high variance nearly doubled under most stochastic scenarios. By contrast, all other regions contracted: the size of the region characterized by stable, low variance (B and D combined) was 8–11% smaller, the region characterized by equilibrium (region A) was 1–28% smaller, and the region characterized by unstable, low variance (region F) was 8–19% smaller. For both regions A and F, the smallest change in size occurred when we implemented stochastic severity. The simulations producing the greatest decrease in region size differed by region. For region A, implementing stochastic disturbance extent caused the greatest decrease in region size. For region F, incorporating stochasticity into all the disturbance parameters resulted in the greatest decrease in region size.

Fig 5

Ecosystem dynamics produced by deterministic and stochastic simulations.

Dynamics are characterized by the mean and variance in occupancy of the mature successional stage under varying disturbance regime characteristics. Colors indicate the values of variance in occupancy of the mature successional stage (Vp) over the duration of the simulation and contours indicate the transitions in ecosystem dynamics associated with spatio-temporal changes in disturbance regimes as summarized by S (expressed as % of the landscape) and T.

Overlaying the classified regions shown in Fig 5, we found that the classification of ecosystem dynamics changed the most for stochastic disturbance regimes at intermediate disturbance frequencies (Fig 6). The largest shifts were from states characterized by stable, low variance (regions B and D) to those characterized by stable, high (region C) or very high variance (region E). For the transition between low to high variance states, the most pronounced differences resulted from stochasticity in single parameters. For the transition between low to very high variance states, the most pronounced differences resulted from stochasticity in all disturbance parameters (Fig 6). Stochasticity in all disturbance parameters resulted in relatively larger shifts from the unstable, low variance state (region F) to stable, low (region D), high (region C) or very high variance (region E) states.

Fig 6

Transition maps indicating where shifts in the classification of ecosystem dynamics resulted from stochastic versus deterministic model configurations.

Ecosystem dynamics were characterized with respect to the mean and variance in occupancy of the mature successional stage over the duration of the simulation and relative to the S (expressed as % of the landscape) and T parameter space. Comparisons are made between the fully deterministic version of the model and (a) stochastic frequency, (b) stochastic spatial extent, (c) stochastic severity, and (d) all stochastic disturbance parameters. Gray shades indicate no change in state, and colors represent transitions between states.

When we explored the sensitivity of the model to different numbers of successional stages, we found minimal differences in the sizes of the regions or domains of the parameter space where qualitatively different ecosystem dynamics emerged (Fig 5e and S2 Fig). The main effect of changing the number of successional stages was that mean severity changed concurrently. The simulations had mean severities of 1.5, 3.5, and 7.5 for models with 4, 8, and 16 successional stages, respectively.

We evaluated a total of 45 published empirical responses; 35 of these were previously examined by Fraterrigo and Rusak (2008) and 10 were from new studies (S3 Table). While disturbance extent varied among the studies, over 60% of had S = 1 when normalized with respect to the size of the system. Comparing empirical results with simulated results, we found that 47% of the responses were correctly classified by both the deterministic and stochastic models (Table 1). Of the remaining 24 responses, 46% (11/24) were correctly classified by one of the stochastic configurations and misclassified by the fully deterministic model, indicating an improvement in data-model agreement and an overall revised model agreement of 71% (32/45; S3 Table). Stochastic configurations of the model correctly classified the variability and dynamics of studied ecosystems more often than expected if cases were evenly divided among the four categories in Table 1 (χ2 = 40.7, p < 0.001).

Table 1

Summary of results from an analysis evaluating the correspondence between empirical results reported in the literature and simulated results from stochastic and deterministic models.

Observed response (region)	% of responses correctly classified by stochastic models and misclassified by deterministic model (observed / predicted)	% of responses correctly classified by deterministic model and misclassified by stochastic models (observed / predicted)	% of responses correctly classified by both stochastic and deterministic models (observed / predicted)	% of responses misclassified by both stochastic and deterministic models (observed / predicted)
All responses	24 (11/45)	2 (1/45)	47 (21/45)	27 (12/45)
Stable, high variance (C or E)	22 (8/37)	3 (1/37)	54 (20/37)	22 (8/37)
Stable, low variance (B or D)	40 (2/5)	0 (0/5)	0 (0/5)	60 (3/5)
Equilibrium (A)	33 (1/3)	0 (0/3)	33 (1/3)	33 (1/3)
Unstable, low variance (F)	-	-	-	-

The stochastic models were successful at classifying both low and high variance dynamics (Table 1). The fully stochastic configuration accounted for most of the improvement in classification (11 cases), followed by configurations with stochastic frequency, extent, severity, or (7 cases each). Improved classification of ecosystem dynamics generally occurred for disturbance regimes that had values of S and T that placed them near the boundaries of regions in the state space. This was especially true for cases where the deterministic model predicted that a disturbance regime would cause unstable dynamics (region F) whereas the stochastic model configuration predicted stable dynamics with high variance (region C or E; S3 Table, section 1). One empirical response was correctly classified by the deterministic model and misclassified by stochastic model, and 27% of the responses were misclassified by all model configurations.

Discussion

Disturbance-driven changes in variability are increasingly recognized as ecologically important, yet our conceptual framework for understanding shifts in variability is incomplete [1, 3]. Here, we found that regimes with stochastic disturbance frequency or extent most often led to an increase in ecosystem variance and altered ecosystem dynamics at intermediate frequencies and large spatial extents (Figs 3 and 5). Stochastic severity contributed to modest increases in variance compared to the baseline deterministic case (Fig 3). The incorporation of stochastic disturbance regimes increased the agreement between empirical and simulated results, as evidenced by an increase in the number of studies whose results were correctly classified by stochastic model configurations (Table 1). This increased agreement was partly the result of an overall increase in the size of the S and T parameter space associated with higher variability (Figs 5 and 6). Given that few disturbance regimes are fully deterministic in nature, these outcomes have important implications for understanding how disturbance will affect ecosystem behavior.

In their seminal work relating characteristics of disturbance regimes to landscape dynamics, Turner et al. [30] concluded that disturbance frequency and extent are critically important for determining system stability. In an extension of this conclusion, we found that the largest increases in the variance of the proportion of the landscape occupied by the mature successional stage over all time steps occurred when stochastic disturbance regimes were characterized by intermediate frequency and large spatial extents (Fig 3). At intermediate frequencies of disturbance, the probability that a landscape will be dominated by the mature successional stage or any single successional stage is low. Most of the system is in different phases of recovery, and this results in large fluctuations in structural attributes over time (Fig 4). This pattern is most prominent when disturbance size is large because more of the system is changing asynchronously. Stochasticity amplifies these effects by increasing the likelihood that disturbances will interact, leading to more patches on distinct successional trajectories. At high or low levels of disturbance frequency, the landscape is dominated by early or late successional stages (Fig 1), respectively, and variance in the occupancy of the mature successional stage declines (Fig 4). Under low disturbance frequency (i.e. long disturbance return interval) however, stochasticity in disturbance extent causes a modest increase in system variance through spatial disturbance interactions, which can convert systems from an equilibrium state to a stable, low or high variance state (Fig 6). Using spatially uniform (analytical) models simulating biomass density, Zelnik et al. [23] likewise found that the combination of spatial and temporal dimensions of disturbances can have large effects on ecosystem stability but spatial interactions between disturbed areas contributed to high variance under high rather than low disturbance frequency. This disparity underscores the fact that different measures of stability may respond differently to disturbance regimes.

Our finding that stochastic frequency and spatial extent increased variance in ecosystem dynamics most when compared to deterministic disturbance parameters agree with a recent meta-analysis of the predictability of community response to disturbance which showed that random disturbances generally led to decreased predictability in abundance [3]. For example, disturbances that create open habitat for new colonists have low response predictability, potentially because variation in dispersal, a largely stochastic process, causes differences in how recovering communities assemble [7]. Likewise, variation in distance to propagule sources can lead to alternate successional pathways [41-43]. Stochastic frequency and resultant clustering of disturbance events can also increase the likelihood that an area will be disturbed again before it has fully recovered, which can reduce the abundance of surviving organisms and propagules, substantially decreasing the predictability of succession [44, 45]. Models that represent these processes are needed to distinguish among potential mechanisms driving variance and recovery patterns [23].

Model configurations with stochastic frequency and extent also altered the spatial and temporal domains of the parameter space where qualitatively different ecosystem dynamics emerged (Figs 5 and 6). In particular, the regions in which stable, high or very high variance is expected were larger when these parameters were stochastic. This is because stochastic disturbance frequency and extent increase temporal fluctuations in the abundance of individual successional stages. Despite this, ecosystem dynamics remained stable, with the mature successional stage dominating across more than half of the parameter space (Fig 1). Moreover, stochasticity did not increase the size of the region where unstable, low variance ecosystem dynamics prevail (region F). This result indicates that stochastic disturbance regimes will not always drive systems to collapse in finite time. Martín et al. [46] showed mathematically that spatial heterogeneity can reduce the possibility of catastrophic state change. While we lack the empirical evidence needed to directly test these hypotheses, recent studies investigating changes in species composition and biome boundary shifts (e.g., [24]) suggest that more frequent disturbance events will lead to higher ecosystem variance. In tropical South America, for example, increased fire frequency is expected to cause a substantial reduction in tree cover and expand the forest-savanna transitional area, under less favorable climate conditions [47]. Similarly, an increased frequency of drought events in a sub-tropical lake was coincident with increased variance in net ecosystem production in the lake [48]. Storm events are also predicted to increase their frequency and extent under a changing climate and to the extent this will affect both terrestrial and aquatic ecosystems, the modelled increase in variance suggests an impairment of our ability to predict ecosystem response [4].

In contrast, simulations with stochastic severity resulted in a smaller increase in variance compared to the deterministic model and shifts in ecosystem dynamics from states characterized by high variance to those characterized by low variance (Figs 3 and 5). As implemented, stochastic severity resulted in disturbances that did not always fully reset the community to the earliest successional stage. Consequently, some disturbed areas recovered to the mature stage more rapidly (i.e., in fewer time steps). This suggests that variation in post-disturbance survival could result in greater temporal stability. In fact, biotic legacies are increasingly viewed as important for promoting ecological resilience to disturbance [49]. However, simulations with stochastic severity also yielded lower severity on average over all time steps compared to deterministic simulations (S1 Table). Although we found that this had only modest effects on patterns of ecosystem variance and recovery dynamics (S2 Fig), additional research is needed to disentangle the effects of stochastic severity from those of lower mean severity on ecosystem variance.

We found that the simulation models incorporating stochasticity correctly classified empirical responses more often than deterministic simulations. Correct classification increased from 47% to 71% when some form of stochasticity was introduced into the purely deterministic model. Improved classification of ecosystem dynamics mainly occurred for disturbances with values of S and T that placed them near the boundaries of regions in the state space, particularly for cases where the deterministic model predicted that a disturbance regime might cause unstable dynamics, but the stochastic model configuration predicted stable dynamics with very high variability. This suggests that accounting for stochasticity better constrains the disturbance conditions that are expected to lead to ecosystem instability. Stochastic simulations may improve predictions of how disturbance regimes affect ecosystem dynamics because they account for reduced predictability of structural properties. Supporting this idea, Murphy and Romanuk (2012) found that disturbance commonly resulted in reduced predictability of structural but not diversity-related responses.

That said, our modelled vs. empirical validation is not without its caveats. Many (over 60%) of the studies used in the set of published results have a disturbance size of 1, where the entire landscape has undergone disturbance. While not ideal from a testing perspective, this is the reality of the published literature on disturbance and variability. Although it is conceivable that researchers have a bias toward studying disturbances of large spatial extent, this may also truly reflect the distribution of disturbance size found in nature. Certainly, aquatic environments are typically disturbed at the whole ecosystem scale due to their spatial homogeneity. Further, given the necessary requirements for inclusion in this comparison, the total number of applicable studies was not large, and the resultant power of our comparison is diminished. Nonetheless, an overall classification improvement of nearly 25% is noteworthy and highlights the need for an increased consideration of variability as a response in studies of ecosystem dynamics.

Ultimately, any model, no matter how complex, offers a simplified view of ecosystems and our refinements to the deterministic model of Turner et al. [30] are no exception. We offer a simple model of ecosystem dynamics based on a step-wise theory of succession that has since been replaced by more sophisticated and nuanced theories [50-52] and no way for competition, dispersal, or legacy effects to operate in the context of succession dynamics [49, 53, 54]. However simple models that can be largely substantiated by empirical results [55] and possess an ability to predict outcomes from observations not used in model design and parameterization [56] can be very powerful tools. Interaction strengths can often be weaker than expected [57, 58], competitive interactions are often not dominant in many ecosystems [59, 60], and the dynamics of species replacement, whether successional or not, certainly operate in ecosystems–both terrestrial [50] and aquatic [61]. In this context, the simplicity of our model allows us to demonstrate that basic attributes of disturbance can generate potentially diverse outcomes in ecosystem variability. Those studies that do not agree with predicted changes in system variability offer the potential for insights into additional mechanisms that may prove to be important in incorporating into an improved framework for forecasting the effects of disturbance on ecosystem dynamics. One promising direction for future studies could include alteration of the degree of stochasticity incorporated into disturbance attributes (i.e., the range over which frequency, extent, and severity can vary). For example, using a spatially explicit simulation model, Fraterrigo et al. [27] evaluated how different levels of environmental stochasticity influenced population dynamics on landscapes that varied in habitat configuration and found that higher levels of stochasticity increased variability in population size. Incorporating a greater range of stochasticity may therefore substantially shift regions in S and T space and capture a larger percentage of studies, thus increasing our understanding of the role that these individual attributes play in altering the ecological resilience of ecosystems.

Understanding the effects of stochastic disturbance regimes on ecosystem variability should provide important practical information to natural resource managers. Managers have appreciated the vital role of disturbances in shaping ecosystem dynamics for some time, relying on natural variability concepts to guide contemporary ecosystem management [14, 15]. Our results provide new insights into the disturbance conditions under which altered dynamics can be expected. Specifically, our results suggest that stochastic disturbance regimes will often yield stable, but high variance dynamics. However, because human activities and climate change are altering environmental conditions, disturbance effects are changing [24]. For example, fire and windthrow have triggered rapid changes in tree community composition in southern boreal and northern hardwood forests across central North America where climate change is shifting the prairie-forest ecotone [62]. Consequently, more sophisticated approaches can build upon these outcomes to improve quantitative predictions about how the characteristics of disturbance regimes will affect ecosystem response variability under anthropogenic environmental change.

Parameters and results for the deterministic and stochastic simulations used to validate the model and explore ecosystem variability and recovery dynamics.

(DOCX)

Click here for additional data file.

Spatial (S) and temporal (T) scales of disturbance used to create state-space diagrams illustrating potential disturbance dynamics.

(DOCX)

Click here for additional data file.

Empirical studies evaluated for data-model agreement.

(DOCX)

Click here for additional data file.

Occupancy by successional stage for simulations with deterministic and stochastic disturbance parameters.

(DOCX)

Click here for additional data file.

State-space diagrams describing potential ecosystem dynamics for fully stochastic (i.e., stochastic spatial extent, frequency, and severity) models with 4, 8, or 16 successional stages with respect to temporal and spatial disturbance regime attributes.

(DOCX)

Click here for additional data file.

Model validation results.

(DOCX)

Click here for additional data file.

References for studies cited in S2 Table.

(DOCX)

Click here for additional data file.

(ZIP)

Click here for additional data file.

22 Oct 2019

PONE-D-19-25590

Stochastic disturbance regimes improve prediction of ecosystem variability and dynamics

PLOS ONE

Dear Dr. Fraterrigo,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

You manuscript has been evaluated by three expert reviewers and myself. Both Reviewer #1 and #2 only raised minor points about your study, whereas Reviewer #3 has more important concerns about your study. I partially agree with Reviewer #3. One the one hand, I understand that your aim is to introduce stochasticity on some parameters that Turner's original work considers deterministic and in that sense your study provides a sufficiently novel contribution to deserve publication in PLOS One (it is not revisit of Turner's work). One the other hand, I think that some parts of the manuscript, such as the heading of the section in which you revisit Turner's model, could be modified to make this distinction as clear throughout the manuscript as it is in the Introduction. In that sense, Reviewer #3 provides very detailed suggestions. Because of this, I am inviting you to resubmit a Minor Revision of the manuscript, but I am planning to send the manuscript back to Reviewer #3 to ensure that the main four points raised in his/her report have been fully addressed.

We would appreciate receiving your revised manuscript by Dec 06 2019 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Ricardo Martinez-Garcia

Academic Editor

PLOS ONE

Journal Requirements:

1. When submitting your revision, we need you to address these additional requirements.

Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors present an investigation of an extension of the nonequilibrium landscape model by Turner et al. [20]. The model assumes that disturbances on the landscape have stochastic features (occurrence, size, etc...) mimicking realistic conditions. Comparing the simulation results with experimental data, the authors find a good agreement when accounting for stochasticity, in contrast with a deterministic (effective equivalent) model, which fails to predict landscape properties in many cases.

The manuscript is well written and the results are sound and well explained. However, I have minor concerns about the current version as I detail below.

1. The authors strongly emphasize that accounting for stochasticity improves prediction of the variability classification. I understand the point, but it is natural that making the model more complex would increase your predictability. This would occur for any pertinent generalization of the Turner's model. The main message is that stochastic regimes significantly impact on ecosystem variability classification. Of course, when considering it, you increase your predictability. Furthermore, the claim that 'stochasticity increases prediction' gives the idea that a stochastic resonance phenomenon is in place, i.e. stochasticity is enhancing a certain feature of the system that you are trying to measure. So, I would suggest the authors to review this terminology used in the manuscript, being careful in order to allow the results to reach with more clarity a broader audience.

2. Simple models indeed can help understand the behavior of ecosystems. But, since the proposed one assumes non-interacting units, two points come to my mind. a) I don't see how external perturbations can affect the ecosystem dynamics if the entities are non-interacting (e.g. Ising model). In my point of view, despite the fact that the variability regimes are being affected, the dynamics itself is the same (system response time properties, for instance). So, I suggest the authors to review the statement that stochastic regimes affect ecosystem dynamics (this is related to comment 1), verifying if it is really applicable in the present context. b) Analytical considerations to support the numerical results seem to be possible. In this sense, I missed basic (mean-field) relations (I would guess that the contour lines in Fig.2 are given by T = const. x S, which with y-log scale cannot be automatically understood). Besides this, it seems that a complete theory for the present model, or for other versions, could be developed, or might be already present in the literature. So, including a discussion or references about this aspect might be helpful to others understand the development on the theoretical side.

Reviewer #2: "Stochastic disturbance regimes improve prediction of ecosystem variability and dynamics" is a very interesting and well-written work.

The study focuses on an important contemporary problem and highlights the importance of considering disturbance fluctuations for a better understanding of the response of actual ecological systems.

Contents are very clear, well organized and exposed. The analysis is done carefully.

Comments:

1. Although already said that extensive analysis is left for future studies, I think it is important a longer discussion about the robustness of the results to other parameter values and degrees of stochasticity.

2. The considered dynamics probably experience transient and stationary regimes. For most of the analysis, 10000 time steps are considered, probably reaching the stationary state. Results would then be averaged over transient and stationary regimes, which have different behaviors. I think it is important to discuss this point and clarify its validity.

3. Concerning the previous point, in section (i) 100 time steps are considered, showing the transient dynamics. I propose to discuss a case with longer times.

Minor comments:

1. In line 122, "spatial extent" has not been defined previously.

2. In line 134, it is not clear how frequency is implemented.

3. In line 143, the concept of "severity" remains slightly unclear.

4. In line 176, I think it should be <=2.

5. In line 184, "width and height" is sometimes changed to "width and length". One of them should be chosen. Maybe "width and length" would be more appropriate for 2D systems.

6. Defined regions are slightly confusing. In particular, "stable" refers to both >50% and <50% of mature stage dominance. For example, In Fig. 3 regimes B and D are both described by "stable, low variance".

7. In line 276, "Comparison WITH OF published..."

8. In line 356, "spatial extent" should be "mature successional stage"?

9. In line 379, as below, the % of the expanded region can be indicated.

10. In line 417, clarify that it is the overlaying of Fig. 3?

11. In Fig. 1, stages 2,3,4,5,6 are not well distinguished with the selected dashed lines.

12. In Fig.2, yellow indicates 100% and black 0% instead of >50% and <50%, respectively.

13. In Fig.2, the deterministic case is not shown for comparison.

14. In Fig. 2, 3, S1 and S2, colorbox values are not defined in the captions.

15. In Fig. 2, 3, S1 and S2, colorbox values are difficult to read because of the small font size.

16. The analysis of Fig 4. could be done for other parameters. Maybe it can be discussed what happens for the rest of the parameters.

Reviewer #3: The manuscript by Fraterrigo et al. looks at how various forms of stochasticity of disturbances affects ecosystem dynamics using a simple theoretical methodology, and use previously published work to test how various ecosystems conform to their model. I believe that the goal of the authors is a very important one, namely better understanding of how stochasticity impacts ecosystem dynamics and predictions related to them, and trying to test this on data from various ecosystems. However, I do not think that the approach the authors took to do this is the most appropriate, and I believe it has several important flaws.

These main issues are:

1) There is some discrepancy in what do the authors aim to achieve in this paper? Is it a revisit of the classic Turner paper, or a test of the role of variability in various disturbance parameters? The way the analysis is done and presented it seems that this is a revisit of the Turner paper (e.g. Fig. 1 and its specific parameters that are the same as in the Turner paper), while the discussion and general claim of the paper is that this is an attempt to understand the role of variability. I think the authors need to choose one of these. Also, assuming that they want the second option, then a better discussion or analysis on the role of the different stochasticity sources, and in particular what happens to their combination, is quite necessary. See also points below: 13, 17 ,18, 20, 21, 25, 26, 32, 37.

2) Much of the assessment/comparison is based on the arbitrary distinction to 6 regions (A-F). Saying that something is significant because it changes the classification between these 6 regions, without explaining how they are qualitatively different or why this classification is objective, is flawed. I believe that this is an issue of the mixture between "two papers" I mention in my previous point. See also points below: 16, 24, 25, 28.

3) The comparison to data is problematic, and seems very skewed/arbitrary to authors choices of how to classify a given system. The way the authors chose to classify each system certainly needs to be better explained, but it also seems that it is quite arbitrary, and therefore why we should believe that this classification and its repercussions makes sense needs to be explained. Regardless, given this issue, and what seems like a problematic data set, it seems that the main claim of the paper, that this addition of stochasticity improves predictions, needs to be toned down. See also points below: 14, 33.

I would like to clarify some things regarding these 3 points. First, I do appreciate the choice of using a simple model, and focusing on qualitative information to test this on data. While I have some issues with the choice of model taken and its applicability, this is not in itself my main concern, but rather the discrepancy between the proposed goal of understanding stochasticity and the actual work done, which in practice mainly focuses on this models's specific results. I also appreciate that finding and collecting the data on various systems is not a simple task, and that most likely there is not enough data out there currently to preform a more detailed analysis. However, if this is the case, then this issue should be better addressed (by explicitly stating the limited power of the test preformed, toning down the strong claim (e.g. title) the authors make, and perhaps mentioning in the discussion that better data is needed).

I list below my other comments and issues, mostly organized by line numbers. While they are not major on their own, many of them make up the main issues I raised above.

4) I think in the introduction some more context is needed on previous work. In particular other works looked at similar issues of disturbance properties (e.g. Miller 2011 PNAS, Miller 2012 Eco Res), many studies looked at other system properties that are ignored here in relation to disturbances (e.g. Leibold 2004 Ecol. Lett. on dispersal, Zelnik 2019 Frontiers on local dynamics), and other models have been used to look at similar issues (e.g. Moloney & Levin 1996 Ecology, McCabe&Dietz 2019 Frontiers).

5) lines 59-71 have some strange phrasings of the topic - line 59: "the response of ecosystem variability", line 61: "disturbance-driven changes in the variability of a response can indicate differences in", line 71: " altered ecosystem variability". I guess you mean to talk about about "changes in variability", how they might matter or inform us, etc. But this is currently unclear.

6) lines 73-75: It is not true that there are 3 relevant properties of disturbances (frequency, severity & extent). For instance, species composition (i.e. how different species are sensitive to a given disturbance) might matter significantly, see for instance Arnoldi et al. 2018 JTB.

7) lines 83-89 - The purpose of this text is not clear, and seems like it is mostly unnecessary.

8) line 152 - The authors themselves defined disturbance frequency in a previous publication "It is widely recognized that disturbance frequency, defined as the mean number of events per time period...", in a way that is not consistent with talking about "stochastic frequency" versus "constant frequency", since both would be called frequency in any case. It would be better to talk about return interval for this definition, if not elsewhere in the manuscript.

9) lines 153-162 this seems ill defined. Are you really always having 16 events? In this case, you should better describe your randomization process (is this a point process?)? Or do you just have 16 disturbances on average? (which would make more sense)

10) lines 196-200 seem counter productive - there is stochasticity in all these processes (even those controlled by humans), it is more the question of how much variation exists, and therefore, how much it might matter that they are stochastic. Trying to untangle the roles of these different forms of stochasticity, and their combination, is certainly a good idea, but the logic of doing so should be stated clearly.

11) lines 207-208 and lines 249-250 give very different runtime numbers for simulations - why?

12) line 286 - what does it mean that "We focus on variance"? If it is simply that the mean doesn't change much between your 4 models, why not say that explicitly?

13) lines 292-304 - is there anything different between what you did for the model runs here and for the T-S parameter space? I guess not, and if so, this whole explanation is cumbersome and unnecessary, since as I understand it, you basically estimated where the different systems (from publications) are on the T-S parameter space, given values of S and T for each system.

14) line 311-313 - How is this classification done? This is a critical point, to which you basically give no information, which is particularly noteworthy considering you give plenty of (often too much) information about the specific model choices and technicalities.

15) lines 319-320 are strange. You say you measure S and T to "facilitate interpretation"? But, without estimating S and T, no comparison can be done - is that not so?

16) lines 331-334, at least on their own, seem unsubstantiated. Why are differences in variance levels indicating "distinct ecosystem dynamics". No qualitative shift, if one exists, is mentioned. Also, why were this specific parameter set chosen, besides the fact that it was used in the original Turner paper? How general is the result that averages do not change much, but variances do?

17) lines 336-341 - Overall I would say that Fig.1 chooses parameter that are not very representative. Given this specific frequency choice, disturbances cannot interact, but this is not an obvious choice. For any frequency that is higher interactions are likely to occur, and this should be mentioned.

18) lines 346 and 352, "indicating high ecosystem stability" and "indicating reduced ecosystem stability" - what is meant by stability here? If it is the lack of variability, then this is correct, except there's no need for "indicating" - it is a fact. If it is something else, then that needs to be better explained.

19) lines 347-364 are mostly stating things that are quite evident. I think they can be shortened or clarified.

20) lines 368-369 - but how much is the variance compared to the "extra" variance from all 3 separately?

21) lines 371-377 are again almost trivial. A multiplication of T*S would tell me the baseline values I would expect for the mature stage. What would be interesting to know is how it drifted from this baseline.

22) lines 378-379 a bit vague, and not something that is apparent from the figure itself. How much actually changed. And can it be accounted for by "adding" the 3 separate components?

23) lines 386-389 appear to be wrong. Did you mean that the areas below and above 50% are separated by the white line?

24) lines 391-396 Amplified by how much? Perhaps a map of this would be more illuminating then some of these other measures that use arbitrary classifications.

25) lines 397-399 comparing the fully stochastic scenario, with extent of frequency scenarios, region C grows less for the fully stochastic, but region E grows more. Is this because they are related to each other, or is something else happening?

26) lines 402-407 tell us things we could see from before, which is fine. But what about an explanation of what is happening? In particular, is the severity low variance related to the specific choice where it can only be lower for the stochastic scenario? What if you compared a more appropriate scenario?

27) line 413 - "major transitions in ecosystem dynamics", again the claim that this represents major transitions seems unsubstantiated.

28) lines 418-420 - It seems that you are saying that near the "boundaries" between your regions (A-F), there are more changes between regions due to addition of different sources of stochasticity. This is to a large extent a form of tautology, since changes near an edge of a region are more likely to lead to a different region due to geometric considerations alone. Or did you mean something else?

29) lines 427-429 - this should be given in context. It is an artifact of the parameter space, e.g. stochasticity means less when disturbances occur very infrequently (no interactions between disturbances) or when they are very frequent (system is almost always at the lowest stage).

30) lines 439-441 - how is this different than before? Throughout the paper you test 4 different stochastic models, and it seems that here as well. So what is different about this? Also, it is not clear what parameter you used here to get Fig. 5.

31) lines 442-445, I guess you meant "f" instead of "T"?

32) lines 439-451 - key insights are missing here. The fact that for T>8 (small f) you get a larger difference between the deterministic and stochastic models is not very surprising, considering that for the deterministic model you cannot get interactions of disturbances, but construction. when T<8 (high f) this changes, and things become more comparable. Interestingly, this is exactly when the extent-only model has larger variance than others, which is consistent with the fact that it doesn't have interactions for T>8 (i.e., the extent-only model "should" get larger values overall, but for T>8 this is not seen). Also, why is the peak around T=1? Is that a general result?

33) lines 467-476 - over 60% of papers were with disturbance size=1. This seems quite problematic for the analysis, and at the very least should be mentioned.

34) lines 502-505 - are you trying to say that stochasticity matters, and that variability and stochasticity are related? That alone seems redundant. If you are saying something else, it should be better clarified.

35) lines 506-519 seem to just repeat results and conclusions from the original Turner paper. And lines 519-527 are quite strange and seem out of place. In particular, talking about dispersal while it is not mentioned at all in the whole paper, would only make sense if this was to discuss how adding dispersal would affect the results - and this is not what the authors do here.

36) lines 543-550 seem to try and make several points, which seem either problematic or that I could not follow properly. In the first sentence you say "in accordance" - of what to what? Are you just saying that the addition of variance due to stochasticity in frequency and extent lead to larger regions ("stable, high or very high variance")? Why does this "suggest" that there is increased fluctuations? And regardless, why not check that explicitly? Also, why did "ecosystem dynamics remained stable..."? Stable in what sense? And how does this indicate that stochasticity does not need to drive systems to alternative states? Certainly, that statement as is (I would rephrase it as "stochasticity will not always drive systems to collapse in finite time") is true, but it is also widely known, and rather, any divergence from it would be surprising. Finally, one point that I suppose you tried to make, but I am not at all sure, is that

37) lines 562-572 point out to a major modeling problem, which is that the comparison of deterministic vs. stochastic severity of disturbances was not done properly. The last sentence calls on more research in this direction, but I think that the minimal thing you needed to do was test this directly, even if in a limited fashion (e.g. compare deterministic model with severity 8 to a stochastic one where severity is randomly chosen between 15 an 1)

38) lines 586-606 - I think it is here that you need to mention the effect of dispersal, that your model does not account for, similarly to how you do this with issues of internal dynamics.

39) Figure S1 does not use colors, and the color axis is not the same for all panels, making visual comparisons difficult.

40) "Turner et al. [20]" is used many many times throughout the manuscript. This can often be shortened, e.g. "deterministic model"

41) line 71 - "Disturbance" should be "A disturbance" or "Disturbances", no?

42) line 253 - typo: "produced by using"

43) line 276 - typo: "Comparison with of published"

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: Yes: Yuval Zelnik

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

22 Jan 2020

See attached file "Response to Reviewers" for responses to specific comments.

Submitted filename: Response to Reviewers_final.docx

Click here for additional data file.

19 Feb 2020

Stochastic disturbance regimes alter patterns of ecosystem variability and recovery

PONE-D-19-25590R1

Dear Dr. Fraterrigo,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Ricardo Martinez-Garcia

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #3: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #3: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #3: N/A

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #3: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #3: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #3: I have read the revised manuscript by Fraterrigo et al. and I find it to be a significant improvement over the original submission. The paper reads more easily, the introduction and discussion are more appropriate to the results, and overall the paper's contribution is more clear. I am therefore happy to recommend it for publication.

I have a few minor suggestions, which the authors may find useful, that I write below.

Line 373, "Stochastic severity generally caused more modest changes in variance compared to other stochastic attributes (Fig 3).", I think you should say that this is result of the decrease of average severity (as you basically later show in lines 465-470)

Comparing Fig.2 and Fig.3, I cannot visually see that there are big changes in variance patterns (Fig.3), and not so much in average cover (Fig.2). I don't think this is necessarily a major point, but I would recommend either clarifying (visually or otherwise) why/how the variance changes are more substantial than those for average cover, or putting less emphasis on this point.

If possible, it would be great if you could clarify further how you classify into the 4 categories. Looking more into this, I realize that in practice you see 3 categories in the studies you found. I think it would be very useful if you therefore clarify the "practical procedure" for each such category (e.g. "If we see that there has been a major change in cover and also X, than we note this as category Y").

In the legend of Fig.6, you write "(e)" instead of "(d)"

A minor note, you spelled my name in the Acknowledgments as Zelnick instead of Zelnik.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #3: Yes: Yuval Zelnik

24 Feb 2020

PONE-D-19-25590R1

Stochastic disturbance regimes alter patterns of ecosystem variability and recovery

Dear Dr. Fraterrigo:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Ricardo Martinez-Garcia

Academic Editor

PLOS ONE