In [1]:
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
using SemanticModels
using SemanticModels.ModelTools.WiringDiagrams
┌ Info: Recompiling stale cache file /Users/jfairbanks6/.julia/compiled/v1.0/Petri/FiiRl.ji for Petri [4259d249-1051-49fa-8328-3f8ab9391c33]
└ @ Base loading.jl:1190
WARNING: using Doctrines.Δ in module Petri conflicts with an existing identifier.
┌ Info: Recompiling stale cache file /Users/jfairbanks6/.julia/compiled/v1.0/SemanticModels/Di9ju.ji for SemanticModels [f5ac2a72-33c7-5caf-b863-f02fefdcf428]
└ @ Base loading.jl:1190
┌ Warning: Error requiring Unitful from DiffEqBase:
│ UndefVarError: AbstractQuantity not defined
│ Stacktrace:
│  [1] getproperty(::Module, ::Symbol) at ./sysimg.jl:13
│  [2] top-level scope at /Users/jfairbanks6/.julia/packages/DiffEqBase/Fs7ts/src/init.jl:90
│  [3] eval at ./boot.jl:319 [inlined]
│  [4] eval at /Users/jfairbanks6/.julia/packages/DiffEqBase/Fs7ts/src/DiffEqBase.jl:1 [inlined]
│  [5] (::getfield(DiffEqBase, Symbol("##453#477")))() at /Users/jfairbanks6/.julia/packages/Requires/9Jse/src/require.jl:67
│  [6] err(::getfield(DiffEqBase, Symbol("##453#477")), ::Module, ::String) at /Users/jfairbanks6/.julia/packages/Requires/9Jse/src/require.jl:38
│  [7] #452 at /Users/jfairbanks6/.julia/packages/Requires/9Jse/src/require.jl:66 [inlined]
│  [8] withpath(::getfield(DiffEqBase, Symbol("##452#476")), ::String) at /Users/jfairbanks6/.julia/packages/Requires/9Jse/src/require.jl:28
│  [9] (::getfield(DiffEqBase, Symbol("##451#475")))() at /Users/jfairbanks6/.julia/packages/Requires/9Jse/src/require.jl:65
│  [10] #invokelatest#1 at ./essentials.jl:697 [inlined]
│  [11] invokelatest at ./essentials.jl:696 [inlined]
│  [12] #3 at /Users/jfairbanks6/.julia/packages/Requires/9Jse/src/require.jl:19 [inlined]
│  [13] iterate at ./generator.jl:47 [inlined]
│  [14] _collect(::Array{Function,1}, ::Base.Generator{Array{Function,1},getfield(Requires, Symbol("##3#4"))}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at ./array.jl:632
│  [15] map at ./array.jl:561 [inlined]
│  [16] loadpkg(::Base.PkgId) at /Users/jfairbanks6/.julia/packages/Requires/9Jse/src/require.jl:19
│  [17] #invokelatest#1 at ./essentials.jl:697 [inlined]
│  [18] invokelatest at ./essentials.jl:696 [inlined]
│  [19] _tryrequire_from_serialized(::Base.PkgId, ::UInt64, ::Nothing) at ./loading.jl:651
│  [20] _require_from_serialized(::String) at ./loading.jl:679
│  [21] _require(::Base.PkgId) at ./loading.jl:967
│  [22] require(::Base.PkgId) at ./loading.jl:858
│  [23] require(::Module, ::Symbol) at ./loading.jl:853
│  [24] top-level scope at In[1]:13
│  [25] eval at ./boot.jl:319 [inlined]
│  [26] softscope_include_string(::Module, ::String, ::String) at /Users/jfairbanks6/.julia/packages/SoftGlobalScope/cSbw5/src/SoftGlobalScope.jl:218
│  [27] execute_request(::ZMQ.Socket, ::IJulia.Msg) at /Users/jfairbanks6/.julia/packages/IJulia/fRegO/src/execute_request.jl:67
│  [28] #invokelatest#1 at ./essentials.jl:697 [inlined]
│  [29] invokelatest at ./essentials.jl:696 [inlined]
│  [30] eventloop(::ZMQ.Socket) at /Users/jfairbanks6/.julia/packages/IJulia/fRegO/src/eventloop.jl:8
│  [31] (::getfield(IJulia, Symbol("##15#18")))() at ./task.jl:259
└ @ Requires /Users/jfairbanks6/.julia/packages/Requires/9Jse/src/require.jl:40
In [2]:
using LabelledArrays
using OrdinaryDiffEq
import OrdinaryDiffEq: solve
using Plots
Petri.N(s) = 1 #sum(s[1:end-3])
In [3]:
(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
Out[3]:
compile (generic function with 3 methods)
In [4]:
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

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")
    g = 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
Out[4]:
odeexpr (generic function with 1 method)
In [5]:
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)))
Out[5]:
display (generic function with 33 methods)
In [6]:
Base.Filesystem.mkpath("img")
X = Petri.X;

Semantic Modeling Software Demo

Agenda

# Layer Type Examples
1 Structured Petri Nets Epidemiological Models (SIR, SIIR)
2 Abstract Wiring Diagrams Ecological Models
3 Executable ODEs Initial Value Problems

Primitives of Petri Net Modeling Framework

We start by defining some primitive models. First is the spontaneous reaction. It takes $X_1\rightarrow X_2$

In [7]:
println("\nSpontaneous reaction: spontaneous = X₁→X₂")
spontaneous = OpenModel([1,2], Model([1,2], [(X[1],X[2])]), [1,2])
display(spontaneous)
Spontaneous reaction: spontaneous = X₁→X₂
Domain: [1, 2]
Codomain: [1, 2]
Model: Model(S=[1, 2], Δ=[X₁()→X₂(), X₁()→X₂()])


Out[7]:
G X1 X1 T1 T1 X1->T1 O1 O1 X1->O1 X2 X2 O2 O2 X2->O2 T1->X2 I1 I1 I1->X1 I2 I2 I2->X2

Parallel reaction:

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$

In [8]:
println("\nParallel reaction: parallel = spontaneous ⊗ spontaneous = X₁→X₂, X₃→X₄")
parallel = spontaneous  spontaneous
display(parallel)
Parallel reaction: parallel = spontaneous ⊗ spontaneous = X₁→X₂, X₃→X₄
Domain: [1, 2, 3, 4]
Codomain: [1, 2, 3, 4]
Model: Model(S=[1, 2, 3, 4], Δ=[X₁()→X₂(), X₃()→X₄()])


