library(regnans)JJ12: the evolution of migratory arrival time
How competition can favour arrival before the seasonal optimum
What the model represents
Migratory birds compete for a fixed number of breeding territories. Arriving early can be costly because conditions are poor, but it also improves an individual’s chance of securing a territory. Arriving at the seasonal optimum would maximise reproduction if there were no competition.
Competition changes that result. Each individual gains by arriving slightly earlier than its competitors, so evolution favours an arrival date before the date that would maximise total population output. This difference between the individual and population optima is an evolutionary tragedy of the commons.
The ecological setting follows Johansson and Jonzén (2012), Game theory sheds new light on ecological responses to current climate change when phenology is historically mismatched. The equations below use the simplified analytic form described by Brännström, Johansson, & von Festenberg (2013), which gives us an exact answer against which to check the numerical solver.
Fitness equations
Competitive ability declines with arrival time x, while reproduction follows a Gaussian curve centred on the seasonal optimum x_opt:
\[ C(x) = e^{-a x}, \qquad R(x) = R_0\,\exp\!\left(-\frac{(x-x_{\mathrm{opt}})^2}{2\sigma^2}\right). \]
With year-to-year survival p, a single resident reaches \(n^\* = K R(x)/(1-p)\). The invasion fitness of a mutant x' (a discrete-time geometric growth rate; \(w = 1\), i.e. log = 0, at the resident) is
\[ w(x') = (1-p)\,\frac{C(x')\,R(x')}{C(x)\,R(x)} + p . \]
The selection gradient is a decreasing straight line, so the model has one continuously stable strategy and no evolutionary branching. Its exact value is
\[ x^\* = x_{\mathrm{opt}} - a\,\sigma^2 , \]
which lies before the population optimum x_opt whenever a > 0.
Example analysis with regnans
a <- 0.125; x_opt <- 0; sigma <- 10
h <- harness_jj12(a = a, x_opt = x_opt, sigma = sigma, R0 = 1, K = 5, p = 0.5)
out <- community_start(bounds(x = c(-40, 40)), harness = h, trait_scale = "linear") |>
community_solve_singularity()
c(numeric = as.numeric(out$traits),
analytic = x_opt - a * sigma^2) # should match: -12.5 numeric analytic
-12.5 -12.5 The numerical solution should match the exact result. That known answer makes this model a useful check on the solver. With these parameters, the evolved strategy arrives a*sigma^2 = 12.5 time units before the seasonal optimum.
Selection gradient and convergence
The gradient is positive for arrival times below x* and negative for arrival times above it. Small evolutionary steps therefore move the population towards x* from either side:
xs <- seq(-30, 10, length.out = 60)
grad <- sapply(xs, function(x) {
community_start(bounds(x = c(-40, 40)), harness = h, trait_scale = "linear") |>
community_add(trait_matrix(x, "x")) |>
community_demography() |>
community_selection_gradient() |>
(\(comm) comm$selection_gradient)()
})
plot(xs, grad, type = "l", lwd = 2,
xlab = "resident arrival time x", ylab = "selection gradient")
abline(h = 0, lty = 3); abline(v = x_opt - a * sigma^2, col = "forestgreen", lty = 2)
legend("topright", legend = "x* = x_opt - a sigma^2", col = "forestgreen",
lty = 2, bty = "n")
Mismatch grows with the competitive advantage a
Increasing a strengthens the benefit of arriving before competitors. The evolved arrival date therefore moves farther ahead of the seasonal optimum:
sapply(c(0, 0.04, 0.125), function(a) {
h <- harness_jj12(a = a, x_opt = 0, sigma = 10)
as.numeric(community_start(bounds(x = c(-40, 40)), harness = h, trait_scale = "linear") |>
community_solve_singularity() |>
(\(comm) comm$traits)())
}) # 0, -4, -12.5 — earlier arrival as competition for territories intensifies[1] 0.0 -4.0 -12.5Following x* as a parameter changes
The calculation above re-solves the model at each value of a. regnans can instead follow the singular strategy along a parameter: community_continue_singularity() steps through a series of values, predicts where \(x^\*\) will be at each from its sensitivity \(dx^\*/dp\), solves from that prediction, and classifies the result. The method is described on the DD99 page.
The parameter enters through a function that takes a value and returns a community built at it. community_parameter_map() writes one for any argument of the model, but you can write your own, and then the parameter can be anything you can build a community from. Here it is the variance \(\sigma^2\) of the reproduction curve, which is not itself an argument of harness_jj12():
by_sigma2 <- function(s2) {
community_start(bounds(x = c(-40, 40)), trait_scale = "linear",
harness = harness_jj12(a = a, x_opt = x_opt, sigma = sqrt(s2)))
}
sigma2 <- seq(25, 200, by = 25)
path <- community_continue_singularity(community_solve_singularity(by_sigma2(sigma2[1])),
by_sigma2, sigma2)
path<singularity_path: 8 points, p from 25 to 200>
classification: CSS (p 25 to 200)
equilibrium solves: 40round(cbind(sigma2 = path$p, x = path$traits[, 1], exact = x_opt - a * path$p,
dx_dsigma2 = path$sensitivity[, 1]), 4) sigma2 x exact dx_dsigma2
[1,] 25 -3.125 -3.125 -0.125
[2,] 50 -6.250 -6.250 -0.125
[3,] 75 -9.375 -9.375 -0.125
[4,] 100 -12.500 -12.500 -0.125
[5,] 125 -15.625 -15.625 -0.125
[6,] 150 -18.750 -18.750 -0.125
[7,] 175 -21.875 -21.875 -0.125
[8,] 200 -25.000 -25.000 -0.125Since \(x^\* = x_{\mathrm{opt}} - a\sigma^2\) is a straight line in \(\sigma^2\), its slope is \(dx^\*/d\sigma^2 = -a\) = -0.125 at every point. The prediction from each point’s slope therefore lands on the next singular strategy, and the solver has nothing left to correct. The strategy is a CSS all along the path:
plot(path$p, path$traits[, 1], type = "n",
xlab = expression("variance of the reproduction curve, " * sigma^2),
ylab = "evolved arrival time x*")
abline(a = x_opt, b = -a, lty = 2, col = "grey40")
points(path$p, path$traits[, 1], pch = 16, col = "#2a78d6")
legend("topright", legend = c("CSS, followed by regnans", expression(x^"*" == x[opt] - a * sigma^2)),
pch = c(16, NA), lty = c(NA, 2), col = c("#2a78d6", "grey40"), bty = "n")
A map from community_parameter_map() can vary several arguments at once. community_parameter_sensitivity() then gives the sensitivity of \(x^\*\) to each, here at the singular strategy solved at the top of this page. Differentiating \(x^\* = x_{\mathrm{opt}} - a\sigma^2\) gives \(-\sigma^2\), \(1\) and \(-2a\sigma\):
S <- community_parameter_sensitivity(out, community_parameter_map(out, c("a", "x_opt", "sigma")))
rbind(numeric = S[1, ], analytic = c(-sigma^2, 1, -2 * a * sigma)) a x_opt sigma
numeric -100 1 -2.500005
analytic -100 1 -2.500000References
Brännström, Åke, Johansson, J., & von Festenberg, N. (2013). The hitchhiker’s guide to adaptive dynamics. Games, 4(3), 304–328. doi:10.3390/g4030304