Metapopulation capacity determines food chain length in fragmented landscapes

Edited by Nils Chr. Stenseth, Universitetet i Oslo, Oslo, Norway, and approved July 8, 2021 (received for review February 9, 2021)
August 20, 2021
118 (34) e2102733118


Understanding the persistence of populations in fragmented landscapes is critical for predicting the consequences of habitat destruction, yet analytical tools are largely lacking. Metapopulation capacity provides one such tool, because it summarizes the influences of habitat area and distribution on population persistence in a single metric. However, surprisingly few efforts have extended this theory to multispecies communities. Our analyses demonstrate the power of metapopulation capacity theory in predicting the persistence of prey–predator pairs and food chains in heterogeneous, fragmented landscapes. Such analytic insights serve as a benchmark to predict the consequences of habitat changes. Our findings thus have broad implications for both ecological research and conservation practices.


Metapopulation capacity provides an analytic tool to quantify the impact of landscape configuration on metapopulation persistence, which has proven powerful in biological conservation. Yet surprisingly few efforts have been made to apply this approach to multispecies systems. Here, we extend metapopulation capacity theory to predict the persistence of trophically interacting species. Our results demonstrate that metapopulation capacity could be used to predict the persistence of trophic systems such as prey–predator pairs and food chains in fragmented landscapes. In particular, we derive explicit predictions for food chain length as a function of metapopulation capacity, top-down control, and population dynamical parameters. Under certain assumptions, we show that the fraction of empty patches for the basal species provides a useful indicator to predict the length of food chains that a fragmented landscape can support and confirm this prediction for a host–parasitoid interaction. We further show that the impact of habitat changes on biodiversity can be predicted from changes in metapopulation capacity or approximately by changes in the fraction of empty patches. Our study provides an important step toward a spatially explicit theory of trophic metacommunities and a useful tool for predicting their responses to habitat changes.
Spatial models offer a critical tool for understanding the maintenance of biodiversity across scales and guiding conservation efforts (1, 2). More than half a century ago, two theoretical frameworks were developed to understand species persistence and diversity in a spatial context, namely, the theory of island biogeography (3) and metapopulation theory (4). These two theories laid the foundation for a new scientific paradigm and have motivated decades of research into the role of space in population dynamics and species interactions (1, 2, 5). While subsequent studies during the 1970s largely focused on the theory of island biogeography, metapopulation theory has been increasingly adopted since the late 1980s because species extinctions in fragmented landscapes have become an important issue (68), particularly for species at high trophic levels (9, 10).
The original metapopulation model highlighted a dynamical equilibrium perspective of species persistence in patchy landscapes: populations undergo continuing local extinction, but recolonization maintains the persistence of the whole metapopulation (4). Despite its simplicity, this model has generated significant insights that are of both theoretical and practical implications. In particular, the model predicts that for a species living in a patchy landscape, the minimum habitat area required for the species to persist equals the observed fraction of empty patches (referred to as Levins’ rule; ref. 8). Such a rule provides a useful practical guide for landscape conservation, which has been validated by some empirical studies (e.g., on the conservation of the northern spotted owl; ref. 11). However, Levins’ model missed two important aspects, namely, species interactions and landscape configuration, which potentially restricts its applications in real-world ecosystems. Consequently, a significant body of work in the past decades has been devoted to extending Levins’ model by incorporating interspecific interactions and explicit landscape configurations (2, 8, 12).
Soon after the inaugural paper that proposed the metapopulation concept (4), Levins and Culver (13) extended metapopulation theory to study the coexistence of two competing species in fragmented landscapes. Later studies further extended the model to a broader context of interspecific relationships, including prey–predator, host–parasitoid, and mutualistic interactions (1416), and more complex species assemblages including multiple competing species (17, 18) and complex food webs (19, 20). These studies have contributed to a new research framework, namely, metacommunity ecology, which has received considerable attention and development during the past two decades (2, 21, 22).
Early metapopulation models were spatially implicit and thus did not account for the area and spatial distribution of habitat patches. This simplification made the models mathematically tractable but prevented them from incorporating real-world landscape heterogeneity (8). Different types of spatially explicit models have been developed to understand the effects of dispersal and landscape configuration on population dynamics and metacommunity patterns (1, 2, 8, 23). While these models often require a simulation approach, the spatially explicit patch-occupancy model of metapopulations (24) provides a remarkable exception. Specifically, Hanski and Ovaskainen (24) formulated the concept of “metapopulation capacity,” which integrates all the information of habitat area and spatial distribution that is relevant to species persistence into a single metric. Mathematically, metapopulation capacity is calculated as the dominant eigenvalue of the “landscape” matrix, the elements of which are determined by species dispersal, habitat area, and arrangement. This analytic approach provides an elegant tool to evaluate species persistence in fragmented, heterogeneous landscapes (24). The concept of metapopulation capacity has thus been widely adopted in empirical studies and conservation biology (2527), and it has motivated new theoretical developments [e.g., predicting how patch network properties and dispersal regimes affect metapopulation capacity and thus species persistence (2830)].
Despite these theoretical advances, surprisingly little progress has been made to extend metapopulation capacity theory to multispecies systems. Based on observational data, van Nouhuys and Hanski (31) showed that as the metapopulation capacity of the host species increased, the occupancy of its parasitoid also increased. Hanski (32) extended their spatially explicit patch-occupancy model to understand the coexistence of competitors in fragmented landscapes but with simulation approaches. It remains unexplored whether the analytical tool of metapopulation capacity can be applied to understand the persistence of interacting species in fragmented landscapes. Such extensions, if applicable, would overcome limitations of Levins’ model and offer novel analytic insights for predicting the responses of community diversity to habitat changes in realistic settings (33, 34).
In this paper, we develop a spatially explicit patch-occupancy model of trophic interactions and examine whether metapopulation capacity can predict the persistence of trophically interacting species. We start with a prey–predator system, for which we derive approximate solutions and show that metapopulation capacity does predict its persistence in fragmented landscapes. We then extend our model to food chains and derive an analytic solution for food chain length as a function of metapopulation capacity, strength of top-down control, and population dynamical parameters. Under certain assumptions, we show that the fraction of empty patches provides a useful indicator for predicting the length of food chains that a fragmented landscape can support, which extends Levins’ rule. We show that our theoretical results shed light on the empirical system described in ref. 31. Our study demonstrates the value of metapopulation capacity theory in understanding the persistence of multitrophic systems in fragmented landscapes and predicting their responses to habitat changes.