Out[8]:
G X1 X1 T1 T1 X1->T1 O1 O1 X1->O1 X2 X2 O2 O2 X2->O2 X3 X3 T2 T2 X3->T2 O3 O3 X3->O3 X4 X4 O4 O4 X4->O4 T1->X2 T2->X4 I1 I1 I1->X1 I2 I2 I2->X2 I3 I3 I3->X3 I4 I4 I4->X4

Model Composition

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.

In [9]:
display(compose(spontaneous,spontaneous))
Domain: [1, 2]
Codomain: [1, 2]
Model: Model(S=[1, 2], Δ=[X₁()→X₂(), X₁()→X₂()])


Out[9]:
G X1 X1 T1 T1 X1->T1 T2 T2 X1->T2 O1 O1 X1->O1 X2 X2 O2 O2 X2->O2 T1->X2 T2->X2 I1 I1 I1->X1 I2 I2 I2->X2

Introducing new variables

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.

In [108]:
display(eye(1))
Domain: [1]
Codomain: [1]
Model: Model(S=[1], Δ=[])

Out[108]:
G X1 X1 O1 O1 X1->O1 I1 I1 I1->X1
In [10]:
display(compose(spontaneous  eye(1), eye(1)  spontaneous))
Domain: [1, 2, 3]
Codomain: [1, 2, 3]
Model: Model(S=[1, 2, 3], Δ=[X₁()→X₂(), X₂()→X₃()])


Out[10]:
G X1 X1 T1 T1 X1->T1 O1 O1 X1->O1 X2 X2 T2 T2 X2->T2 O2 O2 X2->O2 X3 X3 O3 O3 X3->O3 T1->X2 T2->X3 I1 I1 I1->X1 I2 I2 I2->X2 I3 I3 I3->X3

Adding new primitives

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

In [11]:
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)
Infection reaction infect = X₁+X₂→ 2X₂
Domain: [1, 2]
Codomain: [1, 2]
Model: Model(S=[1, 2], Δ=[X₁() + X₂()→2 * X₂(), X₁() + X₂()→2 * X₂()])


Out[11]:
G X1 X1 T1 T1 X1->T1 O1 O1 X1->O1 X2 X2 X2->T1 O2 O2 X2->O2 T1->X2 2 I1 I1 I1->X1 I2 I2 I2->X2

Swapping Arguments

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.

In [83]:
σ() = OpenModel([1,2], NullModel(2), [2,1])
display(σ())
Domain: [1, 2]
Codomain: [2, 1]
Model: Model(S=[1, 2], Δ=[])

Out[83]:
G X1 X1 O2 O2 X1->O2 X2 X2 O1 O1 X2->O1 I1 I1 I1->X1 I2 I2 I2->X2
In [12]:
spontaneousdag = compose(σ(), spontaneous)
display(spontaneousdag)
Domain: [1, 2]
Codomain: [2, 1]
Model: Model(S=[1, 2], Δ=[X₂()→X₁(), X₂()→X₁()])


Out[12]:
G X1 X1 O2 O2 X1->O2 X2 X2 T1 T1 X2->T1 O1 O1 X2->O1 T1->X1 I1 I1 I1->X1 I2 I2 I2->X2

Model Execution

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.

In [13]:
sis = compose(infect, spontaneousdag)
@show odeexpr(sis)
display(sis)
odeexpr(sis) = quote
    ##state#361(du, u, T, t) = begin
            begin
                du.X₁ = (-1 * T[1] * u.X₁ * u.X₂) / N(u) + T[2] * u.X₂
                du.X₂ = (2 * T[1] * u.X₁ * u.X₂) / N(u) + (-1 * T[1] * u.X₁ * u.X₂) / N(u) + -1 * T[2] * u.X₂
            end
        end
end
Domain: [1, 2]
Codomain: [2, 1]
Model: Model(S=[1, 2], Δ=[X₁() + X₂()→2 * X₂(), X₂()→X₁()])


Out[13]:
G X1 X1 T1 T1 X1->T1 O2 O2 X1->O2 X2 X2 X2->T1 T2 T2 X2->T2 O1 O1 X2->O1 T1->X2 2 T2->X1 I1 I1 I1->X1 I2 I2 I2->X2

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

In [14]:
sir = compose(infect⊗eye(1), eye(1)⊗spontaneous)
@show odeexpr(sir)
display(sir)
odeexpr(sir) = quote
    ##state#362(du, u, T, t) = begin
            begin
                du.X₁ = (-1 * T[1] * u.X₁ * u.X₂) / N(u)
                du.X₃ = T[2] * u.X₂
                du.X₂ = (2 * T[1] * u.X₁ * u.X₂) / N(u) + (-1 * T[1] * u.X₁ * u.X₂) / N(u) + -1 * T[2] * u.X₂
            end
        end
end
Domain: [1, 2, 3]
Codomain: [1, 2, 3]
Model: Model(S=[1, 2, 3], Δ=[X₁() + X₂()→2 * X₂(), X₂()→X₃()])


Out[14]:
G X1 X1 T1 T1 X1->T1 O1 O1 X1->O1 X2 X2 X2->T1 T2 T2 X2->T2 O2 O2 X2->O2 X3 X3 O3 O3 X3->O3 T1->X2 2 T2->X3 I1 I1 I1->X1 I2 I2 I2->X2 I3 I3 I3->X3

A modeling language

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.

In [15]:
display(sir⊗sir)
Domain: [1, 2, 3, 4, 5, 6]
Codomain: [1, 2, 3, 4, 5, 6]
Model: Model(S=[1, 2, 3, 4, 5, 6], Δ=[X₁() + X₂()→2 * X₂(), X₂()→X₃(), X₄() + X₅()→2 * X₅(), X₅()→X₆()])


