## Abstract

Disturbance is key to maintaining species diversity in plant communities. Although the effects of disturbance frequency and extent on species diversity have been studied, we do not yet have a mechanistic understanding of how these aspects of disturbance interact with spatial structure of disturbance to influence species diversity. Here we derive a novel pair approximation model to explore competitive outcomes in a two-species system subject to spatially correlated disturbance. Generally, spatial correlation in disturbance favoured long-range dispersers, while distance-limited dispersers were greatly suppressed. Interestingly, high levels of spatial aggregation of disturbance promoted long-term species coexistence that is not possible in the absence of disturbance, but only when the local disperser was intrinsically competitively superior. However, spatial correlation in disturbance led to different competitive outcomes, depending on the disturbed area. Concerning ecological conservation and management, we theoretically demonstrate that introducing a spatially correlated disturbance to the system or altering an existing disturbance regime can be a useful strategy either to control species invasion or to promote species coexistence. Disturbance pattern analysis may therefore provide new insights into biodiversity conservation.

## 1. Introduction

Disturbance has long been recognized as a critical driver of species coexistence and diversity in community ecology [1–6]. Disturbance events can be biotic (e.g. grazing, pest outbreaks) or abiotic (e.g. flooding, fires, drought), and also anthropogenic, such as selective harvesting [7]. The ecological effects of disturbance are complicated and multifaceted, as disturbance events may vary in multiple quantifiable aspects, such as timing, intensity, duration, periodicity and extent [8,9]. To date, most theoretical (e.g. mean-field approximation) and empirical work has focused on exploring how these individual characteristics of disturbance affect species persistence and therefore biodiversity [8–18], while ignoring the spatial structure of disturbance. Yet many natural and anthropogenic disturbances do not occur uniformly across the landscape (e.g. fire, drought, flooding, pest outbreaks, forest logging).

More recently, theoretical modelling has recognized that the spatial characteristics of disturbance, not just disturbance frequency and intensity, also regulate population dynamics [19–25]. In particular, several single-species models have indicated that spatial aggregation in disturbance strongly impedes the growth of distance-limited dispersers, while long-range dispersers are much less influenced [20,24]. This may significantly alter species coexistence and community diversity even if disturbance frequency is not changed [25–27].

Despite intensive studies, whether and how different aspects of the disturbance regime interact in determining species coexistence is much less understood. Miller *et al.* [8,9] theoretically documented the importance of the interaction between disturbance extent and frequency for competitive outcomes, but their work ignored the spatial correlation of disturbance. Banitz *et al.* [25] and dos Santos *et al.* [27], on the other hand, using spatially explicit individual-based models (IBMs), highlighted the individual impact of spatial correlation on biodiversity, independent of disturbance extent, finding that species diversity is very sensitive to spatial correlation when the trade-off in life-history traits includes differences in dispersal distances. However, a mechanistic understanding of the interactive effects of these common aspects of disturbance (including extent and periodicity) and its spatial aggregation on long-term species coexistence is not well developed, especially in regard to the relative importance of these aspects in influencing competitive outcomes among different dispersers. Furthermore, there have been no formal theoretical models to assess under what conditions spatially aggregated disturbances allow species coexistence.

In this study, we model species dynamics in the face of spatially correlated periodic disturbance by using a pair approximation (PA) approach, which has been shown to be useful in characterizing spatial neighbouring correlation between adjacent sites in lattice-structured models [28–35]. Because we expect the spatial aggregation of disturbance to interact with dispersal strategy, we consider two extremes of plant species dispersal: distance-limited, vegetative, clonal spread (local dispersal [36,37]) and long-range, wind-based seed dispersal (global dispersal [38–40]), with local dispersal being restricted to neighbouring sites and global dispersal being uniformly random across the entire landscape. Two mechanisms are included to promote species coexistence: (i) a trade-off between colonization ability (e.g. dispersal distance) and mortality rate [41], and (ii) the intra- and interspecific competition mechanism such that stable species coexistence is possible when heterospecific competition is weaker than intraspecific competition [35,42]. With the model, we systematically investigate how the aggregation level of disturbance interacts with its spatial extent (i.e. area) and periodicity (i.e. frequency) in modifying species coexistence. This sheds more light on the conditions under which spatial correlation in disturbance promotes long-term species coexistence, and allows assessing their relative importance to competitive outcomes.

