MTK version of SIR model using rates

sir_mtk.jl

# Libraries
using ModelingToolkit
using OrdinaryDiffEq

# Define size of the system
nstates = 3 # S,I,R
nrates = 2 # infection, recovery
# Model specific parameters
@parameters β c γ S₀ I₀ R₀

# Define time
@parameters t

# Rates
@variables N[1:nrates](t)

# Model is with respect to time
@derivatives D'~t

# Stoichiometry matrix
A = [[-1 1 0];
     [0 -1 1]]

# Derived variables
x₀ = [S₀, I₀, R₀]
X = x₀ .+ A'N

# Transition rates
(S,I,R) = X
λ = [β*c*I*S/(S+I+R),γ*I]

# ODE for rates
neqs = D.(N) .~ λ
tspan = (0.0,40.0)
nu0 = [λ[1] => 0.0, λ[2] => 0.0]
p = [β => 0.05,
     c => 10.0,
     γ => 0.25,
     S₀ => 990.0,
     I₀ => 10.0,
     R₀ => 0.0]
nsys = ODESystem(neqs)
nprob = ODEProblem(nsys,nu0,tspan,p)
nsol = solve(nprob,Tsit5(),saveat=1.0)