Levins Model

The Levins Model is a foundational continuous-time ordinary differential equation in spatial and population ecology that describes the regional persistence of a metapopulation—a network of spatially discrete local populations connected by dispersal across an uninhabitable landscape matrix. Formulated by American mathematical ecologist Richard Levins in 1969 and 1970 to design biological control strategies against agricultural insect pests, the model shifts the level of demographic tracking from individual organisms (N) to the fraction of occupied habitat patches (p) across an idealized network of infinitely many identical patches.

By treating each local subpopulation as a binary unit (either occupied or unoccupied) subject to stochastic local extinction and recolonization via propagule dispersal, the Levins equation demonstrates that a species can persist indefinitely at the regional scale even when every individual local population is unstable and guaranteed to go extinct over finite timescales. The Levins framework forms the mathematical basis for metapopulation conservation biology, the extinction debt induced by habitat destruction, epidemiological susceptible-infected-susceptible (SIS) patch dynamics, and competitive coexistence via competition-colonization trade-offs.

The Canonical Patch-Occupancy Equation

Consider a landscape containing a large number of discrete, homogeneous habitat patches, of which a proportion p(t) ∈ [0, 1] is occupied by the focal species at time t, while the complementary proportion 1 − p(t) consists of vacant, colonizable patches. Assuming mean-field (global) dispersal—meaning propagules dispersing from any occupied patch have an equal probability of reaching any other patch in the network—and assuming that local population growth within a newly colonized patch reaches local carrying capacity instantaneously relative to the timescale of inter-patch turnover, Levins (1969) wrote the rate of change in patch occupancy as:

dp / dt = c · p ( 1 − p ) − e · p

Here:

• p(t) (dimensionless, 0 ≤ p ≤ 1) is the fraction of suitable habitat patches occupied at time t.
• c (units: time⁻¹) is the per-patch colonization rate parameter, representing the rate at which propagules produced by an occupied patch disperse across the landscape and successfully establish new subpopulations if all encountered patches were empty. Because only the fraction (1 − p) of patches is unoccupied, the total realized regional colonization flux is c · p (1 − p).
• e (units: time⁻¹) is the local patch extinction rate parameter, representing the density-independent per-patch probability per unit time that an occupied patch suffers stochastic demographic or environmental collapse and reverts to an empty state. The reciprocal 1 / e equals the mean lifetime of an individual local population.

Equilibrium Occupancy and the Persistence Threshold

Algebraically rearranging the Levins differential equation into Verhulst logistic form highlights its structural isomorphism with single-species logistic growth:

dp / dt = ( c − e ) p [ 1 − p / ( 1 − e / c ) ]

Setting dp / dt = 0 yields two equilibrium solutions: a trivial extinction equilibrium p₀* = 0 and a non-trivial regional metapopulation equilibrium (p*):

p* = 1 − e / c

Evaluating the derivative f′(p) = c − e − 2*cp* at the two equilibria establishes the Levins Metapopulation Persistence Threshold:

1. Subcritical Extinction Regime (c ≤ e, or c / e ≤ 1): When the colonization rate does not exceed the local extinction rate, the eigenvalue at the origin is non-positive (f′(0) = c − e ≤ 0). The extinction equilibrium p₀* = 0 is globally asymptotically stable, and the metapopulation decays to regional extinction regardless of initial occupancy p(0).
2. Supercritical Persistence Regime (c > e, or R_M = c / e > 1): When the dimensionless basic metapopulation reproduction number (R_M = c / e)—the average number of secondary patches colonized by a single occupied patch during its expected lifetime 1 / e in an otherwise empty landscape—exceeds unity, the origin becomes unstable (f′(0) = c − e > 0) and the interior equilibrium p* = 1 − 1 / R_M is globally asymptotically stable (f′(p*) = −(c − e) < 0). Notably, even in a thriving metapopulation (c > e), a finite fraction of perfectly suitable habitat patches (e / c = 1 − p*) remains unoccupied at any given instant.

Hanski's Rescue Effect and Mainland-Island Propagule Rain

Subsequent work by James H. Brown and Astrid Kodric-Brown (1977), Thomas W. Schoener (1979), and Ilkka Hanski (1982, 1999) relaxed two simplifying assumptions of the classical Levins equation:

• The Internal Rescue Effect (Hanski, 1982): In the classical Levins model, immigration into already-occupied patches is ignored. In reality, when regional occupancy p is high, a continuous influx of immigrants supplements declining local populations and prevents stochastic extinction. Modeling the effective extinction rate as a linearly decreasing function of occupancy, E(p) = e (1 − p) · p, transforms the dynamics into dp / dt = cp(1 − p) − ep(1 − p) = (c − e)p(1 − p). Under a partial rescue effect E(p) = e(1 − a p)p (with 0 < a < 1), the model can generate bistable core-satellite dynamics and metapopulation Allee effects.
• The Mainland-Island (Propagule Rain) Model (Gotelli, 1991): When colonization occurs via a constant external influx of propagules (I) from a permanent mainland source rather than internal patch-to-patch dispersal, the colonization term becomes I(1 − p) and the equation is linear: dp / dt = I(1 − p) − e p, yielding a globally stable equilibrium p* = I / (I + e) that never goes extinct for any I > 0.