A Spatially Explicit Patch-Occupancy Model of Prey–Predator Interactions.

We first develop a spatially explicit model for prey–predator dynamics by extending Hanski and Ovaskainen’s (24) single-species metapopulation model. In our model, the temporal dynamics of probabilities that patch i is occupied by the prey (p1i) and predator (p2i) follow
where Aj is the area of patch j and dij is the distance between patch i and j. c1 and e1 are the prey colonization and extinction parameters, respectively, and α1 captures how fast the colonization rate decreases with distance for the prey. c2, e2, and α2 are similar parameters for the predator. Following ref. 24, we assume that colonization rate increases linearly with patch area and that local extinction rate decreases inversely with patch area. To characterize bottom-up effects, we assume that the predator can recolonize only patches in which the prey is present and the predator is absent (i.e., p1ip2i) and the local extinction of the predator can be caused by either prey or predator extinction. To characterize top-down effects, we assume that the prey has an extinction rate of e1 (in patches with unit area) in the absence of the predator and e1f in the presence of the predator. Thus, the extinction rate for the prey is e11p2ip1i+e1fp2ip1i, where p2ip1i gives the conditional probability of predator’s presence on the patch i provided that the prey is present (i.e., “conditional incidence”; ref. 14). Presence of the predator can either increase (f > 1) or decrease (0 < f < 1) the extinction rate of its prey or have no effect in donor control cases (f = 1) (ref. 14; see SI Appendix, Appendix S1 for discussion).
In this model, the prey and predator may both persist, only the prey can persist, or both go extinct. For instance, when the local extinction rates of the prey and predator are low, both species can persist on the landscape with high occupancies (Fig. 1). As their extinction rates increase, the average occupancies of both species decrease. At certain extinction rates, the predator will go extinct, followed by the extinction of the prey (Fig. 1). Below, we provide the persistence conditions for the prey and predator (Eqs. 2 and 3).
Fig. 1.
The equilibrium average occupancies of prey and predator species along a gradient of local extinction parameter, on a landscape with 20 patches. For simplicity, we assume that prey and predator have same extinction parameter (e). Solid and dashed lines show results from simulation and analytic approximation, respectively. The three upper panels exhibit the distribution of habitat patches and species occupancy corresponding to three values of extinction parameters (e), indicated by the dashed gray lines. The habitat patch is indicated by green points, with size proportional to patch area. The size of blue and red points (relative to green points) indicates the occupancy probability of prey and predator, respectively. See SI Appendix, Table S1 for parameter values.
Because the persistence conditions characterize the thresholds below which the prey or predator species goes extinct, we are essentially interested in the dynamics of Eq. 1 at the brink of species extinction. When the prey species is close to extinction (p1i0 for all i), we expect the predator species to have already gone extinct (i.e., p2i=0 for all i) (e.g., Fig. 1). In such cases, the system reduces to the single prey species, and thus the persistence condition for the prey species is given by (SI Appendix, Appendix S1; ref. 24)
where λ1 denotes the metapopulation capacity of the prey and is calculated as the dominant eigenvalue of the landscape matrix for the prey (see SI Appendix, Appendix S1). When the predator is close to extinction, it has low occupancy probabilities in all patches (i.e., p2i0for alli). In such cases, the predator has weak top-down effects, and thus the prey has similar equilibrium occupancy as in single-species metapopulations. Under a homogeneity assumption (i.e., the prey has the same occupancy probability in all patches), we can derive the persistence condition for the predator as (SI Appendix, Appendix S1)
which is more stringent than Eq. 2a. Here λ2 denotes the metapopulation capacity of the predator. These conditions suggest that top-down control does not affect the prey persistence but modulates predator persistence. SI Appendix, Appendix S1 provides approximate solutions for the average occupancy for both prey and predator.
The above solutions (Eq. 2) are derived under the assumption that the predator goes extinct before the prey, which may often be the case (e.g., refs. 9 and 10). However, if the predator species has a very high or even infinite colonization rate (e.g., c2), the predator should always coexist with or go extinct together with the prey. In such cases, the predator reaches its equilibrium only when p1i=p2i for all i, and we can derive the condition of persistence for both prey and predator species (SI Appendix, Appendix S1):
We conduct numerical simulations to test the accuracy of our analytic approximations of persistence conditions (summarized in Table 1). Simulation results showed that, although derived under the homogeneity assumption, our approximate solutions work remarkably well in predicting prey–predator persistence in heterogeneous landscapes (SI Appendix, Fig. S1). The discrepancy in the threshold extinction rate between our analytic approximation and numerical simulations is generally below 5%. Such discrepancy increases as the strength of top-down control (f) increases, but it does not change with the spatial heterogeneity of patch areas (SI Appendix, Fig. S1). Moreover, our simulations showed that, as the colonization rate of the predator increased from small to large values, the threshold metapopulation capacities for prey and predator persistence exhibited a transition from Eqs. 2 to 3 (SI Appendix, Fig. S2). Interestingly, we found that when the predators substantially reduced the prey extinction rate (i.e., f << 1), the transition could take different pathways for systems starting with different initial predator occupancies (SI Appendix, Fig. S2). In other words, the system can exhibit alternative stable states at intermediate rates of predator colonization (i.e., the prey and predator either both persist or both go extinct) depending on the initial predator occupancy (SI Appendix, Fig. S3 and Appendix S1). Lastly, while our model describes patch-occupancy dynamics in a deterministic manner, simulations of their stochastic counterparts (SI Appendix, Appendix S2) showed that our analytic solutions provide reasonable predictions of the persistence conditions of the prey and predator in stochastic settings, although the average occupancies from stochastic simulations were generally lower than those predicted by our deterministic model (SI Appendix, Fig. S4).
Table 1.
Analytic approximations on the persistence condition of prey–predator systems and the maximum food chain length in spatially explicit metacommunity models
 Finite colonization of predator(s)Infinite colonization of predator(s)
 Species-specific parametersSpecies shared parameters 
 Top-down controlDonor control (f = 1) 
