Metapopulation capacity

Compute metapopulation capacity in R as the leading eigenvalue of the landscape matrix, rank patches by conservation value, and test the persistence threshold.

The Levins model gave a persistence rule for a soup of identical patches: the species holds on when its colonisation outpaces its extinction. Real landscapes are not a soup. Patches differ in size and sit at different distances, so the question becomes whether this particular arrangement of this particular set of patches can support the species. Ilkka Hanski and Otso Ovaskainen answered it in 2000 with a single number: the metapopulation capacity. It collapses the whole geometry of a patch network into one quantity that plays the role of the Levins threshold, and it comes with a bonus: it ranks every patch by how much it contributes to persistence. This post computes it in base R on the landscape from the incidence function tutorial, ranks the patches, and shows how fragile persistence can be when the best patches go first.

From a rate to a landscape number

In the Levins world a species persists when a threshold is beaten. Hanski and Ovaskainen showed that once patches differ, the right generalisation of that threshold is an eigenvalue. Build a matrix M whose entry M_ij measures how strongly patch j can rescue or recolonise patch i: near patches with large areas exchange colonists readily, distant or tiny patches barely at all.

\[M_{ij} = e^{-\alpha\, d_{ij}}\, A_i A_j \quad (i \neq j), \qquad M_{ii} = 0.\]

The distance d_ij enters through the same dispersal kernel as before, and the areas A_i A_j weight the exchange by how much habitat sits at each end. The metapopulation capacity is the leading eigenvalue of this matrix, written lambda_M. A species with extinction-to-colonisation ratio delta persists in the landscape exactly when

\[\lambda_M > \delta.\]

The landscape supplies lambda_M; the species supplies delta. A demanding species (one that goes extinct easily or colonises weakly, so delta is large) needs a landscape with high capacity. The same landscape can hold a hardy species and lose a fragile one.

Building the matrix

We use the clustered network from the incidence function tutorial: six clusters plus scattered isolates, 200 patches, a wide spread of areas, and the same dispersal parameter.

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))
xy <- rbind(clus, iso)
n <- nrow(xy)
A <- rlnorm(n, meanlog = 0, sdlog = 0.9)
D <- as.matrix(dist(xy))
alpha <- 0.6
n

[1] 200

The capacity is one eigendecomposition away. The matrix is symmetric, so the eigenvector is real and, for a connected landscape, single-signed: every patch contributes to the leading mode.

Mmat <- function(D, A, alpha) {
  M <- exp(-alpha * D) * outer(A, A)
  diag(M) <- 0
  M
}
M <- Mmat(D, A, alpha)
ev <- eigen(M, symmetric = TRUE)
lamM <- ev$values[1]
v <- abs(ev$vectors[, 1])
v <- v / sqrt(sum(v^2))
c(lambda_M = round(lamM, 2), eigenvector_single_signed = all(ev$vectors[, 1] > 0) || all(ev$vectors[, 1] < 0))

                 lambda_M eigenvector_single_signed 
                    86.87                      1.00 

This landscape has a capacity of 86.9. On its own that number means nothing; it becomes a verdict only against a species delta. Suppose a species with delta = 58. Since lambda_M = 86.9 exceeds it, the species persists, but the margin is not huge, and the next sections show how quickly it can be eaten away.

Figure 1: The patch network, each patch sized by area and shaded by its contribution to metapopulation capacity. A handful of large, central patches carry most of the capacity; the pale scatter of small isolated patches adds almost nothing.

Which patches matter

The capacity ranks patches for free. The value of a patch is how much capacity the landscape loses when that patch is destroyed: delete row and column i from M, recompute the leading eigenvalue, and take the drop. Do it for all 200.

drop <- sapply(1:n, function(i)
  lamM - eigen(M[-i, -i], symmetric = TRUE, only.values = TRUE)$values[1])
ord <- order(drop, decreasing = TRUE)
top5 <- round(drop[ord[1:5]], 2)
share5 <- 100 * sum(drop[ord[1:5]]) / sum(drop)
trivial <- sum(drop < 0.01 * lamM)
list(top5_values = top5, top5_share_pct = round(share5, 1), patches_under_1pct = trivial)

$top5_values
[1] 20.37 17.05 15.79  6.11  5.19

$top5_share_pct
[1] 68.2

$patches_under_1pct
[1] 185

The distribution of value is savagely unequal. The five most valuable patches account for 68% of the total, and 185 of the 200 patches each contribute less than one per cent of the capacity. The single most valuable patch is number 69, with area 6.61; removing it alone drops the capacity from 86.9 to 66.5, a loss of 23%.

You do not actually need 200 eigendecompositions to find these patches. Hanski and Ovaskainen showed that a patch’s contribution is close to the square of its leading eigenvector element; the drop from removing patch i is well approximated by 2 * lambda_M * v_i^2. The approximation is nearly exact here.

