Checking a metapopulation model

metapopulation

spatial ecology

model checking

ecology tutorial

Three diagnostics for a spatial metapopulation model in R: the equilibrium assumption, a weakly identified dispersal kernel, and stochastic extinction risk.

Author

Tidy Ecology

Published

2026-07-21

The three previous tutorials in this set are convenient to the point of being seductive. The incidence function model fits from a single snapshot; the metapopulation capacity boils a whole landscape down to one eigenvalue; the persistence rule is a clean inequality. Each convenience rests on an assumption that a snapshot cannot check on its own, and when the assumption fails the model does not complain: it returns a confident, wrong answer. This tutorial builds three diagnostics in base R for exactly those failure modes. It reuses the clustered landscape and the fitting machinery from the incidence function post, so the checks apply directly to what came before.

Efun <- function(A, mu, x) pmin(1, mu / A^x)
conn <- function(o, D, A, alpha) {
  M <- exp(-alpha * D) * outer(rep(1, length(A)), A)
  diag(M) <- 0
  as.numeric(M %*% o)
}
Cfun <- function(S, y) S^2 / (S^2 + y^2)

set.seed(2024)
K <- 6; per <- 26
centres <- cbind(runif(K, 3, 25), runif(K, 3, 25))
clus <- do.call(rbind, lapply(1:K, function(k)
  cbind(rnorm(per, centres[k, 1], 1.3), rnorm(per, centres[k, 2], 1.3))))
iso <- cbind(runif(44, 1, 27), runif(44, 1, 27))
patch <- rbind(clus, iso)
n <- nrow(patch)
A <- rlnorm(n, meanlog = 0, sdlog = 0.9)
D <- as.matrix(dist(patch))
mu_t <- 0.8; x_t <- 1.0; y_t <- 2.5; alpha_t <- 0.6
Etrue <- Efun(A, mu_t, x_t)

fitIFM <- function(o, af = 0.6) {
  S <- conn(o, D, A, af)
  nll <- function(p) {
    mu <- exp(p[1]); x <- exp(p[2]); y <- exp(p[3])
    C <- Cfun(S, y); E <- Efun(A, mu, x)
    J <- pmin(pmax(C / (C + E), 1e-9), 1 - 1e-9)
    -sum(dbinom(o, 1, J, log = TRUE))
  }
  f <- optim(log(c(0.5, 1, 2)), nll, method = "Nelder-Mead",
             control = list(maxit = 4000, reltol = 1e-12))
  f <- optim(f$par, nll, method = "BFGS")
  list(par = exp(f$par), ll = -f$value)
}
n

[1] 200

Check one: is the snapshot at equilibrium?

The incidence function model reads the balance of colonisation and extinction off the spatial pattern, on the assumption that the pattern is stationary. A recovering or collapsing network breaks that. To see the damage, fit two snapshots of the same landscape: one at stochastic equilibrium, and one taken while the network is still climbing back from a crash.

The equilibrium snapshot is the one from the incidence function tutorial, run to stationarity. The transient snapshot starts from five per cent occupied and runs only six gentle steps, so the network is a long way from settled. The gentle step (each rate scaled by a small factor) leaves the equilibrium unchanged but slows the approach, making the transient unambiguous.

o <- as.integer(runif(n) < 0.5)
set.seed(303)
for (t in 1:2000) {
  S <- conn(o, D, A, alpha_t); C <- Cfun(S, y_t)
  new <- o; emp <- o == 0L; occ <- o == 1L
  new[emp] <- as.integer(runif(sum(emp)) < C[emp])
  new[occ] <- as.integer(runif(sum(occ)) >= Etrue[occ])
  o <- new
}
o_eq <- o