Out[15]:
G X1 X1 T1 T1 X1->T1 O1 O1 X1->O1 X2 X2 X2->T1 T2 T2 X2->T2 O2 O2 X2->O2 X3 X3 O3 O3 X3->O3 X4 X4 T3 T3 X4->T3 O4 O4 X4->O4 X5 X5 X5->T3 T4 T4 X5->T4 O5 O5 X5->O5 X6 X6 O6 O6 X6->O6 T1->X2 2 T2->X3 T3->X5 2 T4->X6 I1 I1 I1->X1 I2 I2 I2->X2 I3 I3 I3->X3 I4 I4 I4->X4 I5 I5 I5->X5 I6 I6 I6->X6

Merging States

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.

In [16]:
m₁₂ = OpenModel([1], NullModel(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′)
Parallel Infections reactions infect ⊗ infect = X₁+X₂→ 2X₂ && X₃ +X₄ → 2X₄
Domain: [1, 3]
Codomain: [2, 4]
Model: Model(S=[1, 2, 3, 4], Δ=[X₁() + X₂()→2 * X₂(), X₃() + X₄()→2 * X₄()])


Out[16]:
G X1 X1 T1 T1 X1->T1 X2 X2 X2->T1 O1 O1 X2->O1 X3 X3 T2 T2 X3->T2 X4 X4 X4->T2 O2 O2 X4->O2 T1->X2 2 T2->X4 2 I1 I1 I1->X1 I2 I2 I2->X3
In [17]:
display(m₁₂)
Domain: [1]
Codomain: [1, 1]
Model: Model(S=[1], Δ=[])

Out[17]:
G X1 X1 O1 O1 X1->O1 O2 O2 X1->O2 I1 I1 I1->X1
In [94]:
display(compose(m₁₂, parinfect′))
Domain: [1]
Codomain: [2, 3]
Model: Model(S=[1, 2, 3], Δ=[X₁() + X₂()→2 * X₂(), X₁() + X₃()→2 * X₃()])


Out[94]:
G X1 X1 T1 T1 X1->T1 T2 T2 X1->T2 X2 X2 X2->T1 O1 O1 X2->O1 X3 X3 X3->T2 O2 O2 X3->O2 T1->X2 2 T2->X3 2 I1 I1 I1->X1
In [91]:
mutrec = OpenModel([1,2], Model([1,2,3], [(X[1],X[3]), (X[2], X[3])]), [3])
display(mutrec)
Domain: [1, 2]
Codomain: [3]
Model: Model(S=[1, 2, 3], Δ=[X₁()→X₃(), X₂()→X₃()])


Out[91]:
G X1 X1 T1 T1 X1->T1 X2 X2 T2 T2 X2->T2 X3 X3 O1 O1 X3->O1 T1->X3 T2->X3 I1 I1 I1->X1 I2 I2 I2->X2
In [92]:
siir = compose(m₁₂, parinfect′, mutrec)
display(siir)
Domain: [1]
Codomain: [4]
Model: Model(S=[1, 2, 3, 4], Δ=[X₁() + X₂()→2 * X₂(), X₁() + X₃()→2 * X₃(), X₂()→X₄(), X₃()→X₄()])


Out[92]:
G X1 X1 T1 T1 X1->T1 T2 T2 X1->T2 X2 X2 X2->T1 T3 T3 X2->T3 X3 X3 X3->T2 T4 T4 X3->T4 X4 X4 O1 O1 X4->O1 T1->X2 2 T2->X3 2 T3->X4 T4->X4 I1 I1 I1->X1
In [75]:
odeexpr(siir)
Out[75]:
quote
    ##state#462(du, u, T, t) = begin
            begin
                du.X₁ = (-1 * T[1] * u.X₁ * u.X₂) / N(u) + (-1 * T[2] * u.X₁ * u.X₃) / N(u)
                du.X₃ = (2 * T[2] * u.X₁ * u.X₃) / N(u) + (-1 * T[2] * u.X₁ * u.X₃) / N(u) + -1 * T[4] * u.X₃
                du.X₂ = (2 * T[1] * u.X₁ * u.X₂) / N(u) + (-1 * T[1] * u.X₁ * u.X₂) / N(u) + -1 * T[3] * u.X₂
                du.X₄ = T[3] * u.X₂ + T[4] * u.X₃
            end
        end
end

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.

Abstract Knowledge Layer

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.

Ecological Models

In [20]:
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)
Creating food web processes birth, death, predation

Creating food web processes σ, predation†
Out[20]:
compile (generic function with 4 methods)

The Lotka Volterra Model

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.

Wiring Diagrams

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.

In [107]:
bd = b⊗d
lvh = compose(bh⊗idₓ, ph, idₓ⊗dh)
drawhom(lvh, "img/lv_wd")
display(wd(lvh), [:Wolf, :Sheep, :Sheep, :Wolf, :Wolf,:Sheep])
Out[107]:
G n5 birth n1p1:e->n5:w Sheep n3 predation n1p2:e->n3:w Wolf n3:e->n2p1:w Sheep n4 death n3:e->n4:w Wolf n4:e->n2p2:w Wolf n5:e->n3:w Sheep

Compiling a Wiring Diagram

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

In [99]:
lvh
Out[99]:
$\left(\mathrm{birth} \otimes \mathrm{id}_{X}\right) \cdot \mathrm{predation} \cdot \left(\mathrm{id}_{X} \otimes \mathrm{death}\right) : X \otimes X \to X \otimes X$
In [100]:
lv = compile(lvh)
println(lv)
debug_show(lv,"img/lv.svg")
OpenModel:2→2 with Model(∣S∣=3,∣Δ∣=3)
Domain: [1, 2]
Codomain: [1, 2]
Model: Model(S=[1, 2, 3], Δ=[X₁()→2 * X₁(), 1 * X₁() + 1 * X₂()→1.15 * X₂(), X₂()→X₃()])


Out[100]:
G X1 X1 T1 T1 X1->T1 T2 T2 X1->T2 1 O1 O1 X1->O1 X2 X2 X2->T2 1 T3 T3 X2->T3 O2 O2 X2->O2 X3 X3 T1->X1 2 T2->X2 1.15 T3->X3 I1 I1 I1->X1 I2 I2 I2->X2

Three Layers

At this point we can see all three layers of the Lotka-Volterra model. Wiring Diagram -> Petri Net -> ODEs

