Agent based network model using EventQueueABM

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
        ]

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

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

Just wondering if you might know the answer to this related query - I’m trying to change my model so that instead of the infection event being performed through a susceptible agent finding infected agents and potentially becoming infected (receiver model), I want infected agents to find susceptible agents and potentially infect them (transmitter model). I cant make this work however, and my resulting epi curves show agents moving from susceptible to exposed as expected, but no movement from E to I. If start all agents as exposed, they move from E to I as expected, so I think the issue may be that there are two events for the infected agents and issues with the scheduling, rather than an issue with the logic of the disease state transition. Is there a way around this? This is my code:

using ModelingToolkit, DifferentialEquations, Plots  
using DataFrames, StatsPlots
using Graphs, Karnak

#############################################
###### one-mode completely connected graph
#############################################
using Agents, Random
using Random: Xoshiro
#rng = Xoshiro(42)

n_people = 1000

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

    # propensities
    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 == :S && (rand() <= susceptibility)
            contact.status = :E
        end
    end

    function infection_propensity(agent, model)
        if agent.status == :I
            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_)

Thank you in advance!!

Hey @zoek, I have taken a look at the code and I think I understood the problem, I think it stems from the fact that you are updating the contact agent in your new model, not the agent the scheduled the event, that means that for how the EventQueueABM model works, it won’t actually schedule a new event for that agent. Only for the non-contact agent. This explains why the “reverse” version worked. I’m not sure if there is an easy fix for this, but I haven’t thought about deeply about.

Thanks for your insight! I think that makes sense - so you reckon you probably can’t mutate an agent in the action! function, unless it’s the agent that has that event scheduled?
That explains why nothing i was trying was working!

Yes, I mean you can mutate also other agents, just that it doesn’t mean a new event will be scheduled for that one, and so, probably this makes it harder to work as expected because the agent which becomes infected doesn’t schedule a new event

1 Like

Oh actually it could be very simple maybe, you could try adding add_event!(contact, model) for the infected agent.

Seems to work better indeed:

1 Like

However, keep in mind that this is still not equivalent to the one in “reverse mode” because here you also schedule automatically an event for agent in infect!(agent, model), so I believe to make it completely equivalent you would need to schedule events manually as I just suggested but for all the function involved in the model i.e. you should set autogenerate_after_action=false and then use add_event! manually.

1 Like

That looks as expected! Could you please explain where you added add_event!(contact, model) ? Should it be within the model_initialize function? Thanks!!

here:

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

because otherwise no event is newly scheduled for the contact agent when it is exposed

1 Like

Thanks! I’m still having trouble with the autogenerate_after_action=false and manually scheduling, though - is that something i have to include within EventQueueABM?

Sorry, a lot of this stuff is beyond my understanding!

Look at the API page of EventQueueABM, you can see that autogenerate_after_action is a keyword of the type. For manual scheduling I mean to say that you can do the same as with infect! for the other AgentEvent you have in your model.

It’s hard to follow the add_event! implementation without an example in the documentation. Do you mean to add add_event!(agent, model) into the other action functions? I have tried a few things and cant seem to make anything work with autogenerate_after_action = false

Perhaps I should just stick to the receiver model and give up on this one!

This is with autogenerate_after_action=false:

using Agents, Random, Graphs
using Random: Xoshiro
#rng = Xoshiro(42)

n_people = 1000

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

    # propensities
    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 == :S && (rand() <= susceptibility)
            contact.status = :E
            add_event!(contact, model)
        end
        add_event!(agent, model)
    end

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

    function become_infectious!(agent, model)
        agent.status = :I
        add_event!(agent, model)
    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
            add_event!(agent, model)
        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
            add_event!(agent, model)
    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, autogenerate_after_action=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 this is the same as with just the little change I suggested above, because I thought that

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

would have worked, but instead actually to make the model work you need to use:

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

which makes the model equivalent to the other one.

So I think that this is the right solution in the end. Out of curiosity (since I’m one of the designers of EventQueueABM), what do you plan to use it for? :slight_smile:

1 Like

Thanks so much for your help!

The infection peak for this model is still a bit higher than the receiver model though - receiver model get to about 50% infected around time step 50, whereas this model gets to about 70% infected. Do you know why that could be and if there is something else that could be added to the model?

I’m a PhD student and planning to use this to model infectious respiratory diseases in a population, particularly looking at what drives health inequities.