Habitat loss and isolation caused by landscape fragmentation represent a growing threat to global biodiversity. Existing theory suggests that the process will lead to a decline in metapopulation viability. However, since most metapopulation models are restricted to simple networks of discrete habitat patches, the effects of real landscape fragmentation, particularly in stochastic environments, are not well understood. To close this major gap in ecological theory, we developed a spatially explicit, individual-based model applicable to realistic landscape structures, bridging metapopulation ecology and landscape ecology. This model reproduced classical metapopulation dynamics under conventional model assumptions, but on fragmented landscapes, it uncovered general dynamics that are in stark contradiction to the prevailing views in the ecological and conservation literature. Notably, fragmentation can give rise to a series of dualities: a) positive and negative responses to environmental noise, b) relative slowdown and acceleration in density decline, and c) synchronization and desynchronization of local population dynamics. Furthermore, counter to common intuition, species that interact locally (“residents”) were often more resilient to fragmentation than long-ranging “migrants”. This set of findings signals a need to fundamentally reconsider our approach to ecosystem management in a noisy and fragmented world.

The model landscapes were generated by sampling Gaussian random fields with zero mean and a variance-covariance structure with specified parameters. 

Individual-based model:

On each simulated landscape, semelparous individuals mature during every non-overlapping generation and each produces a Poisson-distributed number of propagules, after which the sessile adult dies instantaneously. The propagules disperse according to a bivariate Gaussian kernel with mean zero and a variance-covariance matrix. Propagules that disperse into the matrix die immediately and do not influence local competition. Those that disperse into habitat experience density-dependent competition with other recruits, influencing their probability of becoming established as members of the next generation. We estimated the local propagule densities by convolving the distribution of propagules with a competition kernel. The code for simulating all individual-level processes is vectorized to improve runtimes.

To model regional stochasticity, we assumed that the fecundity rate varies over space and time. Variation was assumed to be independent among generations (no temporal autocorrelation). Fecundity was spatially correlated, such that it is multivariate log-normally distributed. We randomized the spatial distribution of habitat quality over time by resampling a Gaussian random field for fecundity rate at the start of each generation.

In our first simulation script (sim_persistence), we investigate how the classical relationships between metapopulation persistence time (i.e., mean-time to global extinction) and model parameters could be altered on fragmented landscapes. In this analysis, we sample square subregions of a given area from each pre-generated landscape 1, reset the periodic boundary condition, and simulate metapopulation dynamics until no individual remains or the generation limit is reached. This part of our analysis assumes globally synchronized environmental stochasticity, which we describe by randomizing a global value of fecundity rate per generation. In the case of mass-action mixing, propagules disperse uniformly, and offspring competitions are influenced by global rather than local density. We also test the effects of habitat aggregation and, by extension, fragmentation per se by collecting the habitat fragments in each previously sampled subregion into a square patch and rerunning the simulations.

In our second simulation script (sim_abundance), we explore how fragmentation mediates metapopulation abundance (i.e., global population density). Here, each simulation continues for a fixed number of generations long enough for metapopulation abundance to reach quasi-stationarity. Two values of mean fecundity rate are used to represent a contrast in reproductive output: on the habitat-poor landscapes, one gives rise to global persistence and the other to global extinction. We then compare our abundance predictions to results from the stochastic mean-field (Beverton-Holt) model.

In our third simulation script (sim_synchrony), we estimate metapopulation synchrony in the forms of spatial correlations in local occupancy and local abundance, conditional on persistence. Each landscape is divided into a grid of square sampling sites; neighboring sites are separated by one lattice spacing. To measure local occupancy, we determine the fraction of habitat cells that are occupied on sampling sites containing habitat. To measure local abundance, we calculate the number of occupants divided by the number of habitat cells on nonempty sites. Since any inherent spatial synchrony in the system can be confounded if, at the time of measurement, a vast majority of sampling sites is either vacant or fully occupied, we adjust the mean fecundity rates to keep mean global population size always within a standard, ecologically informative range.