In [101]:
@show odeexpr(lv)
f = Petri.mk_function(lv.model)
odeexpr(lv) = quote
    ##state#463(du, u, T, t) = begin
            begin
                du.X₁ = 2 * T[1] * u.X₁ + -1 * T[1] * u.X₁ + (-1 * T[2] * u.X₁ * u.X₂) / N(u)
                du.X₃ = T[3] * u.X₂
                du.X₂ = (1.15 * T[2] * u.X₁ * u.X₂) / N(u) + (-1 * T[2] * u.X₁ * u.X₂) / N(u) + -1 * T[3] * u.X₂
            end
        end
end
Out[101]:
#71 (generic function with 1 method)
In [103]:
ν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
sol[end] = [X₁=44.36119914675906, X₂=5.8649499704357995, X₃=63.78127334096779]
Out[103]:
0 2 4 6 8 0 10 20 30 40 50 60 t Sheep Wolves Dwolves

Putting it all together

The Lotka-Volterra Model is a classic example of dynamical systems in mathematical ecology. We can view it at three layers of knowledge.

  1. Abstract: a wiring diagram with the interaction structure, we have no mathematical or computational knowledge about the dynamics.
  2. Structured: a Petri Net describing reaction kinetics dynamics
  3. Executable: a program for solving the system of ODEs implied by the structured layer

Kind of overkill for such a simple model. But we will benefit as complexity scales up!

Food Chain

The first predator is the second prey

In [78]:
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])
bipredation is (p⊗I)⋅(I⊗p)

foodchain is (bipredation)⋅(bdd). A fish, a bigger fish, and biggest fish
Cannonical form construction proves:  compose(otimes(predation,id(X)),otimes(id(X),predation),otimes(birth,death,death)) == compose(otimes(predation,id(X)),otimes(birth,compose(predation,otimes(death,death))))
Out[78]:
G n5 predation n1p1:e->n5:w fish n1p2:e->n5:w Fish n6 predation n1p3:e->n6:w Shark n3 death n3:e->n2p3:w n4 birth n4:e->n2p1:w n5:e->n4:w n5:e->n6:w n6:e->n3:w n7 death n6:e->n7:w n7:e->n2p2:w
In [26]:
foodchain = compile(foodchainh)
display(foodchain)
Domain: [1, 2, 3]
Codomain: [1, 2, 3]
Model: Model(S=[1, 2, 3, 4, 5], Δ=[1 * X₁() + 1 * X₂()→1.15 * X₂(), 1 * X₂() + 1 * X₃()→1.15 * X₃(), X₁()→2 * X₁(), X₂()→X₄(), X₃()→X₅()])


Out[26]:
G X1 X1 T1 T1 X1->T1 1 T3 T3 X1->T3 O1 O1 X1->O1 X2 X2 X2->T1 1 T2 T2 X2->T2 1 T4 T4 X2->T4 O2 O2 X2->O2 X3 X3 X3->T2 1 T5 T5 X3->T5 O3 O3 X3->O3 X4 X4 X5 X5 T1->X2 1.15 T2->X3 1.15 T3->X1 2 T4->X4 T5->X5 I1 I1 I1->X1 I2 I2 I2->X2 I3 I3 I3->X3
In [27]:
println("As an ordinary differential equation:")
odeexpr(foodchain).args[end].args[end].args[end]
As an ordinary differential equation:
Out[27]:
quote
    du.X₁ = (-1 * T[1] * u.X₁ * u.X₂) / N(u) + 2 * T[3] * u.X₁ + -1 * T[3] * u.X₁
    du.X₅ = T[5] * u.X₃
    du.X₃ = (1.15 * T[2] * u.X₂ * u.X₃) / N(u) + (-1 * T[2] * u.X₂ * u.X₃) / N(u) + -1 * T[5] * u.X₃
    du.X₂ = (1.15 * T[1] * u.X₁ * u.X₂) / N(u) + (-1 * T[1] * u.X₁ * u.X₂) / N(u) + (-1 * T[2] * u.X₂ * u.X₃) / N(u) + -1 * T[4] * u.X₂
    du.X₄ = T[4] * u.X₂
end

Food Star

The first predator is the second predator (with two independent prey species)

In [28]:
foodstarh = compose(ph⊗idₓ, idₓ⊗compose(σh,ph, σh), bh⊗dh⊗bh)
display(wd(foodstarh), [:Rabbit, :Wolf, :Sheep])
Out[28]:
G n5 predation n1p1:e->n5:w Rabbit n1p2:e->n5:w Wolf n6 predation n1p3:e->n6:w Sheep n3 birth n3:e->n2p3:w n4 birth n4:e->n2p1:w n5:e->n4:w n5:e->n6:w n6:e->n3:w n7 death n6:e->n7:w n7:e->n2p2:w
In [29]:
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)
OpenModel:3→3 with Model(∣S∣=4,∣Δ∣=5)
Cannonical form construction proves:  compose(otimes(predation,id(X)),otimes(id(X),compose(braid(X,X),predation,braid(X,X))),otimes(birth,death,birth)) == compose(otimes(predation,id(X)),otimes(id(X),braid(X,X)),otimes(birth,compose(predation,otimes(birth,death))),otimes(id(X),braid(X,X)))
As an ordinary differential equation:
odeexpr(foodstar) = quote
    ##state#367(du, u, T, t) = begin
            begin
                du.X₁ = (-1 * T[1] * u.X₁ * u.X₂) / N(u) + 2 * T[3] * u.X₁ + -1 * T[3] * u.X₁
                du.X₃ = (-1 * T[2] * u.X₃ * u.X₂) / N(u) + 2 * T[5] * u.X₃ + -1 * T[5] * u.X₃
                du.X₂ = (1.15 * T[1] * u.X₁ * u.X₂) / N(u) + (-1 * T[1] * u.X₁ * u.X₂) / N(u) + (1.15 * T[2] * u.X₃ * u.X₂) / N(u) + (-1 * T[2] * u.X₃ * u.X₂) / N(u) + -1 * T[4] * u.X₂
                du.X₄ = T[4] * u.X₂
            end
        end
end
Domain: [1, 2, 3]
Codomain: [1, 2, 3]
Model: Model(S=[1, 2, 3, 4], Δ=[1 * X₁() + 1 * X₂()→1.15 * X₂(), 1 * X₃() + 1 * X₂()→1.15 * X₂(), X₁()→2 * X₁(), X₂()→X₄(), X₃()→2 * X₃()])


