Edit on GitHub

Basic Epidemiology Models

using AlgebraicPetri
using AlgebraicPetri.Epidemiology

using LabelledArrays
using OrdinaryDiffEq
using Plots

using Catlab
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra
using Catlab.Programs.RelationalPrograms

display_uwd(ex) = to_graphviz(ex, box_labels=:name, junction_labels=:variable, edge_attrs=Dict(:len=>".75"));

SIR Model:

define model

sir = @relation (s,i,r) begin
    infection(s,i)
    recovery(i,r)
end
display_uwd(sir)
G n1 infection n6 s n1--n6 n7 i n1--n7 n2 recovery n2--n7 n8 r n2--n8 n3--n6 n4--n7 n5--n8 
p_sir = apex(oapply_epi(sir))
Graph(p_sir)
G S_S S_S T_inf T_inf S_S->T_inf 1 S_I S_I S_I->T_inf 1 T_rec T_rec S_I->T_rec 1 S_R S_R T_inf->S_I 2 T_rec->S_R 1

define initial states and transition rates, then create, solve, and visualize ODE problem

u0 = LVector(S=10, I=1, R=0);
p = LVector(inf=0.4, rec=0.4);

The C-Set representation has direct support for generating a DiffEq vector field

prob = ODEProblem(vectorfield(p_sir),u0,(0.0,7.5),p);
sol = solve(prob,Tsit5())

plot(sol)

SEIR Model:

define model

seir = @relation (s,e,i,r) begin
    exposure(s,i,e)
    illness(e,i)
    recovery(i,r)
end
display_uwd(seir)
G n1 exposure n8 s n1--n8 n9 e n1--n9 n10 i n1--n10 n2 illness n2--n9 n2--n10 n3 recovery n3--n10 n11 r n3--n11 n4--n8 n5--n9 n6--n10 n7--n11 
p_seir = apex(oapply_epi(seir))
Graph(p_seir)
G S_S S_S T_exp T_exp S_S->T_exp 1 S_I S_I S_I->T_exp 1 T_rec T_rec S_I->T_rec 1 S_E S_E T_ill T_ill S_E->T_ill 1 S_R S_R T_exp->S_I 1 T_exp->S_E 1 T_ill->S_I 1 T_rec->S_R 1

define initial states and transition rates, then create, solve, and visualize ODE problem

u0 = LVector(S=10, E=1, I=0, R=0);
p = LVector(exp=.9, ill=.2, rec=.5);

prob = ODEProblem(vectorfield(p_seir),u0,(0.0,15.0),p);
sol = solve(prob,Tsit5())

plot(sol)

SEIRD Model:

define model

seird = @relation (s,e,i,r,d) begin
    exposure(s,i,e)
    illness(e,i)
    recovery(i,r)
    death(i,d)
end
display_uwd(seird)
G n1 exposure n10 s n1--n10 n11 e n1--n11 n12 i n1--n12 n2 illness n2--n11 n2--n12 n3 recovery n3--n12 n13 r n3--n13 n4 death n4--n12 n14 d n4--n14 n5--n10 n6--n11 n7--n12 n8--n13 n9--n14 
p_seird = apex(oapply_epi(seird))
Graph(p_seird)
G S_S S_S T_exp T_exp S_S->T_exp 1 S_I S_I S_I->T_exp 1 T_rec T_rec S_I->T_rec 1 T_death T_death S_I->T_death 1 S_E S_E T_ill T_ill S_E->T_ill 1 S_R S_R S_D S_D T_exp->S_I 1 T_exp->S_E 1 T_ill->S_I 1 T_rec->S_R 1 T_death->S_D 1

define initial states and transition rates, then create, solve, and visualize ODE problem

u0 = LVector(S=10, E=1, I=0, R=0, D=0);
p = LVector(exp=0.9, ill=0.2, rec=0.5, death=0.1);

prob = ODEProblem(vectorfield(p_seird),u0,(0.0,15.0),p);
sol = solve(prob,Tsit5())

plot(sol)