Prey persistence conditione1c1λ1<1ecλ<1U1<1e1fc1λ1<1
Predator persistence conditione1c1λ1+e1f+e2c2λ2<1e(2+f)cλ<13U1<1e1fc1λ1<1
Food chain length (L)k=1Li=1k1eifi+ekckλk<1L<12f8λfce+2f22fL<12(8U1+11)Unlimited if e1fc1λ1<1
ck, ek, and λk represent the colonization rate, local extinction rate, and metapopulation capacities of the species at trophic level k, respectively. fk captures the top-down effect of the trophic level k+1 on the trophic level k. U1 = e/λc represents the average fraction of empty patches for the basal species.

Extensions to Food Chains.

We further extend our spatially explicit model to longer food chains that contain L trophic levels:
where pki denotes the occupancy probability of trophic level k on patch i. The parameters ck, ek, and αk have similar meanings as in Eq. 1. fk captures the top-down effect of the trophic level k+1 on extinction of the trophic level k in a food chain of length L. Because of local dynamics such as trophic cascades, fk could vary with food chain length (e.g., an herbivore might be less able to overexploit its plant resource) if a natural enemy keeps its numbers in check. Based on Eq. 4, we can similarly derive the persistence condition for a food chain of length L (SI Appendix, Appendix S1):
This equation shows that the persistence of the food chain in fragmented landscapes is jointly determined by the metapopulation capacities of all species (λi), the species-specific parameters of colonization (ci) and extinction (ei), and the strength of top-down effects of higher trophic levels on their prey levels (fi). If all species have the same metapopulation capacity (λiλ,i) and dynamical parameters (i.e., e, c, and f), we can derive the upper limit of the food chain length that a landscape can support (SI Appendix, Appendix S1):
The above analytic solution predicts that the maximum food chain length that a landscape can support increases with the metapopulation capacity (λ) and colonization rate (c), but it decreases with the strength of top-down control (f) and the extinction rate (e) (Fig. 2A).
Fig. 2.
(A) The maximum food chain length as a function of metapopulation capacity (λ) and strength of top-down control (f), as predicted by Eq. 6. (B and C) The realized food chain length (points) as functions of the theoretical expectation on the maximum food chain length (L^, Eq. 6) (B) and the average fraction of empty patches (U1) (C). The red lines show the theoretical expected food chain length by Eq. 6 in B and Eq. 7 in C. The blue and green points indicate, respectively, simulated food chains whose realized food chain are higher and lower than the theoretical expectation. The inserted panels in B and C show the relationship of deviation between predicted and realized food chain length with the strength of top-down effects. See SI Appendix, Table S1 for parameter values.
In donor-control systems (f = 1), the above solution (Eq. 6) can be further simplified:
where U1=eλc represents the average fraction of empty patches for the basal species (see SI Appendix, Appendix S1; ref. 24). Eq. 7 suggests that one can use the fraction of empty patches of the basal species to predict the number of trophic levels that a landscape can support. Specifically, the basal species persists (Ld^>1) when U1<1, and a second trophic level can persist (Ld^>2) when U1<1/3, a third trophic level can persist when U1<1/6, and a fourth trophic level can persist when U1<1/10.
Again, we note that the above solutions (Eqs. 57) are derived by assuming finite colonization rates of higher trophic levels. If all predators (i.e., trophic level 2) have infinite colonization rates, the food chain length is unlimited as long as the basal species satisfies the condition given by Eq. 3 (Table 1).
Our simulations confirmed the validity of the solutions given by Eqs. 57. Among 1,000 simulated food chains, more than 90% reached a food chain length as predicted by Eq. 6, while the exceptions exhibited one trophic level higher or lower than predicted (blue and green points, respectively, in Fig. 2B and SI Appendix, Fig. S5). Although it ignores top-down control and assumes identical species parameters, Eq. 7 provided an overall good prediction of realized food chain length for systems with top-down control (Fig. 2C) and species-specific parameters (SI Appendix, Fig. S5). Eqs. 6 and 7, however, tended to overestimate food chain length when f > 1 and underestimate it when f < 1 (inserted panels in Fig. 2 B and C).

Persistence of Trophic Species under Habitat Changes.