Out[29]:
G X1 X1 T1 T1 X1->T1 1 T3 T3 X1->T3 O1 O1 X1->O1 X2 X2 X2->T1 1 T2 T2 X2->T2 1 T4 T4 X2->T4 O2 O2 X2->O2 X3 X3 X3->T2 1 T5 T5 X3->T5 O3 O3 X3->O3 X4 X4 T1->X2 1.15 T2->X2 1.15 T3->X1 2 T4->X4 T5->X3 2 I1 I1 I1->X1 I2 I2 I2->X2 I3 I3 I3->X3

Malaria

In [30]:
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")
Out[30]:
$\left(inf¹³₂₃ \otimes \mathrm{id}_{X}\right) \cdot \left(\mathrm{id}_{X} \otimes \left(\left(\sigma_{X,X} \otimes \mathrm{id}_{X}\right) \cdot \left(\mathrm{id}_{X} \otimes \sigma_{X,X}\right) \cdot \left(\sigma_{X,X} \otimes \mathrm{id}_{X}\right)\right)\right) \cdot \left(\mathrm{id}_{X} \otimes inf¹³₂₃\right) \cdot \left(\mathrm{id}_{X} \otimes \left(\left(\mathrm{id}_{X} \otimes \sigma_{X,X}\right) \cdot \left(\sigma_{X,X} \otimes \mathrm{id}_{X}\right) \cdot \left(\mathrm{id}_{X} \otimes \sigma_{X,X}\right)\right)\right) \cdot \left(\mathrm{cure} \otimes \left(\sigma_{X,X} \cdot \mathrm{cure} \cdot \sigma_{X,X}\right)\right) : X \otimes X \otimes X \otimes X \to X \otimes X \otimes X \otimes X$
In [31]:
display(wd(malariah), ["Sp", "Ip", "Im", "Sm"])
Out[31]:
G n5 inf¹³₂₃ n1p1:e->n5:w Sp n1p2:e->n5:w Ip n1p3:e->n5:w Im n6 inf¹³₂₃ n1p4:e->n6:w Sm n3 cure n3:e->n2p1:w n3:e->n2p2:w n4 cure n4:e->n2p3:w n4:e->n2p4:w n5:e->n3:w n5:e->n6:w n5:e->n6:w n6:e->n3:w n6:e->n4:w n6:e->n4:w
In [32]:
cure = conj(spontaneous, σ())
display(cure)
Domain: [1, 2]
Codomain: [1, 2]
Model: Model(S=[1, 2], Δ=[X₂()→X₁(), X₂()→X₁()])


Out[32]:
G X1 X1 O1 O1 X1->O1 X2 X2 T1 T1 X2->T1 O2 O2 X2->O2 T1->X1 I1 I1 I1->X1 I2 I2 I2->X2
In [33]:
curedag = op(cure)
curedag == spontaneous
Out[33]:
true

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$

In [34]:
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)
Domain: [1, 2, 3]
Codomain: [1, 2, 3]
Model: Model(S=[1, 2, 3], Δ=[X₁() + X₃()→X₂() + X₃(), X₁() + X₃()→X₂() + X₃()])


Out[34]:
G X1 X1 T1 T1 X1->T1 O1 O1 X1->O1 X2 X2 O2 O2 X2->O2 X3 X3 X3->T1 O3 O3 X3->O3 T1->X2 T1->X3 I1 I1 I1->X1 I2 I2 I2->X2 I3 I3 I3->X3
In [35]:
inf′ = compose(σh⊗idₓ, idₓ⊗σh, σh⊗idₓ, inf, conj(σh⊗idₓ, idₓ⊗σh))
display(wd(inf′))
Out[35]:
G n3 inf¹³₂₃ n1p1:e->n3:w n1p2:e->n3:w n1p3:e->n3:w n3:e->n2p1:w n3:e->n2p2:w n3:e->n2p3:w
In [36]:
dualinfecth = compose(inf⊗idₓ, idₓ⊗inf′)
display(wd(dualinfecth), ["Sp", "Ip","Im","Sm"])
Out[36]:
G n3 inf¹³₂₃ n1p1:e->n3:w Sp n1p2:e->n3:w Ip n1p3:e->n3:w Im n4 inf¹³₂₃ n1p4:e->n4:w Sm n3:e->n2p1:w n3:e->n4:w n3:e->n4:w n4:e->n2p2:w n4:e->n2p3:w n4:e->n2p4:w
In [37]:
dualinfect = compile(dualinfecth)
display(dualinfect)
Domain: [1, 2, 3, 4]
Codomain: [1, 2, 3, 4]
Model: Model(S=[1, 2, 3, 4], Δ=[X₁() + X₃()→X₂() + X₃(), X₄() + X₂()→X₃() + X₂()])


Out[37]:
G X1 X1 T1 T1 X1->T1 O1 O1 X1->O1 X2 X2 T2 T2 X2->T2 O2 O2 X2->O2 X3 X3 X3->T1 O3 O3 X3->O3 X4 X4 X4->T2 O4 O4 X4->O4 T1->X2 T1->X3 T2->X2 T2->X3 I1 I1 I1->X1 I2 I2 I2->X2 I3 I3 I3->X3 I4 I4 I4->X4
In [38]:
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")
Cannonical form construction proves:
	  compose(otimes(inf¹³₂₃,id(X)),otimes(id(X),compose(otimes(braid(X,X),id(X)),otimes(id(X),braid(X,X)),otimes(braid(X,X),id(X)))),otimes(id(X),inf¹³₂₃),otimes(id(X),compose(otimes(id(X),braid(X,X)),otimes(braid(X,X),id(X)),otimes(id(X),braid(X,X)))),otimes(cure,compose(braid(X,X),cure,braid(X,X)))) 
	== 
	  compose(otimes(inf¹³₂₃,id(X)),otimes(id(X),braid(X,X),id(X)),otimes(id(otimes(X,X)),braid(X,X)),otimes(id(X),braid(X,X),id(X)),otimes(id(X),inf¹³₂₃),otimes(id(otimes(X,X)),braid(X,X)),otimes(id(X),braid(X,X),id(X)),otimes(cure,cure),otimes(id(otimes(X,X)),braid(X,X)))
