# 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.