The above solutions (Eqs. 2, 3, and 57) clarify that the persistence of prey–predator or food chain systems depends on habitat area and spatial distribution through the metapopulation capacity. Consequently, the effect of habitat changes can be understood from their effects on the metapopulation capacities of all species. Here, we consider three scenarios of habitat changes: 1) habitat deterioration, in which the area of each patch is reduced to a proportion h of its original area (i.e., Ai'=hAi). In this case, we have λ'λ=h2, where λ' and λ represent the new and original metapopulation capacities, respectively (SI Appendix, Appendix S1). 2) Habitat loss, in which a fraction of randomly selected patches is lost, and the remaining fraction is a*100%. In this case, we can derive the change in metapopulation capacity if all patches have the same area and connectivity: λ'λ=na1n1, where n is the initial number of patches (SI Appendix, Appendix S1). 3) Habitat fragmentation, in which the connectivity decreased due to an increased effective distance between all patches (dij'=dij+s, s > 0). In this case, we have λ'λ=eαs (SI Appendix, Appendix S1).
We denote Δλλ'λ, which equals h2,na1n1,andeαs, respectively, under the above three scenarios of habitat changes (i.e., habitat deterioration, loss, and fragmentation). By substituting Δλ into Eqs. 57, we can evaluate how habitat changes will alter the maximum food chain length. In particular, when all species have same parameters, we can substitute λ'=λΔλ to Eq. 6 to predict the food chain length in the changed landscape (L'):
Moreover, in the donor-control case, the fraction of empty patches following habitat changes becomes: U1'=eλ'c=U1Δλ. We then substitute it to Eq. 7 and have
Thus, based on the initially observed fraction of empty patches (U1) and the proportional change in the metapopulation capacity (Δλ) following habitat changes, one can predict the change in food chain length (SI Appendix, Appendix S1).
Our simulations showed that all three types of habitat changes decreased the persistence of higher trophic levels and reduced food chain length (Fig. 3). The change in food chain length was well predicted by our analytic approximation Eq. 8 (Fig. 3). Furthermore, our simulations showed that Eq. 9, though derived under the donor-control case and requiring no detailed information on population dynamical parameters (e.g., c, e, f) also provided a reasonable predictor for food chain changes (Fig. 3).
Fig. 3.
The effect of habitat changes on food chain length: habitat deterioration (A and D), habitat loss (B and E), and habitat fragmentation (C and F). AC illustrate the three types of habitat changes. In A and B, green and gray points represent the remaining and destroyed habitats, respectively. In C, habitats are not reduced, but the distance between patches increases. DF show results in food chains with different strengths of top-down control: f = 0.5 (red), 1 (black), and 2 (blue). Solid lines represent simulated results, and dashed and dotted lines show the theoretical expectations by Eqs. 8 and 9, respectively. See SI Appendix, Table S1 for parameter values.

Empirical Test.

The butterfly–wasp metacommunity in the Åland archipelago provides a unique system to test our model predictions (31). The Glanville fritillary butterfly Melitaea cinxia has two cooccurring parasitoids (i.e., Cotesia melitaearum and Hyposoter horticola), which have low and high rates of dispersal (or colonization), respectively (31, 35). Given that detailed dynamical parameters are unavailable for the two parasitoids, we tested two model predictions: 1) Predators with limited colonization rates (e.g., C. melitaearum) can only persist on landscapes in which the fraction of empty patches U1 < 1/3 (Eq. 7) and their average occupancy could be predicted by p2=13U1 (SI Appendix, Appendix S1); 2) Predators with very high colonization rates (e.g., H. horticola) could persist as long as the prey persists on the landscape (Eq. 3). Our empirical analyses generally supported these two predictions. Among 102 patch networks, C. melitaearum was absent or had low occupancy (only 4 out of 89 networks had an average occupancy >0.05) in those networks in which the fraction of empty patches for the butterfly (i.e., U1) is larger than 1/3; in contrast, it was present in many networks (8 out of 13 networks had an average occupancy >0.05) in those networks in which U1 < 1/3 and its occupancy followed the trend as our model predicted (Fig. 4). Moreover, H. horticola occurred in all patch networks in which the host butterfly occurred (35), consistent with our prediction under high colonization rates.
Fig. 4.
Average occupancy of the parasitoid C. melitaearum as a function of the average fraction of empty patches for the host butterfly (i.e., one minus the average occupancy) across 102 patch networks in the Åland archipelago. Black points represent patch networks consisting of more than 10 patches and circles with no more than 10 patches. The gray area indicates the region where the parasitoid is predicted to persist (i.e., U1 < 1/3). The red line shows the predicted occupancy (SI Appendix, Eq. S4).


Metapopulation capacity offers an elegant, quantitative tool to assess the effects of landscape configuration (e.g., habitat area and spatial distribution) on species persistence in fragmented, heterogeneous landscapes (24). Our theoretical analyses extend metapopulation capacity theory to predict the persistence of multitrophic systems (e.g., prey–predator interactions and food chains). While our analytic solutions are approximate, simulations demonstrate their remarkable accuracy in predicting persistence of species across trophic levels. Although our model was built upon prey–predator interactions, it has broad implications for antagonistic systems including plant–herbivore, host–parasite, and host–parasitoid relationships. As an example, we use similar approaches and derive persistence conditions for a more complex food web module with omnivory and competition (SI Appendix, Appendix S3). As natural ecosystems are inherently diverse and complex, our multitrophic extension of metapopulation capacity theory offers ample opportunities for empirical applications and further theory development.

Persistence of Trophic Systems in Heterogeneous Landscapes.