As an ordinary differential equation:
odeexpr(malaria) = quote
    ##state#368(du, u, T, t) = begin
            begin
                du.X₁ = (-1 * T[1] * u.X₁ * u.X₃) / N(u) + T[3] * u.X₂
                du.X₃ = (T[1] * u.X₁ * u.X₃) / N(u) + (-1 * T[1] * u.X₁ * u.X₃) / N(u) + (T[2] * u.X₄ * u.X₂) / N(u) + -1 * T[4] * u.X₃
                du.X₂ = (T[1] * u.X₁ * u.X₃) / N(u) + (T[2] * u.X₄ * u.X₂) / N(u) + (-1 * T[2] * u.X₄ * u.X₂) / N(u) + -1 * T[3] * u.X₂
                du.X₄ = (-1 * T[2] * u.X₄ * u.X₂) / N(u) + T[4] * u.X₃
            end
        end
end
Domain: [1, 2, 3, 4]
Codomain: [1, 2, 3, 4]
Model: Model(S=[1, 2, 3, 4], Δ=[X₁() + X₃()→X₂() + X₃(), X₄() + X₂()→X₃() + X₂(), X₂()→X₁(), X₃()→X₄()])


Out[38]:
G X1 X1 T1 T1 X1->T1 O1 O1 X1->O1 X2 X2 T2 T2 X2->T2 T3 T3 X2->T3 O2 O2 X2->O2 X3 X3 X3->T1 T4 T4 X3->T4 O3 O3 X3->O3 X4 X4 X4->T2 O4 O4 X4->O4 T1->X2 T1->X3 T2->X2 T2->X3 T3->X1 T4->X4 I1 I1 I1->X1 I2 I2 I2->X2 I3 I3 I3->X3 I4 I4 I4->X4
In [39]:
f = Petri.mk_function(malaria.model)
Out[39]:
#71 (generic function with 1 method)
In [40]:
u0 = @LArray [20,1,0,100.0] (:X₁, :X₂, :X₃, :X₄)
varnames = ["\$S_p\$", "\$I_p\$", "\$I_m\$", "\$S_m\$"];

Endemic Malaria

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.

In [41]:
β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
sol[end] = [X₁=15.732700855242287, X₂=5.267299144757728, X₃=37.20917828829352, X₄=62.79082171170651]
Out[41]:
0 100 200 300 0 25 50 75 100 t Sp Ip Im Sm

Malaria Driven Extinction

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
In [42]:
malvitalh = compose(malariah, bh⊗dh⊗dh⊗bh)
Out[42]:
$\left(inf¹³₂₃ \otimes \mathrm{id}_{X}\right) \cdot \left(\mathrm{id}_{X} \otimes \left(\left(\sigma_{X,X} \otimes \mathrm{id}_{X}\right) \cdot \left(\mathrm{id}_{X} \otimes \sigma_{X,X}\right) \cdot \left(\sigma_{X,X} \otimes \mathrm{id}_{X}\right)\right)\right) \cdot \left(\mathrm{id}_{X} \otimes inf¹³₂₃\right) \cdot \left(\mathrm{id}_{X} \otimes \left(\left(\mathrm{id}_{X} \otimes \sigma_{X,X}\right) \cdot \left(\sigma_{X,X} \otimes \mathrm{id}_{X}\right) \cdot \left(\mathrm{id}_{X} \otimes \sigma_{X,X}\right)\right)\right) \cdot \left(\mathrm{cure} \otimes \left(\sigma_{X,X} \cdot \mathrm{cure} \cdot \sigma_{X,X}\right)\right) \cdot \left(\mathrm{birth} \otimes \mathrm{death} \otimes \mathrm{death} \otimes \mathrm{birth}\right) : X \otimes X \otimes X \otimes X \to X \otimes X \otimes X \otimes X$
In [43]:
display(wd(malvitalh), ["Sp", "Ip","Im","Sm"])
Out[43]:
G n7 inf¹³₂₃ n1p1:e->n7:w Sp n1p2:e->n7:w Ip n1p3:e->n7:w Im n8 inf¹³₂₃ n1p4:e->n8:w Sm n3 death n3:e->n2p2:w n4 death n4:e->n2p3:w n5 cure n5:e->n3:w n9 birth n5:e->n9:w n6 cure n6:e->n4:w n10 birth n6:e->n10:w n7:e->n5:w n7:e->n8:w n7:e->n8:w n8:e->n5:w n8:e->n6:w n8:e->n6:w n9:e->n2p1:w n10:e->n2p4:w
In [44]:
malvital = compile(malvitalh); #compose(malaria, b⊗d⊗d⊗b)
display(malvital)
Domain: [1, 2, 3, 4]
Codomain: [1, 2, 3, 4]
Model: Model(S=[1, 2, 3, 4, 5, 6], Δ=[X₁() + X₃()→X₂() + X₃(), X₄() + X₂()→X₃() + X₂(), X₂()→X₁(), X₃()→X₄(), X₁()→2 * X₁(), X₂()→X₅(), X₃()→X₆(), X₄()→2 * X₄()])


Out[44]:
G X1 X1 T1 T1 X1->T1 T5 T5 X1->T5 O1 O1 X1->O1 X2 X2 T2 T2 X2->T2 T3 T3 X2->T3 T6 T6 X2->T6 O2 O2 X2->O2 X3 X3 X3->T1 T4 T4 X3->T4 T7 T7 X3->T7 O3 O3 X3->O3 X4 X4 X4->T2 T8 T8 X4->T8 O4 O4 X4->O4 X5 X5 X6 X6 T1->X2 T1->X3 T2->X2 T2->X3 T3->X1 T4->X4 T5->X1 2 T6->X5 T7->X6 T8->X4 2 I1 I1 I1->X1 I2 I2 I2->X2 I3 I3 I3->X3 I4 I4 I4->X4
In [45]:
println("As an ordinary differential equation:")
@show fex = odeexpr(malvital)
f = Petri.mk_function(malvital.model)
As an ordinary differential equation:
fex = odeexpr(malvital) = quote
    ##state#371(du, u, T, t) = begin
            begin
                du.X₁ = (-1 * T[1] * u.X₁ * u.X₃) / N(u) + T[3] * u.X₂ + 2 * T[5] * u.X₁ + -1 * T[5] * u.X₁
                du.X₅ = T[6] * u.X₂
                du.X₃ = (T[1] * u.X₁ * u.X₃) / N(u) + (-1 * T[1] * u.X₁ * u.X₃) / N(u) + (T[2] * u.X₄ * u.X₂) / N(u) + -1 * T[4] * u.X₃ + -1 * T[7] * u.X₃
                du.X₂ = (T[1] * u.X₁ * u.X₃) / N(u) + (T[2] * u.X₄ * u.X₂) / N(u) + (-1 * T[2] * u.X₄ * u.X₂) / N(u) + -1 * T[3] * u.X₂ + -1 * T[6] * u.X₂
                du.X₄ = (-1 * T[2] * u.X₄ * u.X₂) / N(u) + T[4] * u.X₃ + 2 * T[8] * u.X₄ + -1 * T[8] * u.X₄
                du.X₆ = T[7] * u.X₃
            end
        end
