Automating adaptive dynamics even more

Julia as a CAS

using Plots
using Symbolics
using Latexify

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. Almost a year later, it’s a good idea to revisit this topic by using Symbolics, as an alternative to Zygote. I strongly suggest you read the previous post, to freshen up on the model and the notation.

As a reminder, in our population growth model, the invasion fitness for a mutant with trait $y$ is

$$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)$$

We can write this as a Symbolics expression:

@variables x y sk sc
sxy = 1 - exp((y^2 - x^2) / (2sk^2) - ((x - y)^2) / (2sc^2))
1 - exp((y^2 - (x^2)) / (2(sk^2)) + (-((x - y)^2)) / (2(sc^2)))

What can we do with this? Well, we can make this into a function, and specifically use a generator to set the values of $\sigma_c$ and $\sigma_k$. So let’s go. We will specifically replace the symbols sc and sk by their numerical values, so that we can build a function of $(x, y)$. This is all done with a call to substitute, build_function, and eval:

function invasion_fitness_generator(; sigma_k::Float64=0.4, sigma_c::Float64=0.2)
    valued_expression = substitute(sxy, Dict(sk => sigma_k, sc => sigma_c))
    ifg = eval(build_function(valued_expression, [x, y]))
    return (x, y) -> ifg([x, y])
end
invasion_fitness_generator (generic function with 1 method)

Let’s try with a first set of conditions – the invasion fitness should be 0 (note that we use the default values of $\sigma_c$ and $\sigma_k$):

invasion_fitness_generator()(0.0, 0.0)
0.0

Good - now let’s plot this around the $x=0$ point. We can define two functions for the invasion fitness, one where coexistence is possible, and one where it isn’t:

coexistence = invasion_fitness_generator(; sigma_k=0.4, sigma_c=0.2)
exclusion = invasion_fitness_generator(; sigma_k=0.2, sigma_c=0.4)
#2 (generic function with 1 method)

The plot becomes:

plot(
    (y) -> coexistence(0.0, y); xlim=(-0.1, 0.1), dpi=600, frame=:origin, lab="Coexistence"
)
plot!((y) -> exclusion(0.0, y); lab="No coexistence")
xaxis!("Mutant trait")
yaxis!("Invasion fitness")

We can very easily define a matrix of values, representing the PIP around the $x=y=0$ point. Let’s remember that $x$ is the rows, and $y$ the columns, so we will need to pay attention to this in our figure:

v = LinRange(-0.5, 0.5, 200)
heatmap(
    v,
    v,
    coexistence.(v, v') .> 0.0;
    dpi=600,
    c=:Blues,
    cbar=false,
    aspectratio=1,
    frame=:grid,
)
xaxis!(extrema(v), "Resident trait")
yaxis!(extrema(v), "Mutant trait")
vline!([0.0]; c=:white, lw=2, ls=:dash, lab="")
hline!([0.0]; c=:white, lw=2, ls=:dash, lab="")

This figure is the PIP for the invasion of $y$ in a resident $x$, so if we want the coexistence region, we need to flip it and overlay it:

v = LinRange(-0.5, 0.5, 200)
pipx = coexistence.(v, v') .> 0.0
pipy = permutedims(pipx)
heatmap(v, v, pipx .* pipy; dpi=600, c=:Greens, cbar=false, aspectratio=1, frame=:grid)
xaxis!(extrema(v), "Resident trait")
yaxis!(extrema(v), "Mutant trait")
vline!([0.0]; c=:white, lw=2, ls=:dash, lab="")
hline!([0.0]; c=:white, lw=2, ls=:dash, lab="")

Fantastic! This coexistence diagram shows that the regions both vertically and horizontally around the singular strategy give positive invasion exponents, and so there is a possible branching point. In order to classify this singularity, we need to evaluate two things: ∂²s/∂x², and ∂²s/∂y². Let’s define a partial differentiation operator using the functions from Symbolics:

∂x = Differential(x)
∂y = Differential(y)
(::Symbolics.Differential) (generic function with 2 methods)

We can apply this to the expression of $s_x(y)$:

∂²s∂x² = (∂x ∘ ∂x)(sxy)
Differential(x)(Differential(x)(1 - exp((y^2 - (x^2)) / (2(sk^2)) + (-((x - y)^2)) / (2(sc^2)))))
∂²s∂y² = (∂y ∘ ∂y)(sxy)
Differential(y)(Differential(y)(1 - exp((y^2 - (x^2)) / (2(sk^2)) + (-((x - y)^2)) / (2(sc^2)))))

This is a big expression (if we pass it to expand_derivatives), but remember that we want to evaluate it at $x = 0, y=0$. This can be done with the right call to substitute:

substitute(expand_derivatives(∂²s∂x²), Dict(x => 0, y => 0))
-(-1 / (sc^2) + -1 / (sk^2))
substitute(expand_derivatives(∂²s∂y²), Dict(x => 0, y => 0))
-(-1 / (sc^2) + 1 / (sk^2))

Sweet! Let’s rewrite this as equations:

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

$$\frac{\partial^2s}{\partial x^2} = \sigma_c^{-2} - \sigma_k^{-2}$$

In order for the equilibrium to lead to a dimorphism, it must be unstable, so the second partial derivative on $x$ must be positive: this is trivially true all the time, as both parameters are positive.

With these two equations, we can also check that whenever there is coexistence (as explained in the previous post), this system will lead to a branching point when $\sigma_k > \sigma_c$.

There is a lot more we could do using Symbolics, but this is a nice illustration. Having both Symbolics and Zygote in our toolkit is a very nice thing to have, and enables a complementarity of approaches in addition to solving differential equations.