The dynamics of interacting species are determined by both abiotic (e.g., landscape configuration) and biotic (e.g., species interactions) conditions they experience. The classic theory of metapopulation capacity captures the effect of abiotic conditions (i.e., habitat area and spatial distribution) on species persistence in a single-species context. In multispecies systems, however, biotic interactions modify the landscapes experienced by any species. For instance, the presence of predators can alter habitat quality and thus local extinction rate for its prey (i.e., top-down effect), whereas the absence of prey makes potential habitats for the predator unavailable (i.e., bottom-up effects). Such modifications add much complexity, making it difficult to analytically predict the persistence of interacting species in heterogeneous landscapes.
Our analytic solutions clarify how landscape configuration (characterized by the metapopulation capacity, λ), population dynamical parameters (e and c), and top-down control (f) jointly determine the persistence of prey–predator interactions and food chains (Table 1). Interestingly, they suggest that if a predator has a finite colonization rate, top-down control (f) impairs the persistence of the predator itself but not that of its prey (Eq. 2). This counterintuitive result can be understood from the weak top-down effect at the brink of species extinction. When the prey species is close to extinction, the predator should have an even lower occupancy or be already extinct and thus have a weak impact on the prey. But we note that although the top-down control by the predator does not affect the persistence condition of the prey, it can significantly decrease the occupancy of the prey when both species are far away from the brink of extinction (Fig. 1). Moreover, if the predator has a very large colonization rate and thus quickly occupies patches wherever its prey occurs, its top-down control can modulate the persistence of prey (Eq. 3). Consequently, when the top-down control increases prey extinction (f > 1), the persistence of the predator could first increase and then decrease with its colonization rate (SI Appendix, Fig. S2). Such a maximum persistence at intermediate colonization rate is in line with previous findings in local prey–predator systems, in which increasing the predators’ attack rate can first increase its own abundance but then decrease it until it reaches extinction (36).
Our model also provides a quantitative tool to predict food chain length. Using spatially homogenous models, Holt (14) predicted that the length of food chains should increase with habitat area and decrease with isolation. Our analytic solutions extend this prediction to heterogeneous landscapes, in which the influence of landscape heterogeneity (e.g., patch areas and spatial distribution) is integrated into metapopulation capacity. Our finding that metapopulation capacity provides an important predictor of species persistence across trophic levels is in line with the recently proposed concept of species-fragmented area relationship (SFAR; ref. 37). But our models extend SFAR by explicitly incorporating trophic interactions and clarifying the effect of top-down control on food chain length (Eq. 6).

An Extended Levins’ Rule.

One interesting insight from classic metapopulation models is Levins’ rule (i.e., the minimum fraction of habitats to be conserved is determined by the fraction of empty patches) (8, 38). Such a rule is particularly useful in practice because the fraction of empty patches can be easily measured, and no detailed knowledge on metapopulation dynamics is needed to predict the consequences of habitat changes (11, 38). Our results extend Levins’ rule to multitrophic systems in heterogeneous landscapes (see SI Appendix, Appendix S4 for similar results in homogeneous landscapes). Specifically, our Eq. 7 provides a quantitative tool to predict food chain length based on the fraction of empty patches for the basal species (U1). As expected, a lower U1 indicates a larger potential for supporting higher trophic levels, and the threshold for increasing one trophic level is given by Eq. 7. The rationale underlying such an extended Levins’ rule is that U1 is an emergent property, which integrates the combined effects of habitat configuration and population dynamical parameters.
Our extended Levins’ rule provides a reasonable explanation for the occupancy of C. melitaearum, a specialist parasitoid of its host butterfly M. cinxia in the Åland archipelago. The butterfly M. cinxia in this system provided the first empirical support for the metapopulation capacity theory (24), which was often studied in a single-species metapopulation context. However, M. cinxia is embedded in a complex species interaction network including parasites and parasitoids (31, 35). A recent study revealed no evidence for a top-down effect of C. melitaearum on the extinction of M. cinxia and showed that they had comparable rates of colonization and extinction (39), which satisfy the assumptions under which our extended Levins’ rule was derived. Consistent with our prediction (Eq. 7 and SI Appendix, Appendix S4), C. melitaearum was found to be absent or have low occupancy in patch networks in which the average fraction of empty patches was larger than 1/3, and it occurred with increasing occupancy as the fraction of empty patches decreased from 1/3 to 0 (Fig. 4). These results support the extended Levins’ rule as a rule of thumb in empirical studies.
However, we make two caveats to our extended Levins’ rule. First, our prediction (Eq. 7) was derived for donor-control interactions and under the assumption that all species have the same dynamical parameters. In the presence of top-down control and variation in species parameters, our numerical simulations suggest that Eq. 7 can either over- or underestimate food chain length (Fig. 2C and SI Appendix, Fig. S5). Moreover, in natural ecosystems, U1 can also be affected by processes that are ignored in our models, for instance, the rescue effect and transient dynamics (15), which may bias our prediction.

Conservation Implications.

Our theoretical results have useful implications for landscape conservation. To determine the critical area or connectivity for food chain persistence, we need to understand how habitat changes affect metapopulation capacity across trophic levels and compare them with their persistence conditions (Eq. 8). These persistence conditions are in turn determined by species-specific parameters of population dynamics (αi, ei, and ci) and trophic interactions (fi), yet estimation of these parameters requires long-term observational data (8) that are not always accessible. Alternatively, our extended Levins’ rule provides a practical way to predict food chain changes based on changes in the observed fraction of empty patches (Eq. 9). For instance, start with a landscape in which the basal species occurs in 90% of all habitat patches (i.e., U110%). Eq. 9 predicts that such a landscape can support up to four trophic levels, but any degree of habitat destruction or fragmentation will cause the extinction of the fourth trophic level. Continuing habitat destruction can further drive the third trophic level to extinction when U1'>1/6 (or Δλ<0.6), which happens when all patches shrink by more than 23% (i.e., h < 0.77), more than 40% patches are lost (i.e., a < 0.6), or the increased distance between patches exceeds half of the characteristic distance (i.e., s<0.5/α) (Eq. 9). Again, we note that Eq. 9 should be used with caution because it ignores top-down effects and variation in species dynamical parameters. Nevertheless, our theoretical predictions are qualitatively consistent with numerical simulations as well as simulations of population dynamics under habitat changes that integrate both bottom-up and top-down effects to show that complex food webs are deconstructed from the top to the bottom (40).
Conservation efforts should thus be targeted at enhancing metapopulation capacity, which can be achieved by increasing the area, number, and/or connectivity of habitat patches (Fig. 3). Conversely, if the objective is to eliminate a parasite species, landscape management should be designed to decrease metapopulation capacity to such an extent that the parasite goes extinct but the host species persists. In this case, characteristics of the parasite species should be taken into account. For instance, if the parasite species has a very high colonization rate (e.g., H. horticola in the M. cinxia system), elimination of such species by reducing metapopulation capacity risks host extinction as well (Eq. 3).


