To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 56bca353 authored by Victor's avatar Victor
Browse files

adding docs

parent d48ba3ed
Pipeline #74670 passed with stage
in 23 minutes and 26 seconds
docs/build docs/build
docss/.git
docs/*.png docs/*.png
docs/uploads docs/uploads
docs/.gitignore
using Documenter, ABMEv
# push!(LOAD_PATH,"/Users/victorboussange/ETHZ/projects/ABMEv/") # not sure this is necessary
pathsrc = joinpath(@__DIR__,"src")
makedocs(sitename="ABMEv.jl",
format = Documenter.HTML(prettyurls = false),
pages = [
"Home" => "index.md",
"Manual" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"manual"),join=true))],
"Mathematics" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"mathematics"),join=true))],
"Examples" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"examples"),join=true))],
"Library" => Any[
"Public" => "lib/public.md",
# "Internals" => map(
# s -> "lib/internals/$(s)",
# sort(readdir(joinpath(@__DIR__, "src/lib/internals")))
# ),
],
"Developping" => [ joinpath(s[end-1:end]...) for s in splitpath.(readdir(joinpath(pathsrc,"dev"),join=true))],
# "contributing.md",
],)
deploydocs(repo= "gitlab.ethz.ch:bvictor/abmev.wiki.git",)
# Developping the code
I recommend to first clone your branch in the directory you like best, and then to
To develop, you ca
```julia
using Pkg
Pkg.dev("path_to_ABMEv_dir")
```
You can also do the same trick with directly the gitlab address, cf [Pkg.jl](https://docs.julialang.org/en/v1/stdlib/Pkg/index.html)
## Future directions
- Try to improve parallelism with the help of Threads.@Threads and @inbounds (cf [tutorial](https://juliagpu.gitlab.io/CUDA.jl/tutorials/introduction/#Introduction-1) )
- Make use of CUDA.jl to accelerate the simulations wih the use of GPU.
## Todo
```julia
abstract type AbstractAgent{A<:Ancestors,R<:Rates} end
```
It seems that we do not need to have the `Rates{}` parameter for `Agent` type.
## Notes
### Numerics
:warning:
[Don’t pass expressions, or strings, pass functions](https://www.youtube.com/watch?v=mSgXWpvQEHE&t=573s)
### Birth and Death mechanisms
<!-- > We are always balanced between taking the integral of the competition and resource kernel as constant, or taking its maximum peak as constant. -->
<!-- :poop: -->
### Geotrait
The geotrait is calculated *a posteriori*, and is not taken into account during the simulation.
> It used to be but for the sake of simplicity we now forget about it.
# Test
## Gillepsie
In Champagnat - HDR Memoire
```
p = 0.1
K = 1000
sigm_mut = 0.01
sigma_alpha = 0.7
sigma_K = 0.9
```
# Gillepsie simulation
## Diversification
### Gaussian birth coefficient, Constant carrying capacity
```julia
a = 0;
sigma_K = .9;
sigma_a = .7;
K0 = 1000;
K(X) = gaussian(X[1],0.,sigma_K)
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
# α(X,Y) = 0.
p_default = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 2000,
"NMax" => Int(2000),
"dt_saving" => 10.0)
na_init = K0
world0 = new_world_G(na_init,p_default,spread = .01)
```
![ABMEv_div_t4000](uploads/3a49ff4fe4db161bf360eea97694ff26/ABMEv_div_t4000.png)
### Constant birth coefficient, Gaussian carrying capacity
```julia
a = 0;
sigma_K = .9;
sigma_a = .7;
K0 = 1000;
K(X) = 1
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0 / gaussian(X[1],0.,sigma_K)
# α(X,Y) = 0.
p_default = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 4000,
"NMax" => Int(2000),
"dt_saving" => 20.0)
na_init = K0
```
![ABMEv_bis_div.ong](uploads/8e1f821923afd74902b3ec6567a1736d/ABMEv_bis_div.ong.png)
> what you could do would be to plot the adaptive dynamics of the monomorphic populations
### Quadratic birth rate
```julia
a = 0.125;
sigma_a = 1.251;
K0 = 1000;
K(X) = 1 - a * X[1]^2
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
# α(X,Y) = 0.
p_default = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 4000,
"NMax" => Int(2000),
"dt_saving" => 20.0)
na_init = K0
```
![Unknown](uploads/e6b6bc135fbd2f48cc5ad8bd18854420/Unknown.png)
#### Equivalence
![PDE_quad_termsol](uploads/c1f5c670a9d82df1349ed473b9954135/PDE_quad_termsol.png)
![ABMEv_quad_time_average_distrib_deep_time](uploads/2613b1f5e919fdbc7ee3386c680a1908/ABMEv_quad_time_average_distrib_deep_time.png)
## No diversification
```Julia
a = 0;
sigma_K = .9;
sigma_a = 1.0;
K0 = 1000;
K(X) = gaussian(X[1],0.,sigma_K)
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
# α(X,Y) = 0.
p_default = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 2000,
"NMax" => Int(2000),
"dt_saving" => 10.0)
na_init = K0
```
![Gillepsie_results_nodiversification](uploads/3a459012508c85cf853246a37537f160/Gillepsie_results_nodiversification.png)
# ABMEv.jl Documentation
This is a suite for simulating an Agent Based Model that captures the evolutionary dynamics of a population in a multidimensional space.
## Getting started
```julia
using ABMEv
```
## The `Agent` structure
This package is an Agent Based Model, where the atomic structure is the `Agent`. It has four attributes
- the ancestors' history of traits, and the corresponding time where the traits have changed,
- a death rate and a birth rate.
```julia
mutable struct Agent{T,U}
# history of traits for geotraits
x_history::Array{U}
# birth time of ancestors
t_history::Array{U,1}
# death rate
d::Float64
#birth rate
b::Float64
end
```
For more understanding of the composite types `T,U`, check the wiki: Gillepsie / Agent Type.
## Parameters of the simulation
Parameters are stored in the parameter dictionary `p`
### General parameters
- ```"reflected"=>true```: if ```true``` then reflection occurs on the first trait -which should stand for geographic position. Depending on the agent type, reflections occurs in the domain $` [-1,1] `$ or between nodes 1 and `p["nodes"]`
- ```"alpha" => α```: is the competition function
- ```"K" => K```: is the birth rate
- ```"tend" => 1.5```: is the time to end simulation
:warning: Check how to define functions α and K in the algorithm section.
### Mutation
If anisotropy in mutation, the following parameters should be declared as arrays where each entry corresponds to a dimension.
- ```mu``` The probability of mutation.
- ```D``` If mutation happens on the agent, the increment follows $`\mathcal{N}_{ 0, D}`$
### Birth
#### Growth
- ```K``` is the birth coefficient ( $`b(x) = K(x)`$ )
### Death
#### Competition
- Competition between agent with trait ```x``` and ```y``` is defined as
```α(x,y)```
- Death coefficient is defined as $`d(x^{(i)}) = \sum_j^{N(t)} \alpha(x^{(i)},x^{(j)})`$
### Fitness
Fitness is defined internally as ```b - d```.
> TODO ```b``` here is confounded with ```K```.
## Launching simulation
Two type of simulation algorithm can be used
```@contents
Pages = ["manual/gillepsie.md",
"manual/wright_fisher.md"]
Depth = 5
```
### Gillepsie algorithm
```julia
using ABMEv
## Defining parameters
sigma_K = 1.; #bandwith of resource
sigma_a = 1.2; #bandwith of competition
K0 = 1000 #carrying capacity
K(X) = gaussian(X[1],0,Float32(sigma_K)) #birth
alpha(X,Y) = gaussian(X[1],Y[1],Float32(sigma_a)) / Float32(K0) #competition
D = [1e-2] #mutation range
mu = [1.] #probability of mutation
NMax = 2000 #number of individual
dt_saving = 1.0 #time step saving
tend = 1000.
using UnPack
p = Dict{String,Any}()
@pack! p = K,alpha,D,mu,NMax,dt_saving,tend
## Initial conditions
agent0 = [Agent(.1 .* randn(Float32,1)) for i in 1:K0]
world0 = vcat(agent0[:],repeat([missing],Int(p["NMax"] - K0)))
## launch simulation
worldall,p["tspan"] = runWorld_store_G(p,world0);
using Plots
Plots.plot(worldall,p)
```
### Wright Fisher algorithm
```julia
sigma_K = .9;
sigma_a = 1.251;
K0 = 1000;
# K(X) = gaussian(X[1],0.,sigma_K)
K(X) = 1 - 0.125 * sum(X.^2)
α(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
# α(X,Y) = 0.
p = Dict(
"alpha" => α,
"K" => K,
"D" => [1e-2],
"mu" => [.1],
"tend" => 10.)
na_init = K0
agents = [Agent( [1e-2] .* randn(1) .- .5) for i in 1:K0]
@time worldall_test,p["tspan"] = runWorld_store_WF(p,agents,mode="std");
```
## References
- [Champagnat and Ferriere founding article](https://linkinghub.elsevier.com/retrieve/pii/S0040580905001632)
- [Champagnat and Ferriere second article - 2008](https://www.tandfonline.com/doi/full/10.1080/15326340802437710)
## Similar packages:
[Agents.jl](https://juliadynamics.github.io/Agents.jl/) It would be worth to have a look! It has been designed by Ali R. Vahdati, from UZH.
# Agent properties
```@autodocs
Modules = [ABMEv]
Pages = ["ABMEv_Agent.jl"]
```
# Diversity measures
```@autodocs
Modules = [ABMEv]
Pages = ["ABMEv_metrics.jl"]
```
# Gillepsie algorithm
## Mathematical foundations
- The original article by Gillepsie:
[**A general method for numerically simulating the stochastic time evolution of coupled chemical reactions**](https://www.sciencedirect.com/science/article/pii/0021999176900413?via%3Dihub)
### Rates
Each individual is assigned a birth $` b_i `$ and death $` d_i `$ rate. The total rate is given by the sum of all individual rates
```math
R(t) = \left[ \sum_i b_i(t) + d_i(t) \right]
```
A particular event, birth or death, is chosen at random with a probability equal to the rate of this event divided by the total rate $`R`$
> This has to be checked, we are still not hundred percent sure
### Time steps
An event is exponentiallly distributed in time, with parameter $`\lambda = U(t)`$. This makes events memoryless, meaning that the probability of having a birth or death event is always the same, no matter when ($`P(X > s_t | X > t) = P(X > s) `$.
> Let $`B(t) = \sum_i b_i(t)`$ and  $`D(t) = \sum_i d_i(t)`$. Let $`T_b, T_d`$ the time for a birth or death event to occur. Then we have $`P(T_b < T_d) = \frac{B(t)}{B(t) + D(t)}`$ (competing exponentials).
#### Inversion method
Let $`U`$ be an $`\mathcal{U}_{(0,1)}`$-distributed random variable and $`F \colon \R \to [0,1]`$ be a distribution function. Then we have
```math
P(I_F(U) \leq x ) = P(U \leq F(x)) = F(x)
```
Thanks to the ***inversion method*** we get the incremental time step $`dt`$, exponentially distributed with parameter $`\lambda = R(t)`$, as
```math
dt(\omega) = -\frac{\log(U(\omega))}{R(t)} \iff X(\omega) = \exp(-U(t)dt(\omega))
```
## Initialize
```@docs
new_world_G
```
## Run
```@docs
runWorld_store_G
```
## Scenarios
As of now, no mode is implemented. For further examples, check the folder `examples` in source code.
## Specific parameters
- ```dt_saving = 10.```
will allow to save the world every 10. time steps. If not specified, the algorithm will return first and last time step world.
- ```NMax``` Maximum number of individuals that can be attained. If attained, then the programm stops.
## Developping
### Efficiency
The simulation are still very long.
:flushed: How to improve it?
- We think it would be more efficient if we found an other way of incrementing mutations
### Parallelism
> For now there is no parallelism implemented for one run
> But we think we should rather set up a pmap or the macro `@Threads` to explore parameter space
```@autodocs
Modules = [ABMEv]
Pages = ["ABMEv_Gillepsie.jl"]
```
# Plotting
ABMEv comes with Plot recipes:
`plot(world::Array{U},p;what=["x","H"],trait = 1) where U <: Union{Missing,Agent{T}} where T`.
An example to use it:
```julia
using Plots;pyplot()
Plots.plot(worldall,p_default,what=["x"],trait=2)
```
You can specify what you want to plot in the array ```what```:
- ```"x"``` returns a scatter plot `(xaxis = time, yaxis = trait)` . Geotrait corresponds to `trait=0`
- `"xs"` only works for agent type `MixedAgent`, because it needs a discrete geographical space. It returns a scatter plot `(xaxis = geographical component, yaxis = trait value)`. Very similar to a `histogram2d` plot, with nicer look.
- ```"gs"``` returns a scatter plot `(xaxis = geotrait, yaxis = trait value)`
- ```"3dgeo"``` plots a 3d diagram with x axis as geotrait and y axis as the second component. :warning: this is probably weekend
- ```"3d"``` plots a 3d diagram with first and second component as x and y axis
- ```"var"``` plots the variance of the component specified by ```trait=2``` :question: with respect to time?
- ```"vargeo"``` plots the variance of the geotrait
```@autodocs
Modules = [ABMEv]
Pages = ["ABMEv_plotting.jl"]
```
# Run the World
For now three algorithms
```@docs
Gillepsie
CFM
WF
```
> Warning : `WF` not implemented yet
```@autodocs
Modules = [ABMEv]
Pages = ["ABMEv_runworld.jl"]
```
# Simulation
```@autodocs
Modules = [ABMEv]
Pages = ["ABMEv_Sim.jl"]
```
# Space
# Utils
```@autodocs
Modules = [ABMEv]
Pages = ["ABMEv_Space.jl"]
```
# Utils
```@autodocs
Modules = [ABMEv]
Pages = ["ABMEv_utils.jl"]
```
# World
```@autodocs
Modules = [ABMEv]
Pages = ["ABMEv_world.jl"]
```
# Wright Fisher algorithm
## Foundations
The Wright Fisher process is an individual based model where the number of agents is constant through time. It is helpful to visualize it through marbles in jar:
![alt text](https://upload.wikimedia.org/wikipedia/commons/0/0b/Random_sampling_genetic_drift.svg)
At each time step, $`N`$ agents are picked up from previous generation to reproduce. Their number of offspring is proportional to their fitness, calculated as usual with **birth and death rates**.
It takes thus **only one time step to go trough one generation**. Thus it is more suit- able for numerical simulations. In practice, the Moran and Wright–Fisher models give qualitatively similar results, but genetic drift runs twice as fast in the Moran model.
From this perspective we can easily get that branches are less stable than in the Gillepsie scenario, for as as time goes to infinity the probability of going extinct is intuitively bigger than 0.
## Initial conditions
One need to construct the world as an array of agents, which will be the ancestors of the following
```julia
agents = [Agent( [2e-1] .* randn(1)) for i in 1:K0]
```
The function
is then called `p["tend"] -1` times.
## Scenarios
You have several options available concerning the resource implemented and competition:
- ``` mode="std"``` is the standard mode
- ``` mode="grad2D"``` corresponds to a an ecological gradient
>We are not sure whether this corresponds to the following two images
- ``` mode="mountain"``` corresponds to a scenario where a mountain arises (with an ecological gradient)
- ``` mode="split"``` corresponds to a scenario where the resource is splitted in two
- ``` mode="graph"``` this guy is probably not working
## Parallelism
You can run your script in parallel, which makes sense for large populations. To do so:
```julia
using Distributed;addprocs()
@everywhere using ABMEv
```
Parallelism only works with Wright Fisher model.
```@autodocs
Modules = [ABMEv]
Pages = ["ABMEv_WF.jl"]
```
# The geotrait
> This note has been taken from `2Y/articles/geotrait/mathematical_notes.md`
## Julia accessors
An agents of type `Agent{Ancestors{true}}` stores the values of its ancestors traits. From it, one can thus study the history of the lineage.
In Landscape Ecology, it is of particular interest to study the geographical position of the lineage through time. This is because richness patterns are thought to arise through allopatric speciation, where populations get separated in space and time. This is the topic of the following section.
## Projection of the geographic history
Consider a set of individuals which lineage geographical history has been stored in a vector $`x^{(i)}(t)`$. If individual $`i`$ was born at time $`t^*`$ then $`x^{(i)}(t>t^*)`$ represents its geographical position, while $`x^{(i)}(t<t^*)`$ represents its ancestors geographical position.
### Isolation in time
We want to account for the isolation in time and possibly in space of the lineage of a given individual. That is, for how long and how distant stay lineages appart?
#### Setting
Consider a discrete setting, where $`N(t)`$ individuals (or populations) evolve over a set of $`M`$ demes disposed in a linear fashion, such that $`x^i(t) \in \{1,2,...,M\}`$. The population is characterised by the counting measure
```math
\nu = \sum_i^{N(t)} \delta_{x^i(t)}
```
We define the geographical history hamming distance $`\mathfrak{g}`$ as
```math
\mathfrak{h}\big(x^{(i)}(t),x^{(j)}(t)\big) = \int_0^t \text{ceil}(\frac{|x^{(i)} - x^{(j)}|(s)}{M-1})ds
```
We extend this definition with the measure $`\mathfrak{h}^*`$ which takes into account geographic distance
```math
\mathfrak{h}^*\big(x^{(i)}(t),x^{(j)}(t)\big) = \int_0^t \Big[x^{(i)} - x^{(j)} \Big]^2(s)\, ds
```
Finally, we introduce the geotrait distance $`\mathfrak{g}`$ as
```math
\mathfrak{g}\big(x^{(i)}(t),x^{(j)}(t)\big) = \Big[\int_0^t ( x^{(i)} - x^{(j)} )(s) \, ds\Big]^2
```
Note that by the triangle inequality we have that
```math
\mathfrak{g}\big(x^{(i)}(t),x^{(j)}(t)\big) \leq \mathfrak{h}^*\big(x^{(i)}(t),x^{(j)}(t)\big)
```
Equality should arises if $`x^{(i)}, x^{(j)}`$ are positively linearly dependent (which should not be the case).
However, what we are eventually interested is a population measure. This measure should be related to the average pairwise distance across the population. Hence we define $`\mathcal{D_d}(\nu,t)`$ such that
```math
\mathcal{D_d}(\nu,t) = \frac{1}{2N^2}\sum_i^N \sum_j^N d(x^{(i)},x^{(j)},t)
```
Let $`g^{(i)}(t) = \int_0^t x^{i}(t) dt`$ and $`G(t) = \{g^{(i)}(t), i \in \{1,2,\dots,N(t)\}\}`$. By observing the following
```math
\frac{1}{2N^2}\sum_{i,j}^N (y_i-y_j)^2 \\=
\frac{1}{2N^2}\sum_{i,j}^N ((y_i-\bar{y}) -(y_j-\bar{y}))^2 \\
= \frac{1}{2N^2} 2N \sum_i^N(y_i-\bar{y})^2 = \text{Var}(Y)
```
Thus we have $`\mathfrak{g}\big(x^{(i)}(t),x^{(j)}(t)\big) = [g^{(i)}(t) - g^{(j)}(t)]^2`$ and hence
```math
\mathcal{D_\mathfrak{g}}(\nu,t) = \text{Var}(G).
```
Also remark that
```math
\mathcal{D_{h^*}}(\nu,t) = \frac{1}{2N^2}\sum_i^N \sum_j^N h^*(x^{(i)},x^{(j)},t) \\
= \frac{1}{2N^2}\sum_i^N \sum_j^N \int_0^t \Big[x^{(i)} - x^{(j)} \Big]^2(s)\, ds \\
= \int_0^t \frac{1}{2N^2}\sum_i^N \sum_j^N \Big[x^{(i)} - x^{(j)} \Big]^2(s)\, ds \\
= \int_0^t \text{Var}(X)(s)ds
```
One could also imagine a value $`h^{(i)}(t) = \frac{1}{2N}\sum_j^{N(t)}\mathfrak{h}^*\big(x^{(i)}(t),x^{(j)}(t)\big)`$ and in this case we would have
```math
\mathcal{D_\mathfrak{h}}(\nu,t) = \frac{1}{N}\sum_i^{N(t)} h^{(i)}(t)
```
### Mobility
How much do lineages move?
Here is a time average of the speed
```math
\frac{1}{N} \sum_i < \partial_t x^{(i)}(t) >_t \\
= \frac{1}{N} \sum_i \int_0^t \partial_s x^{(i)}(s) ds \\
= \frac{1}{N} \sum_i [x^{(i)}(t) - x^{(i)}(0))]
```
But one could also have a moving average, that is, averaging
```math
\frac{1}{N} \sum_i \sum_j \int_{j\tau}^{j(\tau+1)} \frac{1}{\tau} \partial_s x^{(i)}(s) ds
```
\ No newline at end of file
# Mathematical foundations
- [Dieckmann and Law](https://link.springer.com/article/10.1007%2FBF02409751) (probably more accessible than following), and then
- [Nicolas, Ferriere and Meleard](https://www.sciencedirect.com/science/article/pii/S0040580905001632?via%3Dihub)
showed that the canonical equations of adaptive dynamics, describing the evolution in the phenotypic space, can be derived by considering the stochastic individual-based model corresponding to
```math
\partial_t u(t,x) = u(t,x)(1-\frac{\int \alpha(x,y) u(y,t) dy}{K(x)})
\tag{1}
```
in the limit of **rare mutations, small mutational effects, and infinite population sizes**.
Under these assumptions, Dieckmann and Law showed that adaptive dynamics is the first-order approximation of the mean path averaged over infinitely many realizations of the stochastic simulations obtained from the individual-based model.
[[_TOC_]]
## The basics of Adaptive Dynamics
> This section is inspired from *Adaptive Diversification, Doebeli 2011*
Imagine a monomorophic population with trait $`x`$, which follows the dynamics
```math
\partial_t u(t,x) = u(t,x)\big[b(x) - c(x)u(t,x)\big]
```
> Logistic models are thought to be mathematically representative of a large class of models (During2008) :smirk:. Doebeli "Chaos" paper might also be relevant for this purpose
Equilibrium is given by $`u^*(x) \equiv K(x) = b(x)/c(x)`$. Hence one can reformulate equation above by
```math
\partial_t u(t,x) = u(t,x) \left(b(x) - \frac{b(x)u(t,x)}{K(x)}\right )
```
Adaptive dynamics aims at determining the evolutionary trajectory of the trait $`x`$ by considering the fate of rare mutants with traits $`y`$ in resident monomorphic population with trait $`x`$. The resident population is assumed to be at equlibrium $`K(x)`$. Because mutant is rare, the mutant's population dynamics is only affected by the density of the resident, which in turn is unaffected by the mutant's invasion attempt, and hence remains at $`K(x)`$. Hence the effective density that mutant experiences during invasion attempt is detrmined by $`\alpha(x,y)K(x)`$. The function $`\alpha`$ is the **competition kernel** and describes the strength of competition that exerts the phenotype $`x`$ on phenotype $`y`$.
> For our purpose, we consider it as symmetric, which corresponds to the canonical biological examples of birds with different beak size. We think this is also valid for plants, and in general along a genus.
> It is assumed that $`\alpha(x,x) = 1`$ (scaled to unity).
#### Stabilising selection
Stabilising selection ensures that trait $`x`$ is contained, thus avoiding regions of extreme trait values that would be biologically unrealistic. It can appear by limiting birth rate $`b(x)`$ or increasing death rate $`c(x)`$ for extreme x
Now we can determine the population dynamics of a mutant with population size $`N_{mut}`$ and trait $`y`$ experiencing the effective density $`N_{eff} = \alpha(x,y)K(x)`$
```math
\begin{aligned}
\frac{\partial N_{mut}}{\partial x}(t,x) &= N_{mut}(b(y) - c(y)N_{eff})\\
&= N_{mut}(b(y) - c(y)\, \alpha(x,y)K(x))\\
&= N_{mut}(b(y) - \frac{b(y)\, \alpha(x,y)K(x)}{K(y)})
\end{aligned}
```
## Invasion fitness