EcoEvoModelZoo

This model is inspired from Akesson et al. 2021.

Specialized variants are provided in other AbstractModels. Best constructure with init_akesson_model.

This model is inspired from Boussange & Pellissier. 2022

Arguments

• mp: the model parameters
• g: the spatial graph
• phen_space: the discretized phenotypic space
• birth_fn: the birth function. Should be of the form birth_fn(x, s, p)
• competition_fn: the competition function. Should be of the form competition_fn(u_xs2, x, s1, s2, p)

Mathematical model

∂ₜuˣ(s) = uˣ(s)[birth_fn(x, s, p) - uˣ(s)∫competition_fn(uˣ(s₂), x, s, s₂, p) ds₂] + 1/2 μ σ² Δₛuˣ(s) + m L u(s)

where L is the Laplacian matrix of the spatial graph g.

Example 1

using ParametricModels
using UnPack
using DocStringExtensions
using Statistics
using ComponentArrays
using Distributions
using OrdinaryDiffEq
using UnPack
## Defining phenotypic space and spatial graph
M = 7
g = star_graph(M)

rS = 1f0
dS = 0.02f0
phen_space = collect(range(-rS,rS,step=dS)) #grid

##
# Defining birth and competition functions
soptim = 0.5f0 * [-1,1,-1,1,-1,1,-1]
birth_fn(x, s, p) = max(0f0, 1f0 - (soptim[x] - s)^2)
competition_fn(u_xs2, x, s1, s2, p) = u_xs2 ./ p.K

## Defining parameters
σ_mu = 5f-2;
mu = 0.1f0
m = 0.1
K = 1.

p = ComponentArray(σ_mu = σ_mu,
                    mu = mu,
                    m = m,
                    K = K)


