Automating adaptive dynamics some more

In the previous entry, I gave an illustration of how Zygote can help with the basic operations of adaptive dynamics, and notably how it can make pairwise invasibility plot (PIP) easy to produce. Today, I will go over a slightly different topic, which is the trait evolution plot.

Making a trait evolution plot requires (i) to produce the PIP, (ii) to overlap the PIP to its mirror image to identify the region of possible coexistance, and then (iii) to identify the isoclines alongside which the selection gradient vanishes.

Setting up the model

To illustrate this, let’s work on a simple model where species grow logistically and compete with one another. We can write this model as:

$$\dot n_x = r n_x \times \left[1 - \frac{1}{k_0\times f_k(x)}(n_xf_c(x,x))\right]$$

We will assume that the maximal carrying capacity is 1, and define $f_k$ as

$$f_k(x) = \text{exp}\left(-\frac{x^2}{2\sigma_k^2}\right)$$

We will use a competition function such that the intra-specific competition is 1, and the inter-specific competition decreases when the traits become more different:

$$f_c(x,y) = \text{exp}\left(-\frac{(x-y)^2}{2\sigma_c^2}\right)$$

We will also assume that $r = 1$, so we can rewrite our model as

$$\dot n_x = n_x \times \left(1 - \frac{f_c(x,x)}{f_k(x)}n_x\right)$$

Getting the invasion fitness

This gives an invasion fitness for a mutant with trait $y$ of

$$s_x(y) = \frac{\dot n_y}{n_y} = 1 - \frac{f_c(x,y)}{f_k(y)}n_x^\star$$

We can get the value of $n^\star_x$ from above, which is $f_k(x)$, and this finally gives us the invasion fitness:

$$s_x(y) = \frac{\dot n_y}{n_y} = 1 - \frac{f_c(x,y)}{f_k(y)}f_k(x)$$

Let’s re-write this in a more logical way:

$$s_x(y) = 1 - \frac{f_k(x)}{f_k(y)}f_c(x,y)$$

Before we do anything, let’s expand the functions:

$$s_x(y) = 1 - \frac{\text{exp}\left(-\frac{x^2}{2\sigma_k^2}\right)}{\text{exp}\left(-\frac{y^2}{2\sigma_k^2}\right)}\times \text{exp}\left(-\frac{(x-y)^2}{2\sigma_c^2}\right)$$

This actually makes sense, because this is of the form $1-e^{a}/e^{b}\times e^{c}$, so we can rewrite this as $1-e^{a-b+c}$, which in our case is

$$s_x(y) = 1 - \text{exp}\left(-\frac{x^2}{2\sigma_k^2}+\frac{y^2}{2\sigma_k^2}-\frac{(x-y)^2}{2\sigma_c^2}\right)$$

This finally reduces to

$$s_x(y) = 1 - \text{exp}\left(\frac{y^2-x^2}{2\sigma_k^2}-\frac{(x-y)^2}{2\sigma_c^2}\right)$$

Ceci n’est pas une PIP (yet)

This expression is enough to think about the PIP, because we really only care about $s_x(y)>0$, and this happens when the exponent is smaller than 1, which requires the inside of the exponent to be negative, which finally requires

$$\frac{y^2-x^2}{2\sigma_k^2} > \frac{(x-y)^2}{2\sigma_c^2}$$

In other words, the isoclines in the PIP will be given by

$$\frac{y^2-x^2}{2\sigma_k^2} = \frac{(x-y)^2}{2\sigma_c^2}$$

which is trivially true when $y = x$, and somewhat-less-trivially true (but c’mon) when

$$y = x\frac{\sigma_k^2+\sigma_c^2}{\sigma_k^2-\sigma_c^2}$$

So, there is our first important result: if we want to see any coexistence at all, we will need to have a positive isocline, and so $\sigma_k > \sigma_c$, which is to say the species must tolerate being away from its optimum more than it tolerates competition.

Implementing the model

We can write the above as a function, which we will put through various manipulations to figure our when more than one morph can persist.

function s(x,y; σc=0.6, σk=0.7)
  left = (y^2-x^2)/(2σk^2)
  right = ((x-y)^2)/(2σc^2)
  return 1.0 - exp(left - right)
end

s (generic function with 1 method)

Building the PIP

The first step is going to be looking at the PIP, to see which mutants can invade a resident at equilibrium:

traits = LinRange(-4.0, 4.0, 450)
pip = zeros(Float64, (length(traits), length(traits)))
for (i,r) in enumerate(traits)
  for (j,m) in enumerate(traits)
    pip[j,i] = s(r, m)
  end
