ReactiveDynamics.jl

ReactiveDynamics.jl copied to clipboard

I'm trying to reproduce a 'standard' SIR model in ReactiveDynamics, but I get very different results (my hand-coded discrete problem gives the expected result). Can you point me in the direction of where I'm going wrong in defining my model with ReactiveDynamics?

using ReactiveDynamics
using OrdinaryDiffEq
using Random
using Distributions
using Plots

function sir_markov(u,p,t)
    (S, I, R) = u
    (β, γ, N, δt) = p
    infrate = abs(β*I/N*S*δt)
    recrate = abs(γ*I*δt)
    infection = rand(Poisson(infrate))
    recovery = rand(Poisson(recrate))
    S = S-infection
    I = I+infection-recovery
    R = R+recovery
    [S, I, R]
end

tspan = (0.0,40.0)
u0 = [990, 10, 0] # S, I, R
β = 0.5
γ = 0.25
N = 1000.0
δt = 0.1
p = [β, γ, N, δt] # β, γ, N
seed = 1234
Random.seed!(seed)

prob = DiscreteProblem(sir_markov, u0, tspan, p, dt=δt)
sol = solve(prob, FunctionMap())
plot(sol)

# model dynamics
sir_acs = @ReactionNetwork begin
        β*S*I, S+I --> 2I, name=>infection
        γ*I, I --> R, name=>recovery 
end

# simulation parameters
## initial values
@prob_init sir_acs S=990 I=10 R=0
## uncertainty in initial values (Gaussian)
@prob_uncertainty sir_acs S=0.0 I=0.0
## parameters
@prob_params sir_acs β=0.5 γ=0.25
## other arguments passed to the solver
@prob_meta sir_acs tspan=40 dt=0.01
# turn model into a discrete? problem
sir_prob = @problematize sir_acs
# solve the problem over multiple trajectories
sir_sol = @solve sir_prob
# plot the solution
@plot sir_sol