end
Out[45]:
#71 (generic function with 1 method)
In [46]:
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\$"];
In [47]:
β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
sol[end] = [X₁=0.118672823642017, X₂=0.6300807358280642, X₃=15.66598215784487, X₄=115.38905945308686, X₅=25.773026412722405, X₆=152.68431308544328]
Out[47]:
0 10 20 30 40 0 25 50 75 100 t Sp Ip Im Sm

Endemic Malaria with Vital dynamics

In [48]:
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\$"];
In [49]:
β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
sol[end] = [X₁=7.925128113171479, X₂=63.53298752007829, X₃=116.96980173521904, X₄=24.601033560805966, X₅=495.33260872620883, X₆=345.1517256206875]
Out[49]:
0 50 100 150 0 25 50 75 100 t Sp Ip Im Sm

Mosquito Hunting Birds

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

In [50]:
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])
Out[50]:
G n5 malaria n1p1:e->n5:w Sp n1p2:e->n5:w Ip n1p3:e->n5:w Im n1p4:e->n5:w Sm n3 mutualpredation n1p5:e->n3:w B n4 vitals n3:e->n4:w n3:e->n4:w n3:e->n4:w n4:e->n2p1:w n4:e->n2p2:w n4:e->n2p3:w n4:e->n2p4:w n4:e->n2p5:w n5:e->n3:w n5:e->n3:w n5:e->n4:w n5:e->n4:w
In [51]:
birdsh = compose(σh⊗idₓ,idₓ⊗ph, σh⊗idₓ, idₓ⊗ph)
display(wd(birdsh), [:Im, :B, :Sm])
Out[51]:
G n3 predation n1p1:e->n3:w Im n4 predation n1p2:e->n4:w Sm n1p3:e->n3:w B n3:e->n2p1:w n3:e->n4:w n4:e->n2p2:w n4:e->n2p3:w
In [52]:
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])
Cannonical form construction proves:  compose(otimes(braid(X,X),id(X)),otimes(id(X),predation),otimes(braid(X,X),id(X)),otimes(id(X),predation)) == compose(otimes(braid(X,X),id(X)),otimes(id(X),predation),otimes(braid(X,X),id(X)),otimes(id(X),predation))
Out[52]:
G n7 inf¹³₂₃ n1p1:e->n7:w Sp n1p2:e->n7:w Ip n1p3:e->n7:w Im n8 inf¹³₂₃ n1p4:e->n8:w Sm n5 predation n1p5:e->n5:w B n3 birth n3:e->n2p4:w Sm n4 death n4:e->n2p2:w Ip n6 predation n5:e->n6:w B n11 death n5:e->n11:w Sm n6:e->n3:w Sm n12 death n6:e->n12:w B n7:e->n8:w Im n7:e->n8:w Ip n9 cure n7:e->n9:w Sp n8:e->n9:w Ip n10 cure n8:e->n10:w Sm n8:e->n10:w Im n9:e->n4:w Ip n13 birth n9:e->n13:w Sp n10:e->n5:w Im n10:e->n6:w Sm n11:e->n2p3:w Im n12:e->n2p5:w B n13:e->n2p1:w Sp
In [53]:
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");
Domain: [1, 2, 3]
Codomain: [1, 2, 3]
Model: Model(S=[1, 2, 3], Δ=[1 * X₁() + 1 * X₃()→1.15 * X₃(), 1 * X₂() + 1 * X₃()→1.15 * X₃()])


Domain: [1, 2, 3, 4, 5]
Codomain: [1, 2, 3, 4, 5]
Model: Model(S=[1, 2, 3, 4, 5, 6, 7, 8], Δ=[X₁() + X₃()→X₂() + X₃(), X₄() + X₂()→X₃() + X₂(), X₂()→X₁(), X₃()→X₄(), 1 * X₃() + 1 * X₅()→1.15 * X₅(), 1 * X₄() + 1 * X₅()→1.15 * X₅(), X₁()→2 * X₁(), X₂()→X₆(), X₃()→X₇(), X₄()→2 * X₄(), X₅()→X₈()])


The coupled model of Malaria and Predation:

The coupled model of Malaria

In [54]:
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)
As an ordinary differential equation:
odeexpr(birdmal) = quote
    ##state#375(du, u, T, t) = begin
            begin
                du.X₈ = T[11] * u.X₅
                du.X₁ = (-1 * T[1] * u.X₁ * u.X₃) / N(u) + T[3] * u.X₂ + 2 * T[7] * u.X₁ + -1 * T[7] * u.X₁
                du.X₅ = (1.15 * T[5] * u.X₃ * u.X₅) / N(u) + (-1 * T[5] * u.X₃ * u.X₅) / N(u) + (1.15 * T[6] * u.X₄ * u.X₅) / N(u) + (-1 * T[6] * u.X₄ * u.X₅) / N(u) + -1 * T[11] * u.X₅
                du.X₃ = (T[1] * u.X₁ * u.X₃) / N(u) + (-1 * T[1] * u.X₁ * u.X₃) / N(u) + (T[2] * u.X₄ * u.X₂) / N(u) + -1 * T[4] * u.X₃ + (-1 * T[5] * u.X₃ * u.X₅) / N(u) + -1 * T[9] * u.X₃
                du.X₂ = (T[1] * u.X₁ * u.X₃) / N(u) + (T[2] * u.X₄ * u.X₂) / N(u) + (-1 * T[2] * u.X₄ * u.X₂) / N(u) + -1 * T[3] * u.X₂ + -1 * T[8] * u.X₂
                du.X₄ = (-1 * T[2] * u.X₄ * u.X₂) / N(u) + T[4] * u.X₃ + (-1 * T[6] * u.X₄ * u.X₅) / N(u) + 2 * T[10] * u.X₄ + -1 * T[10] * u.X₄
                du.X₆ = T[8] * u.X₂
                du.X₇ = T[9] * u.X₃
            end
        end