## 2. Material and methods

### (a) Pair approximation model

We explore species coexistence in an infinite lattice-structured landscape where each habitat site can be either empty (denoted by 0) or occupied by an individual belonging to species 1 or species 2. The two competitors have different dispersal modes, local for species 1 and global for species 2. Each site has *z* = 4 nearest neighbours (von Neumann neighbourhood). Based on Ying *et al.* [35], the dynamics of the local and global disperser can be described by
2.1*a*

2.1*b*with all parameters defined as in table 1. In equation (2.1*a*), the first term on the right represents the population growth of the local disperser, where is the average probability of vacant sites surrounding the individuals of species 1 in the entire landscape, as three possible states—empty (0) or occupied by species 1 or 2—exist for a randomly chosen neighbouring site of the target individual 1. The second term denotes the sum of the intrinsic mortality (*m _{i}*) and the increased mortality due to both intra- (

*γ*

_{11}

*q*

_{1/1}) and interspecific competition (

*γ*

_{12}

*q*

_{2/1}), in which species competition is determined by both contact frequency (i.e. local density

*q*

_{1/1}or

*q*

_{2/1}) and species sensitivity to neighbours (intra-

*γ*

_{11}or inter-

*γ*

_{12}). In equation (2.1

*b*), the first term on the right side represents the growth of species 2 via seed random dispersal, where (1 −

*ρ*

_{1}−

*ρ*

_{2}) =

*ρ*

_{0}is the fraction of unoccupied sites in the entire landscape. Note that we only consider seeds landing on empty sites and ignore those landing on already occupied sites [33,34]. Similar to species 1, the death of species 2 also includes both intrinsic mortality and the enhanced mortality induced by both intra- and interspecific competition.