tau <- 0.25
stepT <- function(o) {
  S <- conn(o, D, A, alpha_t); C <- Cfun(S, y_t)
  new <- o; emp <- o == 0L; occ <- o == 1L
  new[emp] <- as.integer(runif(sum(emp)) < C[emp] * tau)
  new[occ] <- as.integer(runif(sum(occ)) >= Etrue[occ] * tau)
  new
}
set.seed(505)
o <- as.integer(runif(n) < 0.05)
for (t in 1:6) o <- stepT(o)
o_tr <- o
c(equilibrium_occupied = sum(o_eq), transient_occupied = sum(o_tr))

equilibrium_occupied   transient_occupied 
                 117                   76

The equilibrium snapshot has 117 occupied patches; the transient, still recovering, has 76. Now fit the incidence function model to each and compare with the truth (mu = 0.8, x = 1.0, y = 2.5).

fe <- fitIFM(o_eq); ft <- fitIFM(o_tr)
rbind(equilibrium = round(fe$par, 3),
      transient   = round(ft$par, 3))

             [,1]  [,2]  [,3]
equilibrium 0.441 0.841 3.620
transient   2.344 0.181 2.753

The equilibrium fit lands where the incidence function tutorial left it: x = 0.841, close enough to the true 1.0. The transient fit is wrecked: the area exponent x collapses to 0.181, which says extinction barely depends on area at all. A conservation planner reading that number would conclude, wrongly, that patch size hardly matters. The mechanism is worth seeing directly. Correlate occupancy with area and with connectivity in each snapshot.

S_eq <- conn(o_eq, D, A, alpha_t); S_tr <- conn(o_tr, D, A, alpha_t)
rbind(equilibrium = c(area = cor(o_eq, log(A)), connectivity = cor(o_eq, log(S_eq + 1e-6))),
      transient   = c(area = cor(o_tr, log(A)), connectivity = cor(o_tr, log(S_tr + 1e-6)))) |>
  round(3)

              area connectivity
equilibrium  0.356        0.298
transient   -0.056        0.376

At equilibrium, occupancy tracks area (correlation 0.36): large patches persist, small ones wink out, which is the extinction-area signal the model is built to read. At the transient, that correlation vanishes (-0.056) while the connectivity correlation stays positive (0.376). Recovery is driven by connectivity: the well-connected patches near the surviving core refill first, regardless of their size, so the snapshot carries no area signal for the model to find. That is why x collapses. The tell is simple: a near-zero area exponent from a system you have reason to think is disturbed should be read as a warning about equilibrium, not as biology.

Two panels of occupied and empty points spread across log patch area. Left (equilibrium) shows a logistic trend rising left to right; right (transient) shows an almost horizontal trend. — Figure 1: Occupancy against patch area at equilibrium and during recovery from a crash. The trend line rises with area at equilibrium (extinction hits small patches) but is flat during recovery, when connectivity rather than area decides who is occupied. The flat transient trend is what drives the area exponent to near zero.

Check two: can the snapshot see the dispersal kernel?

The incidence function tutorial fixed alpha by hand and promised an explanation. Here it is: a single snapshot barely constrains the dispersal parameter, yet the metapopulation capacity of the previous tutorial depends on it heavily. Profile the likelihood of the equilibrium snapshot over a range of alpha, re-fitting the other three parameters at each value.

proflik <- function(al) {
  Sa <- conn(o_eq, D, A, al)
  nl <- function(p) {
    mu <- exp(p[1]); x <- exp(p[2]); y <- exp(p[3])
    Cc <- Cfun(Sa, y); Ee <- Efun(A, mu, x)
    Jj <- pmin(pmax(Cc / (Cc + Ee), 1e-9), 1 - 1e-9)
    -sum(dbinom(o_eq, 1, Jj, log = TRUE))
  }
  starts <- list(log(c(0.5, 1, 2)), log(c(1, 0.5, 1)),
                 log(c(0.3, 1.5, 4)), log(c(2, 1, 0.5)))
  best <- Inf
  for (s in starts) {
    r <- optim(s, nl, method = "Nelder-Mead", control = list(maxit = 3000, reltol = 1e-11))
    r <- optim(r$par, nl, method = "BFGS")
    if (r$value < best) best <- r$value
  }
  -best
}
alg <- seq(0.2, 1.4, by = 0.1)
pl <- sapply(alg, proflik)
inb <- range(alg[pl >= max(pl) - 1.92])
c(logLik_span = round(max(pl) - min(pl), 2), ci_low = inb[1], ci_high = inb[2])