end
Out[54]:
#71 (generic function with 1 method)
In [55]:
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\$"];

Malaria Driven Extinction

In the following setting everyone dies, but then the mosquito population becomes locked in a Lotka-Volterra struggle for the earth.

In [56]:
β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
sol[end] = [X₁=0.6518804960409814, X₂=3.9411540896561354e-10, X₃=1.9924967895384285e-8, X₄=56.276722842307045, X₅=0.0011281237872835719, X₆=21.23354569504563, X₇=4.831252597542185, X₈=112.68585807143518]
Out[56]:
0 2 4 6 8 10 0 20 40 60 80 100 t Sp Ip Im Sm B
In [57]:
plt = plot(sol, vars =1:length(varnames)-3, labels=varnames)
#yaxis!(plt, (0,30))
Out[57]:
0 50 100 150 0 20 40 60 80 100 120 t Sp Ip Im Sm B

Endemic Malaria

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.

In [58]:
β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
sol[end] = [X₁=7.712551519101319, X₂=0.6319683940025463, X₃=4.833423048061268, X₄=106.28898842352022, X₅=22.71548262951343, X₆=5139.821778500151, X₇=288.4689530404894, X₈=2948.9339811585596]
Out[58]:
0 25 50 75 100 125 150 175 0 50 100 150 200 250 t Sp Ip Im Sm B

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.

In [59]:
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
Out[59]:
130 140 150 160 170 0 25 50 75 100 125 150 175 t Sp Ip Im Sm B

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.

In [74]:
striplines.(Petri.Δ(birdmal.model, :u))
Out[74]:
11-element Array{Expr,1}:
 :(##δ#451(state) = begin
          begin
              begin
                  u.X₁ > 0 || return nothing
                  u.X₃ > 0 || return nothing
                  u.X₁ -= 1
                  u.X₃ -= 1
              end
              begin
                  u.X₂ += 1
                  u.X₃ += 1
              end
          end
      end)
 :(##δ#452(state) = begin
          begin
              begin
                  u.X₄ > 0 || return nothing
                  u.X₂ > 0 || return nothing
                  u.X₄ -= 1
                  u.X₂ -= 1
              end
              begin
                  u.X₃ += 1
                  u.X₂ += 1
              end
          end
      end)
 :(##δ#453(state) = begin
          begin
              begin
                  u.X₂ > 0 || return nothing
                  u.X₂ -= 1
              end
              u.X₁ += 1
          end
      end)                                                                                                                                               
 :(##δ#454(state) = begin
          begin
              begin
                  u.X₃ > 0 || return nothing
                  u.X₃ -= 1
              end
              u.X₄ += 1
          end
      end)                                                                                                                                               
 :(##δ#455(state) = begin
          begin
              begin
                  u.X₃ > 0 || return nothing
                  u.X₅ > 0 || return nothing
                  u.X₃ -= 1
                  u.X₅ -= 1
              end
              u.X₅ += 1.15
          end
      end)                                                                   
 :(##δ#456(state) = begin
          begin
              begin
                  u.X₄ > 0 || return nothing
                  u.X₅ > 0 || return nothing
                  u.X₄ -= 1
                  u.X₅ -= 1
              end
              u.X₅ += 1.15
          end
      end)                                                                   
 :(##δ#457(state) = begin
          begin
              begin
                  u.X₁ > 0 || return nothing
                  u.X₁ -= 1
              end
              u.X₁ += 2
          end
      end)                                                                                                                                               
 :(##δ#458(state) = begin
          begin
              begin
                  u.X₂ > 0 || return nothing
                  u.X₂ -= 1
              end
              u.X₆ += 1
          end
      end)                                                                                                                                               
 :(##δ#459(state) = begin
          begin
              begin
                  u.X₃ > 0 || return nothing
                  u.X₃ -= 1
              end
              u.X₇ += 1
          end
      end)                                                                                                                                               
 :(##δ#460(state) = begin
          begin
              begin
                  u.X₄ > 0 || return nothing
                  u.X₄ -= 1
              end
              u.X₄ += 2
          end
      end)                                                                                                                                               
 :(##δ#461(state) = begin
          begin
              begin
                  u.X₅ > 0 || return nothing
                  u.X₅ -= 1
              end
              u.X₈ += 1
          end
      end)                                                                                                                                               

Conclusion

  1. SemanticModels.jl represents models at multiple Knowledge Representation Layers.
  2. Composition and Combination are easiest at the Abstract Layer
  3. Structured Layer contains a rigorous mathematical specification of the model
  4. Executable Layer artifacts can be compiled from Structured Layer and exectued to solve complex models
  5. Partially Symbolic Computing + Code Generation $\implies$ A revolution in Scientific Computing!

Future Work

  1. Generalize from Petri Nets to Graphs, DAGs, Relational Ologs, Schemas
  2. Prove stuff
  3. Functional programs as a fallback for code that has dynamics at the structured layer.

Unstable System

The following choice of parameters exhibits an unstable dynamics, without approaching a periodic equilibrium.

In [61]:
β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
sol[end] = [X₁=11.123283077257447, X₂=0.5143960076735042, X₃=0.8427716284564659, X₄=22.263838461449144, X₅=3.5196619156057465, X₆=354.21035812351687, X₇=20.29087996476093, X₈=353.48953852335836]
Out[61]:
0 100 200 300 400 500 0 25 50 75 100 t Sp Ip Im Sm B

The outbreaks exhibit irregular behavior

In [62]:
plt = plot(sol, vars =[2,3], labels=varnames)
#ylims!(plt, (0,10))
#xlims!(plt, (0,450))
Out[62]:
0 100 200 300 400 500 0 5 10 15 20 25 t Sp Ip