In order to construct a closed mathematical system for two competitors with different dispersal strategies in a lattice-structured landscape, we further derive the dynamics of local density *q*_{1/1}, *q*_{2/2} and *q*_{1/2} (note *q*_{2/1} = *ρ*_{2}*q*_{1/2}/*ρ*_{1}; see electronic supplementary material, ‘Supplementary methods’). Therefore, equations (2.1*a,b*) and electronic supplementary material, equations S9–S11 with five variables (*ρ*_{1}, *ρ*_{2}, *q*_{1/1}, *q*_{2/2} and *q*_{1/2}) fully characterize the system's dynamics in a constant environment without involving spatial disturbance. The system stability analysis is given in electronic supplementary material, ‘Supplementary model analysis’, which rigorously demonstrates analytical conditions for model outcomes under some simplifying restrictions.

### (b) Spatially correlated periodic disturbance

We assume that each disturbance event occurs instantaneously (i.e. a pulse disturbance, such as hurricanes, flooding, fire, etc.), killing all plant individuals within the disturbed area in a single time step of the simulation (illustrated in figure 1). However, such a disturbance does not influence habitat quality; disturbed sites thus remain suitable for subsequent species recolonization [44].

In our model, the disturbance regime is controlled by three parameters: extent (*ρ _{d}*, the fraction of disturbed area in the entire landscape;

*d*, disturbed site), spatial correlation (i.e. aggregation of disturbance,

*q*

_{d}_{/d}) and periodicity (a fixed interval between two disturbances,

*T*) (table 1). Similar to local density mentioned above, the parameter

*q*

_{d}_{/d}, which is defined as the conditional probability that a randomly chosen neighbouring site of a disturbed site is also disturbed, describes the spatial clumping degree of disturbed sites, with

*q*

_{d}_{/d}=

*ρ*/

_{dd}*ρ*. The parameter

_{d}*ρ*is the pair density of

_{dd}*d–d*sites in the whole landscape (i.e. the probability when randomly choosing a pair of neighbouring sites that both of them are disturbed).

We apply disturbance in a spatially correlated way as a stochastic dynamic process, meaning that the spatial locations for each disturbance event are random, but with a given degree of spatial aggregation. According to the orthogonal neighbouring correlation algorithm (see more details in papers by Hiebeler [31,43]), the range of *q _{d}*

_{/d}should obey 2.2when we classify habitat sites into two states: disturbed (

*d*) and undisturbed (

*u*). Thus, the allowable aggregation of disturbed sites (

*q*) is determined by the disturbed area (

_{d/d}*ρ*). The disturbed sites are randomly structured when

_{d}*ρ*=

_{d}*q*, while the cases of

_{d/d}*ρ*<

_{d}*q*and

_{d/d}*ρ*>

_{d}*q*denote spatially aggregated and over-dispersed (scattered) disturbance, respectively.

_{d/d}In our model, disturbance occurrence is periodic and deterministic. Let *t* = *nT* () be a time period in which a pulse disturbance occurs, and *T* the time interval between two pulse disturbance events. Just after a pulse disturbance at *t*^{+}, individuals within the disturbed area are removed immediately regardless of their species identity, and a new species pattern is generated (figure 1). A pulse disturbance can reduce population density with (*i* = 1 or 2; *t* = *nT* with ), where *ρ _{i}*(

*t*) and

*ρ*(

_{i}*t*) represent the population density, respectively, just before and after the pulse, and

^{+}*ρ*is the fraction of disturbed area in the whole landscape. When the periodic disturbance is incorporated into the PA model presented in equation (2.1), the system becomes non-autonomous as follows:

_{d}

similar to the SIR (susceptible–infected–recovered) system with pulse vaccination [45,46].

On the other hand, spatial disturbance can change species pattern, so we need to consider how to characterize the dynamics of local population density (*q _{i/j}*) in the face of periodic spatial disturbance. Just after a disturbance pulse occurring at

*t*=

*nT*(), we have (

*i, j*= 1 or 2). If the states of the pair

*i–j*sites are unaltered at

*t*

^{+}, then these sites are both undisturbed, which has a probability

*ρ*(i.e. the probability when randomly choosing a pair of neighbouring sites that both of them are undisturbed;

_{uu}*u*, undisturbed site). Otherwise their states would be changed as

*i*–0, 0–

*j*or 0–0, corresponding with one or both of

*i–j*sites being disturbed, because individuals within the disturbed area would go extinct immediately after the pulse. Thus, we can derive 2.4representing the probability that in a randomly chosen pair of

*i–j*sites, both of them are undisturbed at a disturbance moment (cf. [47]).

Subsequently, the following constraints apply in a disturbed landscape with only two site states—disturbed (*d*) and undisturbed (*u*):
2.5

Thus, we obtain 2.6

Combined with equations (2.3–2.6), the expression of ultimately can be written as 2.7

Based on equation (2.7), we can derive 2.8

Therefore, we can obtain the dynamics of species clumping facing periodic pulse disturbance as shown in equation (S12) (electronic supplementary material, ‘Supplementary methods’).

As a result, equation (2.3) and electronic supplementary material, equation S12 systematically describe the dynamics of population density and spatial pattern for both species 1 and 2 by involving periodic pulse disturbance.

### (c) Numerical simulations

Because of model complexity, we use numerical routines to solve the above system of ordinary differential equations in order to assess the impacts of disturbance aggregation on competitive outcomes. In each simulation, both local and global dispersers are initialized at low density, and each case is run for enough time to achieve a steady state (using *ρ _{i}* < 0.0001 as species extinction threshold;

*i*= 1 or 2). From a variety of preliminary simulations without introducing disturbance, we first determine the range of species mortality rates while keeping both dispersers' birth rates the same in order to yield species coexistence, and then select the intermediate value within the mortality range as a reference case. Next, we choose suitable parameter values for both intra- and interspecific competition that can yield long-term species coexistence. These parameter combinations are used for all subsequent simulations with spatially correlated disturbance. We also vary the parameter values that are kept constant in the following simulations, and determine that they do not alter our general results and conclusions (e.g. electronic supplementary material, figure S1). With the model, we seek to explore whether and how the interactions between different disturbance aspects (

*ρ*,

_{d}*ρ*and

_{d/d}*T*) modify the species coexistence patterns that are shaped in the absence of disturbance.

## 3. Results

We first tested how disturbance extent (*ρ _{d}*) interacts with its spatial correlation (

*q*) in modifying competitive outcomes, while simultaneously varying the mortality rate of global dispersers (

_{d/d}*m*

_{2}) at fixed disturbance periodicity (

*T*= 75) (figure 2). As expected, increasing

