using Petri
using MacroTools
import MacroTools: postwalk, striplines
using ModelingToolkit
import ModelingToolkit: Constant
import Base: ==, ∈, show, ^
using Catlab.Doctrines
import Catlab.Doctrines: ⊗, compose, otimes
using Catlab.WiringDiagrams
using Catlab.Graphics
using Catlab.Graphics.Graphviz
import Catlab.Graphics.Graphviz: Graph, Edge
using SemanticModels
using SemanticModels.ModelTools.WiringDiagrams
using SemanticModels.ModelTools.PetriModels
using LabelledArrays
using OrdinaryDiffEq
import OrdinaryDiffEq: solve
using Plots
Petri.N(s) = 1 #sum(s[1:end-3])
using SemanticModels.ModelTools.OpenPetris
import Petri: NullPetri
import SemanticModels.ModelTools.OpenModels: OpenModel
import SemanticModels.ModelTools: model
⊗(f::OpenModel,g::OpenModel) = otimes(f,g)
⊚(f::OpenModel,g::OpenModel) = compose(f,g)
^(x::Doctrines.FreeSymmetricMonoidalCategory.Ob{:generator}, n::Int) = foldl(⊗, x for i in 1:n)
^(f::Doctrines.FreeSymmetricMonoidalCategory.Hom{:id}, n::Int64) = foldl(⊗, f for i in 1:n)
conj(a::OpenModel, b::OpenModel) = compose(b,a,b)
conj(a, b) = compose(b,a,b)
""" op(f::Model)
return the opposite model you get by reversing the direction of all the transitions
"""
op(f::Model) = Model(f.S, map(f.Δ) do t reverse(t) end)
op(f::OpenModel) = OpenModel(f.dom, op(f.model), f.codom)
wd(ex) = to_wiring_diagram(ex)
function compile(expr::Doctrines.FreeSymmetricMonoidalCategory.Hom{:compose}, ctx::Dict)
# println("∘")
compose(compile.(expr.args, (ctx,)))
end
function compile(expr::Doctrines.FreeSymmetricMonoidalCategory.Hom{:otimes}, ctx::Dict)
# println("⊗")
otimes(compile.(expr.args, (ctx,)))
end
function compile(expr, ctx::Dict)
# @show expr
return ctx[expr]
end
function show(io::IO, z::Petri.Model)
X, Y = z.S, z.Δ
compact = get(io, :compact, true)
if compact
x,y = length(X), length(Y)
print(io,"Model(∣S∣=$x,∣Δ∣=$y)")
else
print(io,"Model(S=$X, Δ=$Y)")
end
end
function show(io::IO, z::OpenModel)
X,Y = z.dom, z.codom
compact = get(io, :compact, true)
if compact
x,y = length(X), length(Y)
print(io,"OpenModel:$x→$y with ")
show(io,z.model)
else
print(io,"Domain: $X\nCodomain: $Y\nModel: ")
show(io, z.model)
end
end
function debug_show(io::IO, f::Model)
X,Y = f.S, f.Δ
print(io,"Model(S=$X, Δ=")
if length(Y) == 0
print("[])")
else
for (i,t) in enumerate(Y)
if i == 1
print(io, "[")
print(io, "$(t[1])→$(t[2]),")
end
if 1 < i < length(Y)
print(io, " $(t[1])→$(t[2]),")
end
if i == length(Y)
print(io, " $(t[1])→$(t[2])])\n")
end
end
end
end
varlookup(x) = begin
string(OpenPetris.X[x].op.name)
end
debug_show(io::IO, f::OpenModel) = begin
X,Y = f.dom, f.codom
print(io,"Domain: $X\nCodomain: $Y\nModel: ")
debug_show(io, f.model)
println(io, "\n")
f′ = OpenModel(f.dom, #varlookup.(f.dom),
Petri.Model(varlookup.(f.model.S),
f.model.Δ), f.codom)#varlookup.(f.codom))
g = SemanticModels.ModelTools.OpenModels.Graph(f′)
end
debug_show(f::OpenModel) = begin
debug_show(stdout, f)
end
function debug_show(io::IO, f::OpenModel, fname::String)
g = debug_show(io, f)
output = run_graphviz(g, prog="dot", format="svg")
write(fname, output)
return g
end
debug_show(f::OpenModel, fname::String) = debug_show(stdout, f, fname)
odeexpr(m::OpenModel) = symbolic_symplify(Petri.odefunc(m.model, :state)) |> striplines
Base.show(io::IO, mime::MIME"text/html", f::OpenModel) = begin
g = debug_show(io, f)
show(io, mime, g)
g
end
import Base.display
display(f::OpenModel) = debug_show(f)
display(f::WiringDiagram) = to_graphviz(f, direction=:horizontal)
pad(f::WiringDiagram, v::Vector) = vcat(v, ["" for i in 1:length(wires(f))-length(v)])
display(f::WiringDiagram, labels::Vector) = SemanticModels.ModelTools.WiringDiagrams.label!(display(f), string.(pad(f, labels)))
Base.Filesystem.mkpath("img")
X = OpenPetris.X;
Agenda
# | Layer | Type | Examples |
---|---|---|---|
1 | Structured | Petri Nets | Epidemiological Models (SIR, SIIR) |
2 | Abstract | Wiring Diagrams | Ecological Models |
3 | Executable | ODEs | Initial Value Problems |
We start by defining some primitive models. First is the spontaneous reaction. It takes $X_1\rightarrow X_2$
OpenModel([1,2], model(PetriModel, Petri.Model([1,2], [(X[1], X[2])])), [1,2])
using SemanticModels.ModelTools.CategoryTheory
OpenPetri([1,2], Petri.Model([1,2], [(X[1],X[2])]), [1,2])
f2 = FinSetMorph(1:2, [1,2])
Decorated(f2, Petri.Model([1,2], [(X[1],X[2])]))
println("\nSpontaneous reaction: spontaneous = X₁→X₂")
spontaneous = OpenModel([1,2],
Petri.Model([1,2], [(X[1],X[2])]),
[1,2])
display(spontaneous)
We can also combine models by stacking them one on top of the other independently $parallel = spontaneous \otimes spontaneous = X_1\rightarrow X_2, X_3\rightarrow X_4$
println("\nParallel reaction: parallel = spontaneous ⊗ spontaneous = X₁→X₂, X₃→X₄")
parallel = OpenPetris.otimes(spontaneous, spontaneous)
display(parallel)
We can compose models by wiring the inputs of one to the outputs of another. We define most of our primitives to have all variables as inputs and outputs.
display(compose(spontaneous,spontaneous))
In order to introduce new variables, you can use the $I$ = eye(1)
model with represents the identity process on 1 variable. This pattern of $a\otimes I; I \otimes b$ will be very important. We will use the $f;g$ operator to denote composition of $f$ followed by $g$ because the mathematical symbol $f\circ g$ is already defined as $g(f)$ which is backwards.
display(eye(1))
display(compose(spontaneous ⊗ eye(1), eye(1) ⊗ spontaneous))
Any Petri Net can be converted into a primitive building block by designating a set of variables to be the input (domain) and output (codomain).
println("\nInfection reaction infect = X₁+X₂→ 2X₂")
infect = OpenModel([1,2], Model([1,2], [(X[1]+X[2], 2*X[2])]), [1,2])
display(infect)
You can reverse the inputs of a model using a special model called $\sigma$. Since all permutations can be factored into swaps (see bubble sort), we can use $\sigma$ to make arbitrary permutations.
σ() = OpenModel([1,2], NullPetri(2), [2,1])
display(σ())
spontaneousdag = compose(σ(), spontaneous)
display(spontaneousdag)
We now have everything we need to make an SIS model. The odeexpr
function converts a petri net into a program for evaluating the corresponding differential equations.
sis = compose(infect, spontaneousdag)
@show odeexpr(sis)
display(sis)
sis.model.S
sis.model.Δ
Petri.fluxes(sis.model)
The SIR model can be built from these primitives too. We use the $a\otimes I; I \otimes b$ pattern to introduce a third variable $X_3=R$
sir = compose(infect⊗eye(1), eye(1)⊗spontaneous)
@show odeexpr(sir)
display(sir)
Every model you build is a new term that you can use to make more complex models. For example $sir\otimes sir$ is two uncoupled $SIR$ diseases.
display(sir⊗sir)
But real world systems don't combine independently all too often, we need a way to couple things. The family of models $m_{i,j}$ combine the $i$th and $j$th input states into one input state, by taking the disjoint union of their transitions. We can use these models to combine systems by "state sharing". This will allow us to build a SIIR model where two strains of a disease compete for susceptible people.
m₁₂ = OpenModel([1], NullPetri(1), [1,1])
println("\nParallel Infections reactions infect ⊗ infect = X₁+X₂→ 2X₂ && X₃ +X₄ → 2X₄")
# set the codomains to be narrower to help with composition
infect′ = OpenModel([1], infect.model, [2])
spontaneous′ = OpenModel([1], spontaneous.model, [2])
parinfect′ = infect′ ⊗ infect′
display(parinfect′)
display(m₁₂)
display(compose(m₁₂, parinfect′))
mutrec = OpenModel([1,2], Model([1,2,3], [(X[1],X[3]), (X[2], X[3])]), [3])
display(mutrec)
siir = compose(m₁₂, parinfect′, mutrec)
display(siir)
odeexpr(siir)
An interesting lemma is that if you have an OpenModel $f$ and a pair of states $x,y$ with no common transition, such that $x\notin dom(f)$ and $y \notin codom(f)$ there is no OpenModel $g$ such that $f;g$ (or $g;f$) has a transition between $x,y$.
This is why we prefer to define primitives such that $dom(f) = codom(f) = States(f)$ when convenient.
The Petri Net layer is useful because it contains a mathematically rigorous formulation of the scientific system under study that can be used to generate simulations at the executable layer. However, they do not cont any capacity for abstraction or encapsulation, which are necessary for human-level intelligent reasoning.
We now introduce Wiring Diagrams through some motivating examples from Ecology.
println("\nCreating food web processes birth, death, predation")
birth = Model([1], [(X[1], 2X[1])])
death = Model([1, 2], [(X[1], X[2])])
pred(α,β,γ) = Model([1, 2], [(α*X[1] + β*X[2], γ*X[2])])
b = OpenModel([1], birth, [1])
d = OpenModel([1], death, [1])
p(α, β, γ) = OpenModel([1,2], pred(α, β, γ), [1,2])
println("\nCreating food web processes σ, predation†")
pdag(α,β,γ) = OpenModel([1,2], Model([1, 2], [(α*X[2] + β*X[1], γ*X[1])]), [1,2])
# Catlab expressions for our variables
Xob = Ob(FreeSymmetricMonoidalCategory, :X)
idₓ = id(Xob)
σh = braid(Xob, Xob)
bh = Hom(:birth, Xob,Xob)
dh = Hom(:death, Xob, Xob)
ph = Hom(:predation, Xob^2, Xob^2)
pdagh = Hom(Symbol("p⋆"), Xob^2, Xob^2)
lookuptable = Dict( idₓ=>eye(1), σh=>σ(), bh=>b, dh=>d, ph=>p(1,1,1.15))
set_primitive!(ctx::Dict, hom, f::OpenModel) = begin
ctx[hom] = f
return ctx
end
set_primitive(hom, f::OpenModel) = set_primitive!(lookuptable, hom, f)
compile(x) = compile(x, lookuptable)
A simple model of predator-prey dynamics that is widely used in ecology is the Lotka-Volterra model. The prey species grows exponentially, and is eaten by a predator species that suffers from contant probability of death (exponential decay). This model has periodic solutions that are non-trigonometric.
In the following diagram of the Lotka-Volterra System, wires represent species and boxes represent interactions. The wires that go into a box are the domain and the wires that go out are the codom. You read it from left to right.
bd = b⊗d
lvh = compose(bh⊗idₓ, ph, idₓ⊗dh)
drawhom(lvh, "img/lv_wd")
display(wd(lvh), [:Wolf, :Sheep, :Sheep, :Wolf, :Wolf,:Sheep])
The wiring diagram approach allows you to specify a model by connecting up a circuit of component models. We can convert the diagramatic representation into a Petri Net with compile
. Note that the state $X_3$ is neither an input nor an output to the model. This implies that no others species can interact with the deceased wolves (ie. vultures are forbidden).
lvh
lv = compile(lvh)
println(lv)
debug_show(lv,"img/lv.svg")
At this point we can see all three layers of the Lotka-Volterra model. Wiring Diagram -> Petri Net -> ODEs
@show odeexpr(lv)
f = Petri.mk_function(lv.model)
νS = 4
ηS = 1.0
δW = 2.0
params = [νS, ηS, δW]
u0 = @LArray [40,2,0.0] (:X₁, :X₂, :X₃)
prob = ODEProblem(f,u0,(0.0,8.0),params)
sol = solve(prob,Tsit5(), tol=1e-5)
@show sol[end]
varnames = [:Sheep, :Wolves, :Dwolves]
plt = plot(sol, vars =1:length(varnames), labels=varnames)
savefig(plt, "img/lv_sol.pdf")
plt
The Lotka-Volterra Model is a classic example of dynamical systems in mathematical ecology. We can view it at three layers of knowledge.
bdd = bd⊗d
bdb = b⊗(d⊗b)
println("bipredation is (p⊗I)⋅(I⊗p)")
bipredation = compose(p(1,1,2)⊗eye(1), eye(1)⊗p(1,1,2))
println("\nfoodchain is (bipredation)⋅(bdd). A fish, a bigger fish, and biggest fish")
foodchainh = compose(ph⊗idₓ, idₓ⊗ph, bh⊗dh⊗dh)
homx = SemanticModels.ModelTools.WiringDiagrams.canonical(FreeSymmetricMonoidalCategory, foodchainh)
println("Cannonical form construction proves: $foodchainh == $homx")
#drawhom(foodchainh, "img/foodchain_wd")
display(wd(foodchainh), [:fish, :Fish, :Shark])
foodchain = compile(foodchainh)
display(foodchain)
println("As an ordinary differential equation:")
odeexpr(foodchain).args[end].args[end].args[end]
foodstarh = compose(ph⊗idₓ, idₓ⊗compose(σh,ph, σh), bh⊗dh⊗bh)
display(wd(foodstarh), [:Rabbit, :Wolf, :Sheep])
foodstar = compile(foodstarh)
println(foodstar)
homx = canonical(FreeSymmetricMonoidalCategory, foodstarh)
println("Cannonical form construction proves: $foodstarh == $homx")
println("As an ordinary differential equation:")
@show odeexpr(foodstar)
display(foodstar)
dualinfect = compose(ph⊗idₓ, idₓ⊗pdagh)
cur = Hom(:cure, Xob⊗Xob, Xob⊗Xob)
curdag = Hom(Symbol("cur⋆"), Xob^2, Xob^2)
inf = Hom(Symbol("inf¹³₂₃"), Xob^3, Xob^3)
#inf′ = Hom(Symbol("inf³¹₂₁"), Xob^3, Xob^3)
#malariah = compose(inf⊗idₓ, idₓ⊗compose(compose(σh⊗idₓ, idₓ⊗σh), inf, compose(idₓ⊗idₓ⊗idₓ, idₓ⊗σh, σh⊗idₓ)),cur⊗conj(cur, σh))
malariah = compose(inf⊗idₓ, idₓ⊗compose(σh⊗idₓ, idₓ⊗σh, σh⊗idₓ), idₓ⊗inf, idₓ⊗conj(σh⊗idₓ, idₓ⊗σh), cur⊗compose(σh, cur, σh))
#drawhom(malariah, "img/malaria_wd")
display(wd(malariah), ["Sp", "Ip", "Im", "Sm"])
cure = conj(spontaneous, σ())
display(cure)
curedag = op(cure)
curedag == spontaneous
In order to build the mosquito borne illness model, we need a new primative. In this case $X_1 + X_3 \rightarrow X_2 + X_3$
trinary = OpenModel([1,2,3], Model([1,2,3], [(X[1]+X[3], X[2]+X[3])]), [1,2,3])
set_primitive(inf, trinary)
display(trinary)
inf′ = compose(σh⊗idₓ, idₓ⊗σh, σh⊗idₓ, inf, conj(σh⊗idₓ, idₓ⊗σh))
display(wd(inf′))
dualinfecth = compose(inf⊗idₓ, idₓ⊗inf′)
display(wd(dualinfecth), ["Sp", "Ip","Im","Sm"])
dualinfect = compile(dualinfecth)
display(dualinfect)
set_primitive(ph, p(1,1,1.15))
set_primitive(cur, cure)
set_primitive(inf, trinary)
#lookuptable = Dict(bh=>b, dh=>d, σh=>σ(), ph=>p(1,1,1.15), inf=>trinary,idₓ=>eye(1), cur=>cure)
malaria = compile(malariah)
homx = canonical(FreeSymmetricMonoidalCategory, malariah)
println("Cannonical form construction proves:\n\t $malariah \n\t== \n\t $homx")
println("As an ordinary differential equation:")
@show odeexpr(malaria)
debug_show(malaria, "img/malaria.svg")
f = Petri.mk_function(malaria.model)
u0 = @LArray [20,1,0,100.0] (:X₁, :X₂, :X₃, :X₄)
varnames = ["\$S_p\$", "\$I_p\$", "\$I_m\$", "\$S_m\$"];
Because we assumed that Malaria is SIS for both People and Mosquitos the only solution is endemic malaria, where there is a constant fraction of infected people at equilibrium.
βpm = 0.0045
βmp = 0.0045
ρp = 0.5
ρm = 0.04
νp = 0.006
νm = 0.1
params = [βpm, βmp,ρp,ρm]
prob = ODEProblem(f,u0,(0.0,365.0),params)
sol = solve(prob,Tsit5(), tol=1e-5)
@show sol[end]
varnames = [:Sp, :Ip, :Im, :Sm, :Dp, :Dm]
plt = plot(sol, vars =1:length(varnames)-2, labels=varnames)
savefig(plt, "img/malaria_sol_endemic.png")
plt
In the case where malaria is fatal to some fraction of the population and you start off with more mosquitos than people:
Blessed are the mosquitos, for they shall inherit the earth
malvitalh = compose(malariah, bh⊗dh⊗dh⊗bh)
display(wd(malvitalh), ["Sp", "Ip","Im","Sm"])
malvital = compile(malvitalh); #compose(malaria, b⊗d⊗d⊗b)
display(malvital)
println("As an ordinary differential equation:")
@show fex = odeexpr(malvital)
f = Petri.mk_function(malvital.model)
u0 = @LArray [20,1,0,100,0.1,0] (:X₁, :X₂, :X₃, :X₄, :X₅, :X₆)
varnames = ["\$S_p\$", "\$I_p\$", "\$I_m\$", "\$S_m\$", "\$Dp\$", "\$Dm\$"];
βpm = 0.045
βmp = 0.045
ρp = 0.1
ρm = 0.1
νp = 0.1
δp = 0.1
δm = 0.1
νm = 0.1
params = [βpm, βmp,ρp,ρm, νp, δp, δm, νm]
prob = ODEProblem(f,u0,(0.0,45.0),params)
sol = solve(prob,Tsit5(), tol=1e-5)
@show sol[end]
varnames = [:Sp, :Ip, :Im, :Sm, :Dp, :Dm]
plt = plot(sol, vars =1:length(varnames)-2, labels=varnames)
# ylims!(plt, (0,100))
# xlims!(plt, (0,50))
# yaxis!("amount", :log10)
savefig(plt, "img/malaria_sol_fatal.png")
plt
u0 = @LArray [20,1,0,100,0.1,0] (:X₁, :X₂, :X₃, :X₄, :X₅, :X₆)
varnames = ["\$S_p\$", "\$I_p\$", "\$I_m\$", "\$S_m\$", "\$Dp\$", "\$Dm\$"];
βpm = 0.024
βmp = 0.024
ρp = 0.3
ρm = 0.3
νp = 0.4
δp = 0.051
δm = 0.02025
νm = 0.1
params = [βpm, βmp,ρp,ρm, νp, δp, δm, νm]
prob = ODEProblem(f,u0,(0.0,150.0),params)
sol = solve(prob,Tsit5(), tol=1e-5)
@show sol[end]
varnames = [:Sp, :Ip, :Im, :Sm, :Dp, :Dm]
plt = plot(sol, vars =1:length(varnames)-2, labels=varnames)
# ylims!(plt, (0,100))
# xlims!(plt, (0,50))
# yaxis!("amount", :log10)
savefig(plt, "img/malaria_sol_fatal.png")
plt
We want to perform cross domain model fusion, by incorporated ecology into our epidemiology. This means adding a species of birds that eat the mosquitos.
This gives rise to the following high level model structure. Where malaria, mutual predation,
and vitals
are submodels
display(wd(compose(Hom(:malaria, Xob^4, Xob^4)⊗idₓ, idₓ^2⊗Hom(:mutualpredation, Xob^3, Xob^3), Hom(:vitals, Xob^5, Xob^5))), [:B, :Sp, :Ip, :Im, :Sm])
birdsh = compose(σh⊗idₓ,idₓ⊗ph, σh⊗idₓ, idₓ⊗ph)
display(wd(birdsh), [:Im, :B, :Sm])
drawhom(birdsh, "img/birds_wd")
homx = canonical(FreeSymmetricMonoidalCategory, birdsh)
println("Cannonical form construction proves: $birdsh == $homx")
# vitals for Sp, Ip, Im, Sm, B are:
# born, die (of malaria), die (of malaria), born, die (starvation)
vitalsh = bh⊗dh⊗dh⊗bh⊗dh
#birdmalh = compose(malvitalh⊗idₓ, idₓ^2⊗compose(birdsh, idₓ^2⊗dh))
birdmalh = compose(malariah⊗ idₓ, idₓ^2⊗birdsh, vitalsh)
drawhom(birdmalh, "img/birdmal_wd")
display(wd(birdmalh), [:B, :Sp, :Ip, :Im, :Sm, :Sm, :Ip,:B,:Sm,:Sm, :B,:Im,:Ip, :Sp,:Ip,:Sm,:Im, :Ip,:Sp,:Im,:Sm,:Im, :B, :Sp])
birds = compose(σ()⊗eye(1), eye(1)⊗p(1,1,1.15), σ()⊗eye(1), eye(1)⊗p(1,1,1.15))
debug_show(birds, "img/birds.svg")
birdmal = compile(birdmalh)
debug_show(birdmal, "img/birdmal.svg");
The coupled model of Malaria and Predation:
println("As an ordinary differential equation:")
fex = @show odeexpr(birdmal)
Petri.N(s) = 1 #sum(s[1:end-3])
f = Petri.mk_function(birdmal.model)
u0 = @LArray [20,1,0,100,0.2,0,0,0.0] (:X₁, :X₂, :X₃, :X₄, :X₅, :X₆, :X₇, :X₈)
varnames = ["\$S_p\$", "\$I_p\$", "\$I_m\$", "\$S_m\$", "\$B\$", "\$Dp\$", "\$Dm\$", "\$D_b\$"];
In the following setting everyone dies, but then the mosquito population becomes locked in a Lotka-Volterra struggle for the earth.
βpm = 0.045
βmp = 0.045
ρp = 0.5
ρm = 0.04
ηI = 0.12
ηS = 0.12
νp = 0.006
δp = 1.0
δm = 0.03
νm = 0.1
δb = 0.85
params = [βpm, βmp,ρp,ρm,ηI, ηS, νp, δp, δm, νm, δb]
prob = ODEProblem(f,u0,(0.0, 160.0),params)
sol = solve(prob,Tsit5(), tol=1e-5)
@show sol[end]
varnames = [:Sp, :Ip, :Im, :Sm, :B, :Dp, :Dm,:Db]
plotvars = [1,2,3,4,5]
plt = plot(sol, vars =plotvars, labels=varnames[plotvars])
ylims!(plt, (0,100))
xlims!(plt, (0,10))
# yaxis!("amount", :log10)
savefig(plt, "img/birdmal_sol_subcritical.png")
plt
plt = plot(sol, vars =1:length(varnames)-3, labels=varnames)
#yaxis!(plt, (0,30))
In the following example, malaria is endemic and cyclical within a population. After a phase of initially wild fluctuation, we settle down into a rythm of cyclic epidemics. The exponential growth and decay cycles are driven by the Lotka Volterra dynamics between the birds and the mosquitos.
βpm = 0.045
βmp = 0.045
ρp = 1.2
ρm = 0.4
ηI = 0.05
ηS = 0.05
νp = 0.600
δp = 1.9
δm = 0.03
νm = 0.6
δb = 0.85
params = [βpm, βmp,ρp,ρm,ηI, ηS, νp, δp, δm, νm, δb]
prob = ODEProblem(f,u0,(0.0,375.0),params)
sol = solve(prob,Tsit5(), tol=1e-5)
@show sol[end]
plt = plot(sol, vars =1:length(varnames)-3, labels=varnames)
ylims!(plt, (0,250))
xlims!(plt, (0,175))
savefig(plt, "img/birdmal_sol_multiepi.pdf")
plt
Zooming in on a cycle you can see that the mosquitos and birds interact with a mixture of Lotka Volterra dynamics and SIS dynamics. The population of birds tracks the population of mosquitos, as the number of birds increases, they drive mosquitos to near extinction. Once there are no birds to eat the mosquitos, the population of mosquitos explodes. At this point the human population has rebounded and is decimated by an epdidemic of malaria. The epidemic ends when the number of human hosts is small enough to all for recovery of the infected mosquitos. Recall that we assumed that Malaria is (SISD) in both people and mosquitos.
You can see from the zoomed in plot that the malaria outbreaks occur approximately half a cycle shifted from the Lotka Volterra cycle in the mosquito population.
plt = plot(sol, vars =1:length(varnames)-3, labels=varnames)
ylims!(plt, (0,175))
xlims!(plt, (130,175))
savefig(plt, "img/birdmal_sol_multiepi_zoom.pdf")
plt
We can even generate the same model as a stochastic simulation. Allow you to describe the model in a declarative high level language and then compile multiple implementations of that model depending on the scientific question at the time.
striplines.(Petri.Δ(birdmal.model, :u))
The following choice of parameters exhibits an unstable dynamics, without approaching a periodic equilibrium.
βpm = 0.045
βmp = 0.045
ρp = 0.5
ρm = 0.05
ηI = 0.6
ηS = 0.12
νp = 0.04
δp = 1.9
δm = 0.03
νm = 0.1
δb = 0.85
params = [βpm, βmp,ρp,ρm,ηI, ηS, νp, δp, δm, νm, δb]
prob = ODEProblem(f,u0,(0.0,560.0),params)
sol = solve(prob,Tsit5(), tol=1e-5)
@show sol[end]
varnames = [:Sp, :Ip, :Im, :Sm, :B, :Dp, :Dm, :Db]
plt = plot(sol, vars =1:length(varnames)-3, labels=varnames)
#ylims!(plt, (0,90))
#xlims!(plt, (0,50))
# yaxis!("amount", :log10)
savefig(plt, "img/birdmal_sol_subcritical.png")
plt
The outbreaks exhibit irregular behavior
plt = plot(sol, vars =[2,3], labels=varnames)
#ylims!(plt, (0,10))
#xlims!(plt, (0,450))