Lande-Nee-May Habitat Destruction Model and Extinction Debt

In 1987, Russell Lande (and subsequently Sean Nee and Robert M. May, 1992) modified the Levins equation to evaluate the impact of permanent habitat destruction on fragmented species, providing a mathematical foundation for umbrella species reserve design. Suppose a fraction D (0 ≤ D < 1) of all habitat patches in a landscape is permanently destroyed, leaving only the fraction h = 1 − D of patches habitable. Because dispersing propagules land randomly across the landscape, only the fraction (1 − D − p) of the total landscape represents vacant, habitable patches available for colonization:

dp / dt = c · p ( 1 − D − p ) − e · p

Setting dp / dt = 0 for p > 0 gives the Lande-Nee-May equilibrium occupancy:

p* = ( 1 − D ) − e / c     (as a fraction of the total original landscape)

or, expressed as the fraction of remaining habitable patches occupied:

p_hab* = p* / ( 1 − D ) = 1 − ( e / c ) / ( 1 − D )

This equation reveals a critical conservation threshold: the metapopulation collapses to regional extinction (p* → 0) before all habitat is destroyed, specifically when the fraction of destroyed habitat reaches the critical destruction threshold (D_c):

D_c = 1 − e / c = p₀*

Remarkably, the maximum fraction of habitat that can be destroyed without causing ultimate regional extinction (D_c) is identical to the fraction of patches (p₀* = 1 − e / c) that the species occupied in the pristine, undisturbed landscape (D = 0). When habitat destruction crosses D > D_c, the species does not vanish overnight; instead, occupancy decays exponentially toward zero over multiple generations—a delayed deterministic trajectory termed extinction debt by David Tilman, Robert May, Clarence Lehman, and Martin Nowak (1994).

Tilman's Hierarchical Competition-Colonization Coexistence

By coupling the Levins patch-occupancy framework with the Competitive Lotka-Volterra Equations, J. G. Skellam (1951), Richard Levins and David Culver (1971), Alan Hastings (1980), and David Tilman (1994) demonstrated how spatial patch turnover allows an unlimited number of species to coexist on a single limiting resource within the biosphere without violating the local competitive exclusion principle. Suppose Species 1 is a superior local competitor (a late-successional K-selected species) that instantaneously displaces Species 2 whenever it colonizes a patch occupied by Species 2, whereas Species 2 is an inferior competitor but a superior colonizer (r-selected fugitive species with c₂ > c₁):

dp₁ / dt = c₁ p₁ ( 1 − p₁ ) − e₁ p₁

dp₂ / dt = c₂ p₂ ( 1 − p₁ − p₂ ) − e₂ p₂ − c₁ p₁ p₂

At equilibrium, the superior competitor occupies p₁* = 1 − e₁ / c₁, and the fugitive species stably coexists in the remaining vacant patches provided its colonization rate exceeds the critical threshold c₂ > c₁ ( c₁ + e₂ − e₁ ) / e₁. Paradoxically, under habitat destruction D, the superior competitor (which typically has the lowest colonization rate c₁) suffers the earliest extinction debt.

Incidence Function Model (IFM) and Real Landscapes

Because real landscapes violate Levins's assumption of identical patch areas and global dispersal—as shown in Ilkka Hanski's long-term empirical study of the Glanville fritillary butterfly (Melitaea cinxia) across ~4,000 dry meadow patches (Plantago lanceolata and Veronica spicata host plants) in the Åland Islands of Finland—Hanski (1994) formulated the Spatially Realistic Incidence Function Model (IFM). In the IFM, patch i has area A_i (scaling local extinction as E_i ∝ A_i^−x) and inter-patch Euclidean distance d_ij (scaling connectivity as S_i = ∑_j≠i p_j exp(−α d_ij) A_j^b). Replacing the homogeneous patch count with the leading eigenvalue (metapopulation capacity λ_M) of the landscape connectivity matrix M restores the exact mathematical form of the Levins persistence condition (λ_M > e / c), bridging idealized analytical theory and real-world reserve design.

1. Levins, R. (1969). Some demographic and genetic consequences of environmental heterogeneity for biological control. Bulletin of the Entomological Society of America, 15(3), 237–240. https://doi.org/10.1093/besa/15.3.237
2. Lande, R. (1987). Extinction thresholds in demographic models of territorial populations. The American Naturalist, 130(4), 624–635. https://doi.org/10.1086/284734
3. Tilman, D., May, R. M., Lehman, C. L., & Nowak, M. A. (1994). Habitat destruction and the extinction debt. Nature, 371(6492), 65–66. https://doi.org/10.1038/371065a0
4. Hanski, I. (1998). Metapopulation dynamics. Nature, 396(6706), 41–49. https://doi.org/10.1038/23876