*m*

_{2}promoted the persistence of local dispersers by weakening the growth rate of global dispersers. However, both disturbance extent and spatial aggregation largely altered the competitive outcome. As shown in figure 2

*a–c*, the local and global dispersers could coexist when disturbance was zero (i.e.

*ρ*= 0 and

_{d}*q*= 0). After introducing disturbance, an increase of the disturbed area resulted in species exclusion and ultimately the extinction of both species. Spatial correlation in disturbance had almost no influence on competitive outcomes at both very high and low disturbance extent, but largely changed the competitive outcomes at intermediate disturbance extent. In general, spatial aggregation in disturbance rapidly excluded local dispersers through limiting their growth, therefore favouring global dispersal. By contrast, spatial over-dispersion in disturbance enhanced the persistence of local dispersers, which could in this case compensate their weakness in dispersal range/distance, thus driving the globally dispersing species to extinction (figure 2

_{d/d}*c*). As a result, coexistence between local and global dispersers was possible at intermediate spatial aggregation under moderate disturbance. When further increasing

*m*

_{2}to 0.014 (figure 2

*d*), local dispersers now outcompeted global dispersers in the absence of disturbance (i.e.

*ρ*= 0 and

_{d}*q*= 0), based on their relatively higher survival rates. After involving disturbance, the spatial correlation interestingly led to species coexistence, and even completely altered the competitive outcome that was observed in a constant environment: global dispersers now outcompeted local dispersers at high

_{d/d}*q*. On the other hand, spatial correlation in disturbance did not influence the extinction threshold of global dispersers (defined as the critical disturbed area above which a species just goes extinct) after excluding local dispersers (figure 2

_{d/d}*a,b*), while it promoted the extinction risk of local dispersers which tolerated less disturbed area at higher spatial correlation (figure 2

*d*).

We then examined the interactive effects of disturbance periodicity (*T*, i.e. time between disturbance events), extent and spatial correlation on species coexistence (figure 3). As expected, increasing the time interval between disturbances increased the region of species coexistence as species have more time to recover to their initial state during each disturbance interval, though the increment became smaller at larger periodicity. As a consequence, both dispersers tolerated more disturbed area with increasing periodicity, hence delaying their competitive exclusion. Similar to figure 2, both disturbance extent and spatial correlation modified species coexistence. Increasing the disturbed area increased species extinction. Spatial aggregation in disturbance strongly favoured long-range dispersers through suppressing the growth of local dispersers, while spatially uncorrelated disturbance selected for distance-limited dispersers. The selection of different dispersers depending on the spatial pattern of disturbance resulted in the possibility for species to coexist at the combination of moderate disturbance with intermediate spatial aggregation.

## 4. Discussion

We developed a new PA model to explore long-term outcomes in a two-species system subject to spatially correlated disturbances at varying periodicities. An important novelty of our model is that conditional probabilities were applied to characterize the spatial correlation of disturbance, a factor which is often overlooked, despite the ubiquity of spatial correlation in real-world disturbance processes. Although considering both neighbouring interactions (species competition and local dispersal) and disturbance aggregation adds complexity, the PA model allows one to identify the effects of disturbance extent and its spatial aggregation on species coexistence. Despite the limitation of PA that spatial correlation beyond nearest neighbours is ignored (e.g. triplet neighbouring correlation; see electronic supplementary material, ‘Supplementary methods’), it yielded competitive outcomes that are generally similar to those of spatially explicit IBMs [24,25,27]. For instance, spatially correlated disturbance favoured long-range dispersers [24], and species coexistence was sensitive to spatial correlation in disturbance [25]. Investigating spatial effects of disturbance on species coexistence with a PA model can be more efficient than using spatially explicit simulations as in IBMs, for two reasons: (i) IBMs require large amounts of computation time and space, as well as hundreds of replicates for each case, to obtain approximate results; and (ii) IBMs produce stochastic fluctuations and, therefore, large deviations. By incorporating the essential spatial aspects, the PA approach can thus bridge the gap between non-spatial models [12–15] (e.g. MFA) and spatially explicit IBMs.