logLik_span      ci_low     ci_high 
       3.99        0.20        0.60

Over a seven-fold range of alpha, the whole profile log-likelihood moves by only 3.99 units. Everything from alpha = 0.2 to 0.6 sits within the usual 1.92-unit interval, so the snapshot cannot tell a short kernel from one three times as wide. The true value, 0.6, is comfortably inside that band, but so is almost everything else.

Figure 2: Profile log-likelihood for the dispersal parameter from a single snapshot. The curve is nearly flat: values across a wide range fall within the dashed interval line, so the snapshot leaves the kernel almost unidentified.

That flatness matters because the metapopulation capacity is not flat in alpha at all. Feeding the two ends of the interval into the capacity calculation shows how much rides on a parameter the data cannot pin down.

Mmat <- function(al) { M <- exp(-al * D) * outer(A, A); diag(M) <- 0; M }
lamOf <- function(al) eigen(Mmat(al), symmetric = TRUE, only.values = TRUE)$values[1]
c(lambda_at_low = round(lamOf(inb[1]), 1),
  lambda_at_high = round(lamOf(inb[2]), 1),
  fold_range = round(lamOf(inb[1]) / lamOf(inb[2]), 1))

 lambda_at_low lambda_at_high     fold_range 
         208.5           86.9            2.4

The capacity ranges from 209 down to 87 across the plausible alpha, a 2.4-fold spread. Every persistence verdict and every patch-value ranking from the capacity tutorial inherits that uncertainty. The honest way to report a capacity from a snapshot is as a range this wide, not a single number, unless alpha has been fixed from independent dispersal data.

Check three: is the threshold a promise?

The capacity rule says persist when lambda_M exceeds the species threshold delta. That is a deterministic statement about an infinite system. A real metapopulation has finitely many patches and turns over by chance, so it can wink out even when the deterministic rule says it is safe. To measure the gap, simulate a spatially explicit stochastic model on the same landscape: each occupied patch goes extinct at a rate that falls with area, each empty patch is colonised at a rate that rises with connectivity, and we ask how often the whole network dies within a fixed horizon.

lamM <- lamOf(alpha_t)
sim_srl <- function(cc, ee, steps, tau, seed) {
  set.seed(seed)
  o <- rep(1L, n)
  for (t in 1:steps) {
    S <- conn(o, D, A, alpha_t)
    pc <- 1 - exp(-cc * S * tau); pe <- 1 - exp(-(ee / A) * tau)
    new <- o; emp <- o == 0L; occ <- o == 1L
    new[emp] <- as.integer(runif(sum(emp)) < pc[emp])
    new[occ] <- as.integer(runif(sum(occ)) >= pe[occ])
    o <- new
    if (sum(o) == 0) return(TRUE)
  }
  FALSE
}
cc <- 1; R <- 250; H <- 100; tau_s <- 0.4
ratios <- c(0.3, 0.5, 0.7, 0.9, 1.1, 1.3)
pext <- sapply(ratios, function(rr) {
  ee <- rr * lamM * cc
  mean(sapply(1:R, function(k) sim_srl(cc, ee, H, tau_s, seed = 5000 + round(rr * 1000) + k)))
})
data.frame(delta_over_lambda = ratios, extinction_prob = round(pext, 3))

  delta_over_lambda extinction_prob
1               0.3           0.044
2               0.5           0.484
3               0.7           0.804
4               0.9           0.932
5               1.1           0.972
6               1.3           0.992

