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. The rule asks whether the capacity clears that threshold, and the next sections take patches out of the landscape and watch the comparison move.

Figure 1: The patch network, each patch sized by area and shaded by its contribution to metapopulation capacity. A handful of large, well-connected 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 a rank correlation of 0.99999, which the four-decimal printout above shows as 1.0000; the two orderings differ only by a few near-adjacent swaps. 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

Habitat loss rarely falls evenly. If development or degradation takes the highest-value patches first, the capacity falls 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 is below the species threshold once 2 of 200 patches have 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 is below the threshold after 126. The capacity depends on which patches remain, not only on how many, so a count of habitat lost says little on its own.

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.

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 can flip within that range: a small alpha is a shallow decay, so the same map can support the species under a long-distance kernel and not under a short-distance one. The reassuring part, which the table above does not show, 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
• Long-distance jumps and stratified spread