## rest of the simulation
tend = 1000f0
tspan = (0f0,tend)
tsteps = (tspan[1]):1.:(tspan[end])
u0 = vcat([K .* pdf.(Normal(so,σ_mu),phen_space') for so in soptim]...)

mp = ModelParams(;p,
                tspan,
                u0,
                alg=Tsit5(),
                saveat = tsteps)

model = EcoEvoGraph(mp, g, phen_space, birth_fn, competition_fn)

@time sol = simulate(model)

# Plotting results!
using PythonCall; plt = pyimport("matplotlib.pyplot")
fig, ax = plt.subplots(1)
ax.plot(model.phen_space, sol[end][1,:], color = "tab:red", label = "Node 1")
ax.plot(model.phen_space, sol[end][2,:], color = "tab:blue", label = "Node 2")
ax.set_xlabel("Phenotype")
ax.set_ylabel("Population number")
display(fig)

Example 2: trait dependent competition

# Trait dependent compeition
rS = 3f0
dS = 0.02f0
phen_space = collect(range(-rS,rS,step=dS)) #grid

p = ComponentArray(σ_mu = σ_mu,
                    mu = mu,
                    m = m,
                    K = K,
                    σ_α = 0.5)
# trait dependent competition
competition_fn(u_xs2, x, s1, s2, p) = u_xs2 * exp(- 0.5f0 * (s1 - s2)^2 ./ p.σ_α^2) ./ p.K

u0 = vcat([K .* pdf.(Normal(so,σ_mu),phen_space') for so in soptim]...)

mp = ModelParams(;p,
                tspan,
                u0,
                alg=Tsit5(),
                saveat = tsteps)
model = EcoEvoGraph(mp, g, phen_space, birth_fn, competition_fn)

sol = simulate(model)
fig, ax = plt.subplots(1)
ax.set_title("Trait-dependent competition")
ax.plot(model.phen_space, sol[end][1,:], color = "tab:red", label = "Node 1")
ax.plot(model.phen_space, sol[end][2,:], color = "tab:blue", label = "Node 2")
ax.set_xlabel("Phenotype")
ax.set_ylabel("Population number")
display(fig)

EcosystemModelMcKann

This model is inspired from McCann 1994. Similar to EcosystemModelOmnivory, but without monivory.

Model parameters

• x_c, x_p: mass-specific metabolic rate of consumers and predators
• y_c, y_p: ingestion rate per unit metabolic rate of consumers and predators.
• R_0, C_0: half saturation densities for the type II functional responses of the consumers and predators

Example

```julia alg = BS3() abstol = 1e-6 reltol = 1e-6 tspan = (0., 800) tsteps = 550:4:800 ptrue = (xc = [0.4], xp = [0.08], yc = [2.01], yp = [5.00], R0 = [0.16129], C_0 = [0.5])

u0_true = [0.5,0.8,0.5]

model = EcosystemModelMcCann(ModelParams(;p = ptrue, tspan, u0 = u0true, alg, reltol, abstol, saveat = tsteps, verbose = false, # suppresses warnings for maxiters maxiters = 50000, )) sol = simulate(model, u0 = u0true)

EcosystemModelOmnivory

This model is inspired from McCann 1997.

Model parameters

• x_c, x_p: mass-specific metabolic rate of consumers and predators
• y_c, y_pr, y_pc: ingestion rate per unit metabolic rate of consumers and predators.
• R_0, R_02, C_0 half saturation densities for the type II functional responses of the consumers and predators
• ω: omnivory strength

Example

alg = BS3()
abstol = 1e-6
reltol = 1e-6
tspan = (0., 800)
tsteps = 550:4:800

p_true = (x_c = [0.4], 
        x_p = [0.08], 
        y_c = [2.01], 
        y_pr = [2.00], 
        y_pc = [5.0], 
        R_0 = [0.16129], 
        R_02 = [ 0.5], 
        C_0 = [0.5],
        ω =[ 0.4])

u0_true = [0.5,0.8,0.5]


model = EcosystemModelOmnivory(ModelParams(;p = p_true,
                                        tspan,
                                        u0 = u0_true,
                                        alg,
                                        reltol,
                                        abstol,
                                        saveat = tsteps,
                                        verbose = false, # suppresses warnings for maxiters
                                        maxiters = 50_000,
                                        ))
sol = simulate(model, u0 = u0_true)

```

Returns landscape parameters

The ResourceCompetition model is a type of ecological model that simulates competition for resources between different species of plankton. This model is inspired by Huisman et al. 1999 Nature.

Arguments

• mp is a ModelParams object that contains information about the time span, the algorithm to use for numerical integration, and the tolerances to use.
• nN is the number of plankton species.
• nR is the number of resources.
• S is a vector of length nR containing the supply rates of each resource.
• mu: function for the ResourceCompetition model is set to the default function:

μ(R, r, K) = r .* R ./ (K .+ R)

where:

• R is a vector of length nR containing the current concentration of each resource.
• r is a vector of length nN containing the maximum uptake rate of each plankton species.
• K is a vector of length nN containing the half-saturation constant of each plankton species.

This model is inspired from Huisman et al. 1999 Nature., but where Leibig's law is replaced by imperfect substituable resources (smooth minimum). The smooth min function is parametrized by s, which is a trainable parameter.

Arguments

• mp::MP: ModelParams struct containing parameter values for the model.
• nN::NN: Number of plankton species in the model.
• nR::NR: Number of resource species in the model.
• S::SS: Supply rates.

Example

alg = BS3()
abstol = 1e-6
reltol = 1e-6
tspan = (0.,1000.)
step = 2.
nN = 5
nR = 5
r = ones(nN)
m = D = 0.25
S = [6., 10., 14., 4, 9]
K = [0.39 0.34 0.30 0.24 0.23;
    0.22 0.39 0.34 0.30 0.27;
    0.27 0.22 0.39 0.34 0.30;
    0.30 0.24 0.22 0.39 0.34;
    0.34 0.30 0.22 0.20 0.39]'

C = [0.04 0.04 0.07 0.04 0.04;
    0.08 0.08 0.08 0.10 0.08;
    0.10 0.10 0.10 0.10 0.14;
    0.05 0.03 0.03 0.03 0.03;
    0.07 0.09 0.07 0.07 0.07]'
s = -7.

p_true = (r = r, m = [m], D = [D], K = K, C = C, s = Float64[s])

N0 = [0.1 + i / 100 for i in 1:5]
R0 = S
u0 = [N0;R0]
tsteps = tspan[1]:step:tspan[2]

model = ResourceCompetitionSmoothMin(ModelParams(;p = p_true,
                                        tspan,
                                        u0,
                                        alg,
                                        reltol,
                                        abstol,
                                        saveat = tsteps
                                        ), 
                                        nN, nR, S)

# Species 1:5
sol = simulate(model, u0 = u0)

General multitrophic ecosystem model that models the population dynamics of multiple species in an ecosystem.

Arguments

• mp: Model parameters of type ModelParams.
• intinsic_growth_rate: A function that takes in the model parameters p and time t and returns the intrinsic growth rate of the population.
• carrying_capacity: A function that takes in the model parameters p and time t and returns the carrying capacity of the population.
• competition: A function that takes in the state u, model parameters p, and time t and returns the effect of population density on the growth rate.
• resource_conversion_efficiency: A function that takes in the model parameters p and time t and returns the efficiency of converting resources into population growth.
• feeding: A function that takes in the state u, model parameters p, and time t and returns the feeding rates of the population.

Example 1

Simulating a chaotic 5 compartment food web model proposed in Post et al. 2000.

alg = BS3()
abstol = 1e-6
reltol = 1e-6
tspan = (0., 600)
tsteps = (tspan[1]):1.:(tspan[end])
N = 5 # number of compartment

x_c = 0.15, 
x_p = 0.08
y_c = 2.3
y_p = 1.7
R_0 = 0.25
C_0 = 0.5
ω = 0.2

# constructing simple chain foodweb
# Population vector is organised in such a way
# N = [R1, C1, P, R2, C2]
foodweb = SimpleWeightedDiGraph(N)
add_edge!(foodweb, 2, 1, 1.) # C1 to R1
add_edge!(foodweb, 5, 4, 1.) # C2 to R2
add_edge!(foodweb, 3, 2, ω) # P to C1
add_edge!(foodweb, 3, 5, 1-ω) # P to C2

H = sparse(zeros(N,N))
H[2,1] = 1 / (x_c * y_c); H[5,4] = 1 / (x_c * y_c); 
H[3,2] = 1 / (x_p * y_p); H[3,5] = 1 / (x_p * y_p)

q = sparse(zeros(N,N))
q[2,1] = x_c * y_c / R_0; q[5,4] = x_c * y_c / R_0;
q[3,2] = x_p * y_p / C_0; q[3,5] = x_p * y_p / C_0

p = (r = vcat(1., -x_c, -x_p, 1., -x_c,), 
        K = ones(N), 
        A = diagm(vcat(1,0,0,1,0)), 
        ϵ = ones(N), 
        q = q,
        H = H, 
        W = adjacency_matrix(foodweb))

u0_true = rand(N)


model = SimpleEcosystemModel(;mp = ModelParams(;p,
                                        tspan,
                                        u0 = u0_true,
                                        alg,
                                        reltol,
                                        abstol,
                                        saveat = tsteps,
                                        verbose = false, # suppresses warnings for maxiters
                                        maxiters = 50_000,
                                        ))
sol = simulate(model)

Temperature as a function of space (x), time (t), and some climate parameters

funcresp!(F, n, p, _)

Type II functional response Input:

• n: Vector of population densities of all species in a given patch
• Th: Vector of handling times (with dummy values for resource species)
• arate: Vector of attack rates (with dummy values for resource species)
• W: Adjacency matrix of trophic network; W(i,j)=1 if i eats j and 0 otherwise

Output:

• A matrix F(i,j), the feeding rate of consumer i on resource j

generate_network(SR::Int, SC::Int)

return matrix W[i,j], which is nonzero if consumer i eats resource j SR is the number of resource, SC the number of consumer species.

The bipartite network is generated as follows. First, both resources and consumers are labeled consecutively, based on their initial temperature adaptations: resource 1 / consumer 1 are the most cold-adapted, and resource S / consumer S the most warm-adapted. Next, we always put a feeding link between consumer i and resource i. Finally, each consumer is randomly linked to four other resource species. This yields a feeding network where every consumer is connected to five resources altogether.

Note: it seems that the definition of W_ij is that it determines which resource i is eat by consumer j

init_akesson_model(
;
    mp,
    land,
    temp,
    width_growth,
    competition,
    trophic,
    SR,
    SC,
    kwargs...
)

Arguments

• width_growth::WidthGrowth: an instance of the WidthGrowth struct representing the width growth function.
• competition::Competition: an instance of the Competition struct representing the competition function.
• trophic::Tr: an instance of the Trophic struct representing the trophic interaction. Defaults to Trophic{false}() representing a single trophic level.
• SR::Int: the number of species in each patch. Defaults to 50.
• SC::Int: the number of patches. Defaults to 0 (only one patch).
• kwargs: additional keyword arguments to pass to the Base.ODEProblem constructor.

Example

Here we simulate Akesson model in the most simple version with one trophic level, one patch, and where competition is temperature dependent. The code starts by setting the simulation time frame, choosing a differential equation solver algorithm, and defining the tolerance levels for the solver. The model is then initialized with various parameters, including the landscape size, carrying capacity, and temperature dependency of competition. The temperature over a period of 6500 years is plotted using matplotlib. The model is then simulated before and after climate change, with the population density and population mean trait (optimal temperature) plotted for each time step.

using EcoEvoModelZoo
using PythonCall
using OrdinaryDiffEq
using ParametricModels
plt = pyimport("matplotlib.pyplot")

tstart = -4000 # starting time (relative to start of climate change at t = 0)
tend = 2500 # time at which integration ends
alg = Tsit5()
abstol = 1e-6
reltol = 1e-6
# integrate prob with Julia
tspan = (tstart, 0)
tsteps = tspan[1]:200:tspan[2]
L = 10
SR = 50
temp = Temperature()
land = Landscape(L)

model = init_akesson_model(;mp = ModelParams(;tspan, alg, reltol, abstol, saveat = tsteps),
                            land = Landscape(L),
                            temp,
                            SR = SR,
                            width_growth = WidthGrowth{:TraitDep}(), 
                            competition = Competition{:TraitDep}(), 
                            trophic= Trophic{false}(), 
                            )

# Plotting temperature over 6500 years, from -4000 to 2500
fig, ax = plt.subplots(1)
ts = -100:1:400
p = get_p(model)
temp_ts = [temp.(land.x, t)[1] for t in ts]
ax.plot(ts, temp_ts)
display(fig)

#################
### BEFORE Climate Change ###
#################
sol = simulate(model)
ode_data = Array(sol)

# initial temperature
temp.(land.x, 100.)

# plotting evolution before CC
fig, ax = plt.subplots(1)
l = 1 # patch index
plotting_N_through_time(ax, ode_data, tsteps, l, model)
# ax.set_yscale("log")
display(fig)


#################
### After Climate Change  ###
#################
tspan = (0, 300)
step = 1
tsteps = tspan[1]:1.:tspan[2]
sol = simulate(model, tspan = tspan, saveat = tsteps, u0 = ode_data[:,:,end])
ode_data = Array(sol)

# plotting evolution after CC
fig, axs = plt.subplots(1, 2)
ts = 1:length(tsteps)
l = 1 # patch index
plotting_N_through_time(axs[0], ode_data, tsteps, l, model)
axs[0].set_ylabel("Population density")
# axs[0].set_yscale("log")
plotting_mu_through_time(axs[1], ode_data, tsteps, l, model)
axs[1].set_ylabel("Population mean trait
(Optimal temperature)")
fig.tight_layout()

Apply twice continuously differentiable smoothed step function to a number x Input: x: Distance from pole, measured in units of the pole-to-equator distance Output: 0 if x < 0; 10x^3-15x^4+6*x^5 if 0 <= x <= 1; otherwise 1