The deterministic rule is a step: safe below a ratio of one, doomed above it. The stochastic reality is a smooth ramp. At a ratio of 0.5 (that is, lambda_M twice the threshold, deep in the supposedly safe zone), the network still dies out 48% of the time within the horizon; at a ratio of 0.7 the extinction probability is already 80%. Comfort well below the threshold is not comfort at all.

A smooth S-shaped curve of extinction probability rising from near zero to near one as the ratio increases, overlaid on a step function that jumps from zero to one at a ratio of one. A vertical dashed line marks the threshold. — Figure 3: Stochastic extinction probability against the ratio of species threshold to metapopulation capacity, with the deterministic verdict as a step. The deterministic rule flips from safe to doomed at a ratio of one; the simulated extinction risk rises smoothly and is already substantial well inside the safe side.

The lesson generalises the Levins tutorial’s stochastic detour to a real landscape: a threshold is necessary, not sufficient. A finite metapopulation needs its capacity to sit well above the threshold, not just barely above it, before persistence is anything like assured.

What the checks can and cannot do

The three diagnostics constrain the model but none of them can confirm it, and the reason is the same in each case: a snapshot is a single frame of a moving system. The equilibrium check can flag a pattern that looks disturbed, but it cannot prove a quiet-looking snapshot really is at rest rather than drifting slowly. The kernel check honestly reports that one snapshot leaves alpha loose, but it cannot tighten it; only independent dispersal data can. And the stochastic check shows that persistence is a probability, not a fact, so no single observation can certify a metapopulation as safe. Persistence is a claim about the future, and the only real test is more snapshots through time: watch the same network across years, and the equilibrium assumption, the turnover rates, and the persistence verdict all become checkable in a way that one map never allows. These tutorials fit what a single snapshot can give; knowing its limits is what keeps the tidy numbers honest.

References

Hanski I 1994. Journal of Animal Ecology 63(1):151-162 (10.2307/5591).

Ovaskainen O, Hanski I 2002. Theoretical Population Biology 61(3):285-295 (10.1006/tpbi.2002.1586).

Moilanen A 2002. Oikos 96(3):516-530 (10.1034/j.1600-0706.2002.960313.x).

Lande R 1993. The American Naturalist 142(6):911-927 (10.1086/285580).

Hanski I, Ovaskainen O 2000. Nature 404(6779):755-758 (10.1038/35008063).