In our two-species system, spatial correlation of disturbance favoured the long-range dispersers, while the local dispersers were greatly suppressed (figures 2 and 3). This effect is rooted in the longer time required by local dispersers to reach and recolonize the centres of large openings, relative to global dispersers with random establishment. Essentially, by creating clustered openings, increasing local correlation in dynamically disturbed systems increases the temporal variability of population dynamics (i.e. reducing population growth rate [48]). The stronger population growth rate of long-range dispersers disappeared under spatially uncorrelated disturbance, where local dispersal allowed populations to quickly exploit the surrounding small gaps (e.g. via clonal growth). From an evolutionary perspective, our model indirectly suggests that the spatial aggregation of disturbance selects for longer-range dispersal, or more generally, that the disturbance pattern selects the optimal dispersal distance [24]. Similarly, spatially correlated disturbance would favour species with *r*-strategy, as *r*-selected species are strong colonizers that easily exploit newly available habitat, while *K*-selected species can outcompete the more transient colonizers but are slow to expand their range. However, spatial correlation in disturbance did not promote the persistence of global dispersers, as their extinction thresholds (i.e. the maximally tolerated disturbed area above which a species just goes extinct) were not influenced after they had excluded local dispersers (figure 2*a,b*). In contrast with global dispersers, the extinction risk of locally dispersing species was aggravated (i.e. they tolerated much less disturbed area) when disturbance was aggregated (figure 2*c,d*). This indicates that models based on randomly distributed disturbances [9] might underestimate real species extinction rates, as many disturbances exhibit some degree of spatial aggregation. Furthermore, interspecific competition between local and global dispersers can modify both their extinction thresholds (figure 2*c,d*), which could have ramifications for management of plant species even with sympatric species within the same region. Note that, compared with the spatial aggregation of disturbance, the spatial aggregation of habitat fragments has completely opposite (i.e. positive) effects on species persistence, because connected habitats provide more opportunities for local recolonization [33,34,49].

Spatial correlation in disturbance led to different competitive outcomes, depending on how much area was disturbed (figures 2 and 3). When the disturbance extent was very small, then spatial correlation had almost no influence, obviously because species can in this case quickly recover to the initial steady state in between two disturbances. When the disturbed area was very large, both local and global dispersers went extinct, likewise nullifying the impact of spatial correlation. Only at intermediate disturbance extent did spatial correlation strongly alter the competitive outcomes, with high correlation favouring global dispersers and low correlation selecting for local dispersers (figures 2*c* and 3). As a consequence, species exclusion would accelerate at both very high and very low spatial correlation in disturbance, making species coexistence possible only at intermediate levels. These results contradict the earlier conclusion of Moloney & Levin [50] that spatial aggregation of disturbance only plays a minor role in long-term species coexistence. Possibly, the three interacting plant species used in their manipulated grassland experiments did not significantly differ in dispersal distance/range, while interspecific variation in dispersal ability was a key factor in the current model. On the other hand, disturbance periodicity did not change the general competitive patterns shaped by a combination of disturbance extent and spatial aggregation (figure 3; electronic supplementary material, figure S1). Yet increasing periodicity enlarged the region of initial competition patterns formed in the absence of disturbance (i.e. *ρ _{d}*

*=*

*q*= 0), and reduced the region where both species went extinct. In other words, the location where effects of spatial correlation on competitive outcomes were strong in the diagram was displaced towards higher disturbance extent as periodicity increased. From an evolutionary perspective, species with

_{d/d}*r*-strategy would be favoured if disturbance frequency is high, because of their strong reproduction and colonization ability.