approx <- 2 * lamM * v^2
c(correlation = round(cor(drop, approx), 4),
  rank_correlation = round(cor(drop, approx, method = "spearman"), 4))

     correlation rank_correlation 
          0.9992           1.0000 

The eigenvector approximation matches the exact removal drop with a correlation of 0.999 and an identical ranking (rank correlation 1). So the leading eigenvector is a ready-made map of conservation priority: patches with large v_i are the ones to protect.

Fragility: when the best patches go first

Persistence with a margin of 86.9 against 58 sounds comfortable. It is not, because habitat loss rarely falls evenly. If development or degradation takes the highest-value patches first, the capacity collapses far faster than the patch count suggests. Compare two removal orders: worst case (remove patches in order of value) against neutral (remove them at random).

delta <- 58
lam_seq <- numeric(n + 1); lam_seq[1] <- lamM; keep <- rep(TRUE, n)
for (k in 1:n) {
  keep[ord[k]] <- FALSE
  lam_seq[k + 1] <- if (sum(keep) >= 2)
    eigen(M[keep, keep], symmetric = TRUE, only.values = TRUE)$values[1] else 0
}
set.seed(11)
lam_rnd <- numeric(n + 1); lam_rnd[1] <- lamM; keepr <- rep(TRUE, n); rnd <- sample(n)
for (k in 1:n) {
  keepr[rnd[k]] <- FALSE
  lam_rnd[k + 1] <- if (sum(keepr) >= 2)
    eigen(M[keepr, keepr], symmetric = TRUE, only.values = TRUE)$values[1] else 0
}
cross_val <- which(lam_seq < delta)[1] - 1
cross_rnd <- which(lam_rnd < delta)[1] - 1
c(by_value = cross_val, at_random = cross_rnd)

 by_value at_random 
        2       126 

Removing the most valuable patches first, the capacity falls below the species threshold after only 2 of 200 patches are gone; the first three removals take it from 86.9 to 66.5, 55.9, then 50.2. Removing patches at random, the same landscape survives the loss of 126 patches before crossing. Two versus a hundred and twenty-six: the fate of the metapopulation depends far less on how much habitat is lost than on which habitat.

Figure 2: Metapopulation capacity as patches are removed, worst case (highest-value patches first) against random order. The dashed line is the species persistence threshold. Targeted loss crosses it almost immediately; random loss takes many more patches.

The honest limit: it all rides on the kernel

The capacity is only as trustworthy as the dispersal parameter fed into it. That parameter, alpha, is the one the incidence function tutorial fixed by hand because a snapshot cannot pin it down, and the capacity is sensitive to it. Recompute lambda_M across a range.

lamOf <- function(a) eigen(Mmat(D, A, a), symmetric = TRUE, only.values = TRUE)$values[1]
als <- c(0.3, 0.6, 1.0, 1.5)
lams <- sapply(als, lamOf)
persists <- lams > delta
data.frame(alpha = als, lambda_M = round(lams, 1), persists = persists)

  alpha lambda_M persists
1   0.3    158.5     TRUE
2   0.6     86.9     TRUE
3   1.0     48.9    FALSE
4   1.5     29.2    FALSE

As alpha moves from 0.3 to 1.5, the capacity slides from 158 down to 29. Against delta = 58 the verdict flips: the species persists for a short-dispersal kernel but is doomed for a long-distance one, on the same map. The reassuring part is that the ranking of patches is far steadier than the verdict: the top few patches stay near the top across kernels, so a conservation priority list is less sensitive to getting alpha wrong than a yes-or-no persistence call is. Still, any absolute claim from lambda_M inherits the uncertainty in alpha, which is exactly why the checking tutorial insists that parameter be reported with its full width, not as a point.

That coupling is the theme of the last tutorial in this set: the persistence threshold is a clean line in theory, but every ingredient that feeds it (the kernel, the equilibrium assumption, the finiteness of a real population) blurs it. A number as tidy as lambda_M earns its keep only alongside the checks that say how much to trust it.

References

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

Ovaskainen O, Hanski I 2001. Theoretical Population Biology 60(4):281-302 (10.1006/tpbi.2001.1548).

Hanski I 1999. Metapopulation Ecology. Oxford University Press. ISBN 978-0-19-854065-6.

Hanski I, Gilpin M 1991. Biological Journal of the Linnean Society 42(1-2):3-16 (10.1111/j.1095-8312.1991.tb00548.x).

Ovaskainen O, Hanski I 2003. Theoretical Population Biology 64(4):481-495 (10.1016/S0040-5809(03)00102-3).

Related tutorials

• The Levins metapopulation model
• The incidence function model
• Checking a metapopulation model
• Population viability and extinction risk