--- title: "Checking a metapopulation model" description: "Three diagnostics for a spatial metapopulation model in R: the equilibrium assumption, a weakly identified dispersal kernel, and stochastic extinction risk." date: "2026-07-21 12:00" categories: [R, metapopulation, spatial ecology, model checking, ecology tutorial] image: thumbnail.png image-alt: "Three small panels: a flat occupancy-versus-area trend, a nearly flat likelihood curve, and a smooth extinction-probability curve crossing a step function." --- The three previous tutorials in this set are convenient to the point of being seductive. The incidence function model fits from a single snapshot; the metapopulation capacity boils a whole landscape down to one eigenvalue; the persistence rule is a clean inequality. Each convenience rests on an assumption that a snapshot cannot check on its own, and when the assumption fails the model does not complain: it returns a confident, wrong answer. This tutorial builds three diagnostics in base R for exactly those failure modes. It reuses the clustered landscape and the fitting machinery from the incidence function post, so the checks apply directly to what came before. ```{r} #| label: setup Efun <- function(A, mu, x) pmin(1, mu / A^x) conn <- function(o, D, A, alpha) { M <- exp(-alpha * D) * outer(rep(1, length(A)), A) diag(M) <- 0 as.numeric(M %*% o) } Cfun <- function(S, y) S^2 / (S^2 + y^2) set.seed(2024) K <- 6; per <- 26 centres <- cbind(runif(K, 3, 25), runif(K, 3, 25)) clus <- do.call(rbind, lapply(1:K, function(k) cbind(rnorm(per, centres[k, 1], 1.3), rnorm(per, centres[k, 2], 1.3)))) iso <- cbind(runif(44, 1, 27), runif(44, 1, 27)) patch <- rbind(clus, iso) n <- nrow(patch) A <- rlnorm(n, meanlog = 0, sdlog = 0.9) D <- as.matrix(dist(patch)) mu_t <- 0.8; x_t <- 1.0; y_t <- 2.5; alpha_t <- 0.6 Etrue <- Efun(A, mu_t, x_t) fitIFM <- function(o, af = 0.6) { S <- conn(o, D, A, af) nll <- function(p) { mu <- exp(p[1]); x <- exp(p[2]); y <- exp(p[3]) C <- Cfun(S, y); E <- Efun(A, mu, x) J <- pmin(pmax(C / (C + E), 1e-9), 1 - 1e-9) -sum(dbinom(o, 1, J, log = TRUE)) } f <- optim(log(c(0.5, 1, 2)), nll, method = "Nelder-Mead", control = list(maxit = 4000, reltol = 1e-12)) f <- optim(f$par, nll, method = "BFGS") list(par = exp(f$par), ll = -f$value) } n ``` ## Check one: is the snapshot at equilibrium? The incidence function model reads the balance of colonisation and extinction off the spatial pattern, on the assumption that the pattern is stationary. A recovering or collapsing network breaks that. To see the damage, fit two snapshots of the same landscape: one at stochastic equilibrium, and one taken while the network is still climbing back from a crash. The equilibrium snapshot is the one from the incidence function tutorial, run to stationarity. The transient snapshot starts from five per cent occupied and runs only six gentle steps, so the network is a long way from settled. The gentle step (each rate scaled by a small factor) leaves the equilibrium unchanged but slows the approach, making the transient unambiguous. ```{r} #| label: snapshots o <- as.integer(runif(n) < 0.5) set.seed(303) for (t in 1:2000) { S <- conn(o, D, A, alpha_t); C <- Cfun(S, y_t) new <- o; emp <- o == 0L; occ <- o == 1L new[emp] <- as.integer(runif(sum(emp)) < C[emp]) new[occ] <- as.integer(runif(sum(occ)) >= Etrue[occ]) o <- new } o_eq <- o tau <- 0.25 stepT <- function(o) { S <- conn(o, D, A, alpha_t); C <- Cfun(S, y_t) new <- o; emp <- o == 0L; occ <- o == 1L new[emp] <- as.integer(runif(sum(emp)) < C[emp] * tau) new[occ] <- as.integer(runif(sum(occ)) >= Etrue[occ] * tau) new } set.seed(505) o <- as.integer(runif(n) < 0.05) for (t in 1:6) o <- stepT(o) o_tr <- o c(equilibrium_occupied = sum(o_eq), transient_occupied = sum(o_tr)) ``` The equilibrium snapshot has `r sum(o_eq)` occupied patches; the transient, still recovering, has `r sum(o_tr)`. Now fit the incidence function model to each and compare with the truth (`mu` = 0.8, `x` = 1.0, `y` = 2.5). ```{r} #| label: transient-fit fe <- fitIFM(o_eq); ft <- fitIFM(o_tr) rbind(equilibrium = round(fe$par, 3), transient = round(ft$par, 3)) ``` The equilibrium fit lands where the incidence function tutorial left it: `x` = `r round(fe$par[2],3)`, close enough to the true 1.0. The transient fit is wrecked: the area exponent `x` collapses to `r round(ft$par[2],3)`, which says extinction barely depends on area at all. A conservation planner reading that number would conclude, wrongly, that patch size hardly matters. The mechanism is worth seeing directly. Correlate occupancy with area and with connectivity in each snapshot. ```{r} #| label: mechanism S_eq <- conn(o_eq, D, A, alpha_t); S_tr <- conn(o_tr, D, A, alpha_t) rbind(equilibrium = c(area = cor(o_eq, log(A)), connectivity = cor(o_eq, log(S_eq + 1e-6))), transient = c(area = cor(o_tr, log(A)), connectivity = cor(o_tr, log(S_tr + 1e-6)))) |> round(3) ``` At equilibrium, occupancy tracks area (correlation `r round(cor(o_eq,log(A)),2)`): large patches persist, small ones wink out, which is the extinction-area signal the model is built to read. At the transient, that correlation vanishes (`r round(cor(o_tr,log(A)),3)`) while the connectivity correlation stays positive (`r round(cor(o_tr,log(S_tr+1e-6)),3)`). Recovery is driven by connectivity: the well-connected patches near the surviving core refill first, regardless of their size, so the snapshot carries no area signal for the model to find. That is why `x` collapses. The tell is simple: a near-zero area exponent from a system you have reason to think is disturbed should be read as a warning about equilibrium, not as biology. ```{r} #| label: fig-transient #| fig-cap: "Occupancy against patch area at equilibrium and during recovery from a crash. The trend line rises with area at equilibrium (extinction hits small patches) but is flat during recovery, when connectivity rather than area decides who is occupied. The flat transient trend is what drives the area exponent to near zero." #| fig-alt: "Two panels of occupied and empty points spread across log patch area. Left (equilibrium) shows a logistic trend rising left to right; right (transient) shows an almost horizontal trend." #| echo: false #| warning: false library(ggplot2) dft <- rbind( data.frame(logA = log(A), occ = o_eq, panel = "equilibrium"), data.frame(logA = log(A), occ = o_tr, panel = "transient (recovering)")) ggplot(dft, aes(logA, occ)) + geom_point(aes(colour = panel), position = position_jitter(height = 0.04), alpha = 0.5, size = 1.4, show.legend = FALSE) + geom_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE, colour = "#275139", linewidth = 1) + scale_colour_manual(values = c("equilibrium" = "#2f8f63", "transient (recovering)" = "#b5534e")) + facet_wrap(~panel) + labs(x = "log patch area", y = "occupied") + theme_minimal(base_size = 13) + theme(panel.grid.minor = element_blank()) ``` ## Check two: can the snapshot see the dispersal kernel? The incidence function tutorial fixed `alpha` by hand and promised an explanation. Here it is: a single snapshot barely constrains the dispersal parameter, yet the metapopulation capacity of the previous tutorial depends on it heavily. Profile the likelihood of the equilibrium snapshot over a range of `alpha`, re-fitting the other three parameters at each value. ```{r} #| label: alpha-profile proflik <- function(al) { Sa <- conn(o_eq, D, A, al) nl <- function(p) { mu <- exp(p[1]); x <- exp(p[2]); y <- exp(p[3]) Cc <- Cfun(Sa, y); Ee <- Efun(A, mu, x) Jj <- pmin(pmax(Cc / (Cc + Ee), 1e-9), 1 - 1e-9) -sum(dbinom(o_eq, 1, Jj, log = TRUE)) } starts <- list(log(c(0.5, 1, 2)), log(c(1, 0.5, 1)), log(c(0.3, 1.5, 4)), log(c(2, 1, 0.5))) best <- Inf for (s in starts) { r <- optim(s, nl, method = "Nelder-Mead", control = list(maxit = 3000, reltol = 1e-11)) r <- optim(r$par, nl, method = "BFGS") if (r$value < best) best <- r$value } -best } alg <- seq(0.2, 1.4, by = 0.1) pl <- sapply(alg, proflik) inb <- range(alg[pl >= max(pl) - 1.92]) c(logLik_span = round(max(pl) - min(pl), 2), ci_low = inb[1], ci_high = inb[2]) ``` Over a seven-fold range of `alpha`, the whole profile log-likelihood moves by only `r round(max(pl)-min(pl),2)` units. Everything from `alpha` = `r inb[1]` to `r inb[2]` sits within the usual 1.92-unit interval, so the snapshot cannot tell a short kernel from one three times as wide. The true value, 0.6, is comfortably inside that band, but so is almost everything else. ```{r} #| label: fig-alpha #| fig-cap: "Profile log-likelihood for the dispersal parameter from a single snapshot. The curve is nearly flat: values across a wide range fall within the dashed interval line, so the snapshot leaves the kernel almost unidentified." #| fig-alt: "A gently sloping, nearly flat curve of profile log-likelihood against the dispersal parameter, with a horizontal dashed line close below the peak intersecting a wide span of the curve." #| echo: false #| warning: false dfa <- data.frame(alpha = alg, ll = pl) ggplot(dfa, aes(alpha, ll)) + geom_hline(yintercept = max(pl) - 1.92, linetype = "dashed", colour = "#b5534e") + geom_line(colour = "#275139", linewidth = 1) + geom_point(colour = "#275139", size = 1.6) + annotate("text", x = 1.1, y = max(pl) - 1.92 + 0.25, label = "1.92 below peak", colour = "#b5534e", size = 3.3) + labs(x = "dispersal parameter alpha", y = "profile log-likelihood") + theme_minimal(base_size = 13) + theme(panel.grid.minor = element_blank()) ``` That flatness matters because the metapopulation capacity is not flat in `alpha` at all. Feeding the two ends of the interval into the capacity calculation shows how much rides on a parameter the data cannot pin down. ```{r} #| label: alpha-capacity Mmat <- function(al) { M <- exp(-al * D) * outer(A, A); diag(M) <- 0; M } lamOf <- function(al) eigen(Mmat(al), symmetric = TRUE, only.values = TRUE)$values[1] c(lambda_at_low = round(lamOf(inb[1]), 1), lambda_at_high = round(lamOf(inb[2]), 1), fold_range = round(lamOf(inb[1]) / lamOf(inb[2]), 1)) ``` The capacity ranges from `r round(lamOf(inb[1]),0)` down to `r round(lamOf(inb[2]),0)` across the plausible `alpha`, a `r round(lamOf(inb[1])/lamOf(inb[2]),1)`-fold spread. Every persistence verdict and every patch-value ranking from the capacity tutorial inherits that uncertainty. The honest way to report a capacity from a snapshot is as a range this wide, not a single number, unless `alpha` has been fixed from independent dispersal data. ## Check three: is the threshold a promise? The capacity rule says persist when `lambda_M` exceeds the species threshold `delta`. That is a deterministic statement about an infinite system. A real metapopulation has finitely many patches and turns over by chance, so it can wink out even when the deterministic rule says it is safe. To measure the gap, simulate a spatially explicit stochastic model on the same landscape: each occupied patch goes extinct at a rate that falls with area, each empty patch is colonised at a rate that rises with connectivity, and we ask how often the whole network dies within a fixed horizon. ```{r} #| label: extinction-curve lamM <- lamOf(alpha_t) sim_srl <- function(cc, ee, steps, tau, seed) { set.seed(seed) o <- rep(1L, n) for (t in 1:steps) { S <- conn(o, D, A, alpha_t) pc <- 1 - exp(-cc * S * tau); pe <- 1 - exp(-(ee / A) * tau) new <- o; emp <- o == 0L; occ <- o == 1L new[emp] <- as.integer(runif(sum(emp)) < pc[emp]) new[occ] <- as.integer(runif(sum(occ)) >= pe[occ]) o <- new if (sum(o) == 0) return(TRUE) } FALSE } cc <- 1; R <- 250; H <- 100; tau_s <- 0.4 ratios <- c(0.3, 0.5, 0.7, 0.9, 1.1, 1.3) pext <- sapply(ratios, function(rr) { ee <- rr * lamM * cc mean(sapply(1:R, function(k) sim_srl(cc, ee, H, tau_s, seed = 5000 + round(rr * 1000) + k))) }) data.frame(delta_over_lambda = ratios, extinction_prob = round(pext, 3)) ``` The deterministic rule is a step: safe below a ratio of one, doomed above it. The stochastic reality is a smooth ramp. At a ratio of `r ratios[2]` (that is, `lambda_M` twice the threshold, deep in the supposedly safe zone), the network still dies out `r round(100*pext[2],0)`% of the time within the horizon; at a ratio of `r ratios[3]` the extinction probability is already `r round(100*pext[3],0)`%. Comfort well below the threshold is not comfort at all. ```{r} #| label: fig-extinction #| fig-cap: "Stochastic extinction probability against the ratio of species threshold to metapopulation capacity, with the deterministic verdict as a step. The deterministic rule flips from safe to doomed at a ratio of one; the simulated extinction risk rises smoothly and is already substantial well inside the safe side." #| fig-alt: "A smooth S-shaped curve of extinction probability rising from near zero to near one as the ratio increases, overlaid on a step function that jumps from zero to one at a ratio of one. A vertical dashed line marks the threshold." #| echo: false #| warning: false dfe <- data.frame(ratio = ratios, pext = pext) step_df <- data.frame(ratio = c(0.3, 1, 1, 1.3), y = c(0, 0, 1, 1)) ggplot() + geom_vline(xintercept = 1, linetype = "dashed", colour = "#93a87f") + geom_line(data = step_df, aes(ratio, y), colour = "#93a87f", linewidth = 0.9) + geom_line(data = dfe, aes(ratio, pext), colour = "#b5534e", linewidth = 1) + geom_point(data = dfe, aes(ratio, pext), colour = "#b5534e", size = 2) + annotate("text", x = 0.62, y = 0.85, label = "deterministic verdict", colour = "#6f836a", size = 3.3) + annotate("text", x = 1.18, y = 0.35, label = "stochastic risk", colour = "#b5534e", size = 3.3) + labs(x = "threshold / capacity (delta / lambda_M)", y = "extinction probability") + theme_minimal(base_size = 13) + theme(panel.grid.minor = element_blank()) ``` The lesson generalises the Levins tutorial's stochastic detour to a real landscape: a threshold is necessary, not sufficient. A finite metapopulation needs its capacity to sit well above the threshold, not just barely above it, before persistence is anything like assured. ## What the checks can and cannot do The three diagnostics constrain the model but none of them can confirm it, and the reason is the same in each case: a snapshot is a single frame of a moving system. The equilibrium check can flag a pattern that looks disturbed, but it cannot prove a quiet-looking snapshot really is at rest rather than drifting slowly. The kernel check honestly reports that one snapshot leaves `alpha` loose, but it cannot tighten it; only independent dispersal data can. And the stochastic check shows that persistence is a probability, not a fact, so no single observation can certify a metapopulation as safe. Persistence is a claim about the future, and the only real test is more snapshots through time: watch the same network across years, and the equilibrium assumption, the turnover rates, and the persistence verdict all become checkable in a way that one map never allows. These tutorials fit what a single snapshot can give; knowing its limits is what keeps the tidy numbers honest. ## References Hanski I 1994. Journal of Animal Ecology 63(1):151-162 (10.2307/5591). Ovaskainen O, Hanski I 2002. Theoretical Population Biology 61(3):285-295 (10.1006/tpbi.2002.1586). Moilanen A 2002. Oikos 96(3):516-530 (10.1034/j.1600-0706.2002.960313.x). Lande R 1993. The American Naturalist 142(6):911-927 (10.1086/285580). Hanski I, Ovaskainen O 2000. Nature 404(6779):755-758 (10.1038/35008063). ## Related tutorials - [The Levins metapopulation model](../the-levins-metapopulation-model/) - [The incidence function model](../incidence-function-model/) - [Metapopulation capacity](../metapopulation-capacity/) - [Checking a dynamic occupancy model](../checking-a-dynamic-occupancy-model/)