Interestingly, a high spatial correlation in disturbance promoted long-term species coexistence, but only in case the local dispersers competitively excluded the global dispersers in the absence of disturbance (figure 2*d*; electronic supplementary material, figure S1). Essentially, spatial correlation in disturbance balances the trade-off between dispersal ability and mortality rate. As global dispersers are superior in dispersal range/distance relative to local dispersers, their exclusion was realized by weakening the survival rate (hence local dispersers become superior species). In this case, introducing spatially correlated disturbance reduces the growth rate of local dispersers, which consequently compensates for the mortality disadvantage of global dispersers to some extent, ultimately leading to species coexistence. Using a spatially explicit IBM, Banitz *et al.* [25] proposed that over-dispersed disturbances lead to species exclusion in species-rich communities. However, in our model, it depends on the competitive outcomes in a constant environment. For instance, if two different dispersers coexist when disturbance is absent (figure 2*a,b*), then spatially over-dispersed disturbance can maintain the coexistence pattern (as opposed to aggregated disturbance) by promoting local recolonization opportunities for short-range dispersers. However, we stress that stable species coexistence is less likely by introducing spatially structured disturbance, if global dispersers outcompete local dispersers in a constant environment. Three effects explain this in our model: (i) spatial aggregation in disturbance had almost no influence on the extinction thresholds of global dispersers; (ii) the colonization ability of global dispersers prevailed because of random establishment; and (iii) spatially correlated disturbance aggravated the extinction risk of local dispersers.

Our PA model is most applicable to pulse perturbations that occur nearly instantaneously—at least relative to the lifespan of individuals—such as extreme events (e.g. fire, flooding, extreme drought, pest outbreaks). So-called ‘press’ disturbances that gradually deteriorate the environment (e.g. persistent pollutants) would not qualify, because we assumed no decline in habitat suitability for species recolonization. For model simplicity and mathematical tractability, we only considered two extremes of dispersal—local dispersal to nearest neighbouring sites and global dispersal with random establishment across the entire landscape. This assumption is relatively restrictive as most species in nature exhibit a specific dispersal kernel (e.g. Gaussian kernel). However, our results should be considered as representative for two extreme dispersal modes that delimit the continuum in between, because both types of dispersal have long been used to simplify ecological models, yielding approximately the same results as more realistic dispersal [33,34,51–55].

In this study, we elaborated a mechanistic framework to demonstrate that the spatial structure of disturbance has the potential to alter competitive outcomes. Based on a single population model, Johst & Drechsler [20] proposed that lowering spatial correlation in disturbance may cause a species to benefit, or at least not to suffer. However, application to multispecies communities should be executed with caution, as reducing spatial correlation of disturbance in our model can accelerate competitive exclusion (figure 2*c,d*), thereby decreasing species diversity. With respect to conservation, our results suggest that introducing spatial disturbance to an ecosystem or altering the aggregation of an existing disturbance regime (e.g. implementing prescribed burning, human-controlled flooding, grazing or mowing) could be a useful management tool either to promote species coexistence or to control species invasion. For instance, the current model might apply to the Konza Prairie ecosystem in Kansas, USA, where most of the native species exhibit vegetative reproduction (e.g. big bluestem), while many invasive species (e.g. nodding thistle) have very long-range seed dispersal [56]. In this case, our model predicted that introducing spatially over-dispersed dynamic disturbance (via controlled fire) would promote the local vegetative species but suppress the widely dispersed invaders. Further study could thus focus on model validation in systems such as these. Moloney & Levin [50] have already explored the impacts of the spatial and temporal structure of disturbance on the population dynamics of three interacting plant species reproducing mainly via seed dispersal in an annual, serpentine grassland. Similarly, Zhang & Shea [57] experimentally manipulated the temporal autocorrelation of disturbance occurrence to impede the invasion of the highly dispersive thistle, as theoretically confirmed by Garrison *et al.* [58]. Extending this type of experiment to test the role of spatial correlation of disturbance in the coexistence between two significantly different dispersers (one with seed dispersal and another with clonal growth) would be fairly straightforward. Overall, analysing spatial disturbance patterns with PA modelling can provide new insights into species coexistence, and therefore biodiversity conservation.

## Authors' contributions

J.L. and Z.L. conceived the study; J.L. and Z.Y. conducted simulations; J.L. performed the analysis; J.L., D.A.W, A.D.M and I.N. discussed the results; J.L. wrote the manuscript; and D.A.W, A.D.M and I.N. revised the manuscript.

## Competing interests

We have no competing financial interests.

## Funding

Z.L. was supported by the National Key Basic Research Program of China (no. 2013CB429905) and the National Natural Science Foundation of China (no. 41571505). A.D.M. was supported in part by an NSF grant (no. DEB-1353301).

## Acknowledgement

This paper is contribution no. 72 of the Central Michigan University Institute for Great Lakes Research.

- Received March 8, 2016.
- Accepted April 11, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.