end
heatmap(traits, traits, pip.>0, aspectratio=1, frame=:box, c=:Greens, dpi=300, legend=false)
xaxis!(extrema(traits), "Resident")
yaxis!(extrema(traits), "Mutant")

So far, so good – this PIP goes through positive values both vertically and horizontally around the singular strategy, so it is convergence stable (left to its own devices, a single morph will evolve towards $x=0$), but not evolutionary stable. Let’s look at the region of coexistence.

Defining the region of coexistence

We can now flip this plot to see if the resident can invade the mutant back - this is an important step, because coexistence requires mutual invasibility from rarity.

heatmap(traits, traits, permutedims(pip).>0, aspectratio=1, frame=:box, c=:Oranges, dpi=300, legend=false)
xaxis!(extrema(traits), "Resident")
yaxis!(extrema(traits), "Mutant")

Finally, the region of possible coexistence is given by the overlap between both plots, i.e. the values of traits for which both morphs can invade one another from rarity.

coex = (pip.>0).*(permutedims(pip).>0)
heatmap(traits, traits, coex, aspectratio=1, frame=:box, c=:Purples, dpi=300, legend=false)
xaxis!(extrema(traits), "Resident")
yaxis!(extrema(traits), "Mutant")

But from here, where do we go?

Well, the trait evolution plot requires that we look at the points where the gradients vanish. So let’s do this. As an important note, we do not need to evaluate the gradient outside of the coexistence region.

gr1 = zeros(Float64, size(pip))
gr2 = zeros(Float64, size(pip))
for (i,r) in enumerate(traits)
  for (j,m) in enumerate(traits)
    if coex[j,i]
      gr1[j,i], gr2[j,i] = gradient((r,m) -> s(r,m), r, m)
    else
      gr1[j,i] = gr2[j,i] = NaN
    end
  end
end

pl_r = contour(traits, traits, gr1, fill=true, c=:PRGn, frame=:box, aspectratio=1)
xaxis!(pl_r, extrema(traits), "Resident")
yaxis!(pl_r, extrema(traits), "Mutant")
title!(pl_r, "Resident trait")

pl_m = contour(traits, traits, gr2, fill=true, c=:PiYG, frame=:box, aspectratio=1)
xaxis!(pl_m, extrema(traits), "Resident")
yaxis!(pl_m, extrema(traits), "Mutant")
title!(pl_m, "Mutant trait")

plot(
  pl_r,
  pl_m,
  dpi=300
)

Let’s finally show the direction of trait change, by measuring the norm of the two vectors:

plot(frame=:box, aspectratio=1, dpi=300)
for r in LinRange(minimum(traits), maximum(traits), 32)
  for m in LinRange(minimum(traits), maximum(traits), 32)
    if (s(r,m)>0) & (s(m,r)>0)
      dr, dm = gradient((r,m) -> s(r,m), r, m).*0.02
      plot!([r, r+dr], [m, m+dm], lab="", arrow=:arrow, c=:darkgrey)
    end
  end
end
xaxis!(extrema(traits), "Resident")
yaxis!(extrema(traits), "Mutant")

The PIP (and other visualisations) tells us that the system is convergence stable (we can follow the grey arrows to the $(0,0)$ point in the above figure), but is not evolutionary stable, which may be a branching point. We can also see that after some limit, the gradient vanishes (there are no arrows drawn), which suggests that the branching process will stop somewhere along these lines.

Deep dive into the singularity

At this point, we need to check three things. First, evolutionary stability, which requires $\partial^2/\partial y^2 s < 0$:

f = (x) -> s(0.0, x)
@assert f''(0.0) > 0

Our system is not evolutionary stable, and this is what we expected. Second, we need to check convergence stability, which requires $\partial^2/\partial x^2 s > \partial^2/\partial y^2 s$:

g = (x) -> s(x, 0.0)
@assert g''(0.0) > f''(0.0)

And finally, the protected dimorphism requires $\partial^2/\partial y^2 s> - \partial^2/\partial x^2 s$:

@assert f''(0.0) > - g''(0.0)

And there we go – whenever our system allows coexistence, the singular strategy can be reached by convergence, is unstable, and leads to a protected dimorphism: branching will take place and two new morphs will be created.

Conclusion

In conclusion, Zygote shreds. The last bit, i.e. the evaluation of the second order partial derivatives, is simply mind blowing. That we can define a function f and then have magically f'' be a function that gives the corect result is… game changing. This is a really good tool to sense-check some simulations as well – implementing models of adaptive dynamics has all sort of tricks, and having a tool to facilitate the model analysis gives another way to check the code.