Our study provides an important step toward understanding and predicting trophic interactions in heterogeneous, fragmented landscapes. By extending metapopulation capacity theory to multitrophic systems, we derive analytical insights into the persistence of trophically interacting species in heterogeneous landscapes and their responses to habitat changes. Our findings have implications for both ecological studies and conservation practices. In particular, our extended Levins’ rule offers a practical tool for managing habitat configuration to maintain food chain length, an important indicator of vertical diversity and ecosystem functioning (41). Our model framework may serve as a benchmark for future theoretical research that incorporates other scenarios of area and context dependency of colonization and extinction processes (30, 42, 43) and more complex trophic interactions to understand how metapopulation capacity may predict the spatial scaling of trophic structure (15, 44). While we have assumed constant dynamical parameters following habitat changes, future models should take into account potential habitat change–induced alterations in these parameters (e.g., the strength of trophic interactions [f]) (45). Such theoretical developments, together with empirical work designed to test theoretical predictions, will contribute to a more operational and predictive science of conservation biology.


Analytical Investigation of the Models.

Our analyses of the spatially explicit models of prey–predator interactions (Eq. 1) and food chains (Eq. 4) were built upon the classic results by ref. 24. Although deriving the exact persistence conditions for higher trophic levels was difficult (except under donor-control cases; see SI Appendix, Appendix S1), we could provide approximate solutions under a homogeneity assumption, in which species at the lower trophic level had the same occupancy in all patches. For brevity, all details on the derivation of our analytic solutions (summarized in Table 1) are presented in SI Appendix, Appendix S1.


We conducted numerical simulations to evaluate the accuracy of our approximate solutions for persistence conditions. Our simulated landscapes consisted of n patches (n = 20 to 100), with their areas (Ak) sampled from a log-uniform distribution (i.e., log2AkU[2,2]) and locations sampled from a two-dimensional uniform distribution on a square landscape of [0, 10] × [0, 10]. These simulated landscapes provided habitats for prey–predator pairs or food chains, for which we assigned dynamical parameters (ci, ei, αi, and fi) for each species. We simulated the dynamics of these trophic systems using Eq. 1 (for prey–predator pairs) or Eq. 4 (for food chains) until they reached equilibrium (10,000 time steps) and record status of each species (extinction or persistence). For prey–predator pairs, we gradually changed species (e.g., the extinction parameter in Fig. 1 and SI Appendix, Fig. S1) or landscape parameters (e.g., patch area in SI Appendix, Fig. S2) and determined numerically the threshold conditions for the persistence of each species. We then compared these numerical conditions to our analytical predictions (i.e., Eqs. 2 and 3). For food chains, each simulation was initialized with eight trophic levels, and the realized food chain length was recorded at the end of simulations. We compared the realized food chain length with the prediction by Eqs. 57 (Fig. 2 and SI Appendix, Fig. S5). In all simulations, we calculated the weighted average occupancy for species i by pi¯=kpikAk2kAk2 (24) and then derived the fraction of empty patches of the basal species as U1=1p1¯.
We simulated three scenarios of habitat changes and examined their effects on the persistence of food chains (i.e., Fig. 3). Our simulations started with a landscape consisting of 50 patches and hosting a food chain of eight trophic levels all with equal parameters (c, e, f, α). Three ways of habitat changes were then simulated: habitat deterioration (i.e., decreasing the area of all patches proportionally by 1-h), habitat loss (i.e., randomly removing a fraction of patches by 1-a), or habitat fragmentation (i.e., increasing the distance between each pair of patches by s). For each scenario and level (i.e., 1-h, 1-a, and s) of habitat changes, we ran the dynamics of food chains (Eq. 4) to equilibrium and recorded the realized food chain length. In simulating habitat loss, we randomly removed patches from 50 to 5 and repeated the process 100 times, from which we took the median food chain length at each level of habitat loss. We compared the simulation results with our analytic predictions by Eqs. 8 and 9.
See SI Appendix, Table S1 for model parameters used for generating each figure.

Empirical Data.

The butterfly M. cinxia in the Åland archipelago is a model system for the study of metapopulation dynamics (24). It inhabits a 50 × 70 km2 landscape made up of ∼4500 small meadow habitat patches, each surrounded by unsuitable habitat. The patches are surveyed annually to count the number of butterfly larval nests (since 1993) as well as the cocoons of the specialist parasitoid wasp, C. melitaearum (since 1998) (39, 46). The butterfly occupies several hundred meadows annually, with a turnover of between 50 and 200 local extinctions and colonizations (46). C. melitaearum occupies 5 to 15% of the local butterfly populations and also has a high rate of local population turnover (31, 39). Systematic sampling has also confirmed the presence of another parasitoid wasp, H. horticola, in virtually all M. cinxia populations (35).
The 4,500 habitat patches occur in clusters in the landscape. We have divided them into 131 networks using the software SPOMSIM (47) based on the connectivity of patches as a function of the distances between them, their areas, and the characteristic dispersal distance of the butterfly. Thus, butterflies inhabiting a network easily colonize patches within it, but movement among networks is rare. For this study, we used 10 y of survey data, from 2005 to 2014. We omitted those networks located on isolated islands where the wasps rarely appear (17/131) and containing patches with missing records during the survey period (16/131). This left 102 networks (note that four networks are on islands and contain missing values). For each network, we first calculated the temporal average of butterfly occupancy for each patch k (pk) and then derived the average fraction of empty patches by U1=1kpkωk, where the weight ωk was given by the squared area of the patch k (following ref. 24). We similarly calculated the temporal average of patch occupancy for the wasp and derived its average network-level occupancy (again weighted by ωk).

