Agent based network model using EventQueueABM

Specific Domains Modelling & Simulations
1

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
        ]
2

Maybe this thread can help you: EventQueueABM example giving incorrect results · Issue #1110 · JuliaDynamics/Agents.jl · GitHub

edit: actually not, you do it the correct way in your propensity function

3

I think that link did help, thank you so much! - the author of that had a contact rate set which i have set and this fixed my scaling issue (happening when i increased the population between 100, 1000 and 10000. I copied the format for my propensities and i seem to be able to match my compartmental model now. Here is what i have come up with:

using Plots  
using DataFrames, StatsPlots
using Graphs, Karnak

using Agents, Random

n_people = 10000

    # agent struct
    @agent struct Person(GraphAgent)
        status::Symbol  # S E I or R
    end

    susceptibility = 0.84
    exposed_duration = 0.5
    infected_duration = 0.125
    recovered_duration = 0.05
    contact_rate = 1.0

    function infect!(agent, model)
        contact = random_nearby_agent(agent, model)
        if contact.status == :I && (rand() <= susceptibility)
            agent.status = :E
        end
    end

    function infection_propensity(agent, model)
        if agent.status == :S 
            return contact_rate
        else
            return 0.0
        end
    end

    function become_infectious!(agent, model)
        agent.status = :I
    end
    
    function become_infectious_propensity(agent, model)
        if agent.status == :E
            return exposed_duration
        else 
            return 0.0
        end
    end

    function recover!(agent, model)
            agent.status = :R
        end

    function recover_propensity(agent, model)
        if agent.status == :I
            return infected_duration
        else
            return 0.0
        end
    end

     function become_susceptible!(agent, model)
            agent.status = :S
    end 

    function become_susceptible_propensity(agent, model)
        if agent.status == :R
            return recovered_duration
        else 
            return 0.0
        end
    end

    # events used
    infect_event = AgentEvent(action! = infect!, propensity = infection_propensity)
    become_infectious_event = AgentEvent(action! = become_infectious!, propensity = become_infectious_propensity)
    recover_event = AgentEvent(action! = recover!, propensity = recover_propensity)
    become_susceptible_event = AgentEvent(action! = become_susceptible!, propensity = become_susceptible_propensity)

    # 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_)

but please let me know if this seems wrong!

1 Like