Hello, I’m trying to use EventQueueABM from agents.jl on a completely connected network to model disease spread (susceptible, exposed, infected and recovered states). Since it’s completely connected, I believe I should be getting the same epicurve output (counts per time step for each of my disease states) as a stochastic compartmental model that uses continuous time. But I must be doing something wrong, or missing something, for my EventQueueABM set up. Any suggestions? - also i have tried to scale the susceptibility and i’ve definitely done this the wrong way but it was giving me a closer curve than without it.
using Plots
using DataFrames, StatsPlots
using Graphs, Karnak
using Agents, Random
n_people = 100
# agent struct
@agent struct Person(GraphAgent)
status::Symbol # S E I or R
end
#Infection propensity function for S agents
function infection_propensity(agent, model)
if agent.status == :S
neighbors = nearby_agents(agent, model, 1)
n_neighbors = length(neighbors)
n_infected = count(n -> n.status == :I, neighbors)
return (1 - ((1 - (susceptibility))^n_infected))
end
return 0.0
end
# propensities
susceptibility = (0.84*10/(n_people))
exposed_duration = 0.5
infected_duration = 0.125
recovered_duration = 0.05
# events
function infect!(agent, model)
if agent.status == :S
agent.status = :E
end
return
end
function become_infectious!(agent, model)
if agent.status == :E
agent.status = :I
end
return
end
function recover!(agent, model)
if agent.status == :I
agent.status = :R
end
return
end
function become_susceptible!(agent, model)
if agent.status == :R
agent.status = :S
end
return
end
# events used
infect_event = AgentEvent(action! = infect!, propensity = infection_propensity)
become_infectious_event = AgentEvent(action! = become_infectious!, propensity = exposed_duration)
recover_event = AgentEvent(action! = recover!, propensity = infected_duration)
become_susceptible_event = AgentEvent(action! = become_susceptible!, propensity = recovered_duration)
# create tuple of events
events = (infect_event, become_infectious_event, recover_event, become_susceptible_event)
#= Initialize model =#
function initialize_model()
# create complete graph with 100 nodes
complete_network = complete_graph(n_people)
network2 = copy(complete_network)
space2 = GraphSpace(network2)
n_infects = n_people * 0.01 #
model = EventQueueABM(Person, events, space2; #rng,
warn = false)
n_infects_int = Int(round(n_infects))
statuses = vcat(fill(:I, n_infects_int), fill(:S, n_people - n_infects_int))
# Shuffle statuses randomly
shuffle!(statuses)
# Add agents with shuffled statuses
for i in 1:n_people
add_agent!(i, Person, model, statuses[i])
end
return model
end
###############################
# to run once
###############################
# Run simulation
@time model = initialize_model()
using DataFrames, StatsBase, Plots
# Define statuses
all_statuses = [:S, :E, :I, :R]
# Define adata
adata = [(a -> a.status == X, count) for X in all_statuses]
# Run simulation
adf,_ = run!(model, 100.0; adata, when = 0.5, dt = 0.01)
# Ensure missing values are replaced with 0s
for col in names(adf, Not(:time))
replace!(adf[!, col], missing => 0)
end
# Plotting
plt_ = plot(xlabel = "Time", ylabel = "Count", lw = 2, legend = :topright,
size = (1000, 800), title = "IBM on completely connected network")
for col in names(adf, Not(:time))
plot!(plt_, adf.time, adf[!, col], label = col)
end
display(plt_)
And my ODE was using population proportions and the same parameters as propensities from the above code:
u0 = [
S => 0.99,
E => 0.0,
I => 0.01,
R => 0.0
]
p = [
β => 0.84,
γ => 0.125,
σ => 0.5,
ω => 0.05
]