Data Availability

R code and Excel data have been deposited in GitHub (


We thank Jordi Bascompte, Andrew Gonzalez, and Mingyu Luo for helpful discussions. This work was supported by the National Natural Science Foundation of China (Grants 31988102 and 31870505). U.B. was supported by the German Research Foundation (DFG) through the Research Unit FOR-2716 (DynaCom) and German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig (Grant DFG: FZT 118, 202548816). M.L. was supported by the TULIP Laboratory of Excellence (Grant ANR-10-LABX-41).

Supporting Information

Appendix (PDF)


D. Tilman, P. Kareiva, Eds., Spatial Ecology: The Role of Space in Population Dynamics and Interspecific Interactions (Princeton University Press, 1997).
M. Holyoak, M. A. Leibold, R. D. Holt, Eds., Metacommunities: Spatial Dynamics and Ecological Communities (University of Chicago Press, 2005).
E. O. Wilson, R. H. MacArthur, The Theory of Island Biogeography (Princeton University Press, 1967).
R. Levins, Some demographic and genetic consequences of environmental heterogeneity for biological control. Am. Entomol. (Lanham Md.) 15, 237–240 (1969).
J. B. Losos, R. E. Ricklefs, The Theory of Island Biogeography Revisited (Princeton University Press, 2009).
L. Fahrig, Habitat fragmentation: A long and tangled tale. Glob. Ecol. Biogeogr. 28, 33–41 (2019).
N. M. Haddad et al., Habitat fragmentation and its lasting impact on Earth’s ecosystems. Sci. Adv. 1, e1500052 (2015).
I. Hanski, Metapopulation Ecology (Oxford University Press, 1999).
A. Kruess, T. Tscharntke, Habitat fragmentation, species loss, and biological control. Science 264, 1581–1584 (1994).
K. F. Davies, C. R. Margules, J. F. Lawrence, Which traits of species predict population declines in experimental forest fragments? Ecology 81, 1450–1461 (2000).
J. H. Lawton, S. Nee, A. J. Letcher, P. H. Harvey, “Animal distributions: Patterns and processes” in Large-Scale Ecology and Conservation Biology, P. J. Edwards, R. May, N. R. Webb, Eds. (Blackwell Scientific Publications, 1994), pp. 41–58.
I. Hanski, M. E. Gilpin, Metapopulation Biology: Ecology, Genetics, and Evolution (Academic Press, London, 1997).
R. Levins, D. Culver, Regional coexistence of species and competition between rare species. Proc. Natl. Acad. Sci. U.S.A. 68, 1246–1248 (1971).
R. D. Holt, “From metapopulation dynamics to community structure: Some consequences of spatial heterogeneity” in Metapopulation Biology (Academic Press, 1997), pp. 149–164.
R. D. Holt, J. H. Lawton, G. A. Polis, N. D. Martinez, Trophic rank and the species–area relationship. Ecology 80, 1495–1504 (1999).
S. Nee, R. M. May, Dynamics of metapopulations: Habitat destruction and competitive coexistence. J. Anim. Ecol. 61, 37–40 (1992).
A. Hastings, Disturbance, coexistence, history, and competition for space. Theor. Popul. Biol. 18, 363–373 (1980).
D. Tilman, C. L. Lehman, Habitat destruction and the extinction debt. Nature 371, 65–66 (1994).
D. Gravel, F. Massol, E. Canard, D. Mouillot, N. Mouquet, Trophic theory of island biogeography. Ecol. Lett. 14, 1010–1016 (2011).
P. Pillai, A. Gonzalez, M. Loreau, Metacommunity theory explains the emergence of food web complexity. Proc. Natl. Acad. Sci. U.S.A. 108, 19293–19298 (2011).
M. A. Leibold et al., The metacommunity concept: A framework for multi‐scale community ecology. Ecol. Lett. 7, 601–613 (2004).
M. A. Leibold, J. M. Chase, Metacommunity Ecology (Princeton University Press, 2018).
T. Gross et al., Modern models of trophic meta-communities. Philos. Trans. R. Soc. Lond. B Biol. Sci. 375, 20190455 (2020).
I. Hanski, O. Ovaskainen, The metapopulation capacity of a fragmented landscape. Nature 404, 755–758 (2000).
J. K. Schnell, G. M. Harris, S. L. Pimm, G. J. Russell, Estimating extinction risk with metapopulation models of large-scale fragmentation. Conserv. Biol. 27, 520–530 (2013).
M. Strimas-Mackey, J. F. Brodie, Reserve design to optimize the long-term persistence of multiple species. Ecol. Appl. 28, 1354–1361 (2018).
R. Huang, S. L. Pimm, C. Giri, Using metapopulation theory for practical conservation of mangrove endemic birds. Conserv. Biol. 34, 266–275 (2020).
L. J. Gilarranz, J. Bascompte, Spatial network structure and metapopulation persistence. J. Theor. Biol. 297, 11–16 (2012).
J. Grilli, G. Barabás, S. Allesina, Metapopulation persistence in random fragmented landscapes. PLOS Comput. Biol. 11, e1004251 (2015).
S. Wang, F. Altermatt, Metapopulations revisited: The area-dependence of dispersal matters. Ecology 100, e02792 (2019).
S. Van Nouhuys, I. Hanski, Colonization rates and distances of a host butterfly and two specific parasitoids in a fragmented landscape. J. Anim. Ecol. 71, 639–650 (2002).
I. Hanski, Spatial patterns of coexistence of competing species in patchy habitat. Theor. Ecol. 1, 29–43 (2008).
H. M. Martinson, W. F. Fagan, Trophic disruption: A meta-analysis of how habitat fragmentation affects resource consumption in terrestrial arthropod systems. Ecol. Lett. 17, 1178–1189 (2014).
J. M. Tylianakis, T. Tscharntke, O. T. Lewis, Habitat modification alters the structure of tropical host-parasitoid food webs. Nature 445, 202–205 (2007).
C. Couchoux, P. Seppä, S. van Nouhuys, Strong dispersal in a parasitoid wasp overwhelms habitat fragmentation and host population dynamics. Mol. Ecol. 25, 3344–3355 (2016).
K. S. McCann, Food Webs (Princeton University Press, Princeton, NJ, 2012).
I. Hanski, G. A. Zurita, M. I. Bellocq, J. Rybicki, Species-fragmented area relationship. Proc. Natl. Acad. Sci. U.S.A. 110, 12715–12720 (2013).
S. Nee, How populations persist. Nature 367, 123–124 (1994).
Ø. H. Opedal, O. Ovaskainen, M. Saastamoinen, A. L. Laine, S. van Nouhuys, Host-plant availability drives the spatiotemporal dynamics of interacting metapopulations across a fragmented landscape. Ecology 101, e03186 (2020).
R. Ryser et al., The biggest losers: Habitat isolation deconstructs complex food webs from top to bottom. Proc. Biol. Sci. 286, 20191177 (2019).
S. Wang, U. Brose, Biodiversity and ecosystem functioning in food webs: The vertical diversity hypothesis. Ecol. Lett. 21, 9–20 (2018).
O. Ovaskainen, Long-term persistence of species and the SLOSS problem. J. Theor. Biol. 218, 419–433 (2002).
E. A. Fronhofer et al., Bottom-up and top-down control of dispersal across major organismal groups. Nat. Ecol. Evol. 2, 1859–1863 (2018).
U. Brose, A. Ostling, K. Harrison, N. D. Martinez, Unified spatial scaling of species and their trophic interactions. Nature 428, 167–171 (2004).
J. van de Koppel et al., The effects of spatial scale on trophic interactions. Ecosystems (N. Y.) 8, 801–807 (2005).
S. P. Ojanen, M. Nieminen, E. Meyke, J. Pöyry, I. Hanski, Long-term metapopulation study of the Glanville fritillary butterfly (Melitaea cinxia): Survey methods, data management, and long-term population trends. Ecol. Evol. 3, 3713–3737 (2013).
A. Moilanen, SPOMSIM: Software for stochastic patch occupancy models of metapopulation dynamics. Ecol. Modell. 179, 533–550 (2004).

Information & Authors


Published in

Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 118 | No. 34
August 24, 2021
PubMed: 34417316


Data Availability

R code and Excel data have been deposited in GitHub (

Submission history

Published online: August 20, 2021
Published in issue: August 24, 2021


  1. fragmentation
  2. habitat changes
  3. heterogeneous landscapes
  4. trophic interactions


We thank Jordi Bascompte, Andrew Gonzalez, and Mingyu Luo for helpful discussions. This work was supported by the National Natural Science Foundation of China (Grants 31988102 and 31870505). U.B. was supported by the German Research Foundation (DFG) through the Research Unit FOR-2716 (DynaCom) and German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig (Grant DFG: FZT 118, 202548816). M.L. was supported by the TULIP Laboratory of Excellence (Grant ANR-10-LABX-41).


This article is a PNAS Direct Submission.



Shaopeng Wang1 [email protected]
Institute of Ecology, Key Laboratory for Earth Surface Processes of the Ministry of Education, College of Urban and Environmental Sciences, Peking University, 100871 Beijing, China;
Ulrich Brose
EcoNetLab, German Centre for Integrative Biodiversity Research Halle-Jena-Leipzig, 04103 Leipzig, Germany;
Institute of Biodiversity, Friedrich Schiller University Jena, 07743 Jena, Germany;
Organismal and Evolutionary Biology Research Program, University of Helsinki, FI-00014 Helsinki, Finland;
Department of Ecology and Evolutionary Biology, Cornell University, Ithaca, NY 14853;
Department of Biology, University of Florida, Gainesville, FL 32611;
Michel Loreau
Theoretical and Experimental Ecology Station, CNRS, Moulis 09200, France


To whom correspondence may be addressed. Email: [email protected].
Author contributions: S.W. designed research; S.W. performed research; S.W. contributed new reagents/analytic tools; S.W. and S.v.N. analyzed data; S.W. wrote the paper; and S.W., U.B., S.v.N., R.D.H., and M.L. interpreted the results and revised the manuscript.

Competing Interests

The authors declare no competing interest.

Metrics & Citations


Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.

Citation statements



If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by


    View Options

    View options

    PDF format

    Download this article as a PDF file


    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to access the full text.

    Single Article Purchase

    Metapopulation capacity determines food chain length in fragmented landscapes
    Proceedings of the National Academy of Sciences
    • Vol. 118
    • No. 34







    Share article link

    Share on social media