Continuous space social distancing for COVID-19

This is a model similar to our SIR model for the spread of COVID-19. But instead of having different cities, we let agents move in one continuous space and transfer the disease if they come into contact with one another. This model is partly inspired by this article, and can complement the SIR graph model. The graph model can model virus transfer between cities, whilst this model can be used to study what happens within a city.

The example here serves additionally as an introduction to using continuous space, modelling billiard-like collisions in that space, and animating the agent motion in the space. Notice that a detailed description of the basics of the model regarding disease spreading exists in the SIR example, and is not repeated here.

It is also available from the Models module as Models.social_distancing.

Moving agents in continuous space

Let us first create a simple model where balls move around in a continuous space. We need to create agents that comply with ContinuousSpace, i.e. they have a pos and vel fields, both of which are tuples of float numbers.

using Agents, Random

mutable struct Agent <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    mass::Float64
end

The mass field will come in handy later on, when we implement social isolation (i.e. that some agents don't move and can't be moved).

Let's also initialize a trivial model with continuous space

function ball_model(; speed = 0.002)
    space2d = ContinuousSpace((1, 1), 0.02)
    model = ABM(Agent, space2d, properties = Dict(:dt => 1.0), rng = MersenneTwister(42))

    # And add some agents to the model
    for ind in 1:500
        pos = Tuple(rand(model.rng, 2))
        vel = sincos(2π * rand(model.rng)) .* speed
        add_agent!(pos, model, vel, 1.0)
    end
    return model
end

model = ball_model()
AgentBasedModel with 500 agents of type Agent
 space: periodic continuous space with 50×50 divisions
 scheduler: fastest
 properties: Dict(:dt => 1.0)

We took advantage of the functionality of add_agent! that creates the agents automatically. For now all agents have the same absolute speed, and mass.

The agent step function for now is trivial. It is just move_agent! in continuous space

agent_step!(agent, model) = move_agent!(agent, model, model.dt)

dt is our time resolution, but we will talk about this more later! Cool, let's see now how this model evolves.

using InteractiveDynamics
import CairoMakie

abm_video(
    "socialdist1.mp4",
    model,
    agent_step!;
    title = "Ball Model",
    frames = 50,
    spf = 2,
    framerate = 25,
)

As you can see the agents move in a straight line in periodic space. There is no interaction yet. Let's change that.

Billiard-like interaction

We will model the agents as balls that collide with each other. To this end, we will use two functions from the continuous space API:

  1. interacting_pairs
  2. elastic_collision!

We want all agents to interact in one go, and we want to avoid double interactions (as instructed by interacting_pairs), so we define a model step and re-run the animation.

function model_step!(model)
    for (a1, a2) in interacting_pairs(model, 0.012, :nearest)
        elastic_collision!(a1, a2, :mass)
    end
end

model2 = ball_model()

abm_video(
    "socialdist2.mp4",
    model2,
    agent_step!,
    model_step!;
    title = "Billiard-like",
    frames = 50,
    spf = 2,
    framerate = 25,
)

Alright, this works great so far!

Agents.jl is not a billiards simulator!

Please understand that Agents.jl does not accurately simulate billiard systems. This is the job of Julia packages HardSphereDynamics.jl or DynamicalBilliards.jl. In Agents.jl we only provide an approximating function elastic_collision!. The accuracy of this simulation increases as the time resolution dt decreases, but even in the limit dt → 0 we still don't reach the accuracy of proper billiard packages.

Also notice that the plotted size of the circles representing agents is not deduced from the interaction_radius (as it should). We only eye-balled it to look similar enough.

Immovable agents

For the following social distancing example, it will become crucial that some agents don't move, and can't be moved (i.e. they stay "isolated"). This is very easy to do with the elastic_collision! function, we only have to make some agents have infinite mass

model3 = ball_model()

for id in 1:400
    agent = model3[id]
    agent.mass = Inf
    agent.vel = (0.0, 0.0)
end

let's animate this again

abm_video(
    "socialdist3.mp4",
    model3,
    agent_step!,
    model_step!;
    title = "Billiard-like with stationary agents",
    frames = 50,
    spf = 2,
    framerate = 25,
)

Adding Virus spread (SIR)

We now add more functionality to these agents, according to the SIR model (see previous example). They can be infected with a disease and transfer the disease to other agents around them.

mutable struct PoorSoul <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    mass::Float64
    days_infected::Int  # number of days since is infected
    status::Symbol  # :S, :I or :R
    β::Float64
end

Here β is the transmission probability, which we choose to make an agent parameter instead of a model parameter. It reflects the level of hygiene of an individual. In a realistic scenario, the actual virus transmission would depend on the β value of both agents, but we don't do that here for simplicity.

We also significantly modify the model creation, to have SIR-related parameters. Each step in the model corresponds to one hour.

const steps_per_day = 24

using DrWatson: @dict
function sir_initiation(;
    infection_period = 30 * steps_per_day,
    detection_time = 14 * steps_per_day,
    reinfection_probability = 0.05,
    isolated = 0.0, # in percentage
    interaction_radius = 0.012,
    dt = 1.0,
    speed = 0.002,
    death_rate = 0.044, # from website of WHO
    N = 1000,
    initial_infected = 5,
    seed = 42,
    βmin = 0.4,
    βmax = 0.8,
)

    properties = @dict(
        infection_period,
        reinfection_probability,
        detection_time,
        death_rate,
        interaction_radius,
        dt,
    )
    space = ContinuousSpace((1,1), 0.02)
    model = ABM(PoorSoul, space, properties = properties, rng = MersenneTwister(seed))

    # Add initial individuals
    for ind in 1:N
        pos = Tuple(rand(model.rng, 2))
        status = ind ≤ N - initial_infected ? :S : :I
        isisolated = ind ≤ isolated * N
        mass = isisolated ? Inf : 1.0
        vel = isisolated ? (0.0, 0.0) : sincos(2π * rand(model.rng)) .* speed

        # very high transmission probability
        # we are modelling close encounters after all
        β = (βmax - βmin) * rand(model.rng) + βmin
        add_agent!(pos, model, vel, mass, 0, status, β)
    end

    return model
end

Notice the constant steps_per_day, which approximates how many model steps correspond to one day (since the parameters we used in the previous graph SIR example were given in days).

To visualize this model, we will use black color for the susceptible, red for the infected infected and green for the recovered, leveraging InteractiveDynamics.abm_plot.

sir_model = sir_initiation()

sir_colors(a) = a.status == :S ? "#2b2b33" : a.status == :I ? "#bf2642" : "#338c54"

fig, abmstepper = abm_plot(sir_model; ac = sir_colors)
fig # display figure

We have increased the size of the model 10-fold (for more realistic further analysis)

To actually spread the virus, we modify the model_step! function, so that individuals have a probability to transmit the disease as they interact.

function transmit!(a1, a2, rp)
    # for transmission, only 1 can have the disease (otherwise nothing happens)
    count(a.status == :I for a in (a1, a2)) ≠ 1 && return
    infected, healthy = a1.status == :I ? (a1, a2) : (a2, a1)

    rand(model.rng) > infected.β && return

    if healthy.status == :R
        rand(model.rng) > rp && return
    end
    healthy.status = :I
end

function sir_model_step!(model)
    r = model.interaction_radius
    for (a1, a2) in interacting_pairs(model, r, :nearest)
        transmit!(a1, a2, model.reinfection_probability)
        elastic_collision!(a1, a2, :mass)
    end
end

Notice that it is not necessary that the transmission interaction radius is the same as the billiard-ball dynamics. We only have them the same here for convenience, but in a real model they will probably differ.

We also modify the agent_step! function, so that we keep track of how long the agent has been infected, and whether they have to die or not.

function sir_agent_step!(agent, model)
    move_agent!(agent, model, model.dt)
    update!(agent)
    recover_or_die!(agent, model)
end

update!(agent) = agent.status == :I && (agent.days_infected += 1)

function recover_or_die!(agent, model)
    if agent.days_infected ≥ model.infection_period
        if rand(model.rng) ≤ model.death_rate
            kill_agent!(agent, model)
        else
            agent.status = :R
            agent.days_infected = 0
        end
    end
end

Alright, now we can animate this process for default parameters

sir_model = sir_initiation()

abm_video(
    "socialdist4.mp4",
    sir_model,
    sir_agent_step!,
    sir_model_step!;
    title = "SIR model",
    frames = 100,
    ac = sir_colors,
    as = 10,
    spf = 1,
    framerate = 20,
)

Exponential spread

We can all agree that these animations look interesting, but let's do some actual analysis of this model. The quantity we wish to look at is the number of infected over time, so let's calculate this, similarly with the graph SIR model.

infected(x) = count(i == :I for i in x)
recovered(x) = count(i == :R for i in x)
adata = [(:status, infected), (:status, recovered)]

Let's do the following runs, with different parameters probabilities

r1, r2 = 0.04, 0.33
β1, β2 = 0.5, 0.1
sir_model1 = sir_initiation(reinfection_probability = r1, βmin = β1)
sir_model2 = sir_initiation(reinfection_probability = r2, βmin = β1)
sir_model3 = sir_initiation(reinfection_probability = r1, βmin = β2)

data1, _ = run!(sir_model1, sir_agent_step!, sir_model_step!, 2000; adata)
data2, _ = run!(sir_model2, sir_agent_step!, sir_model_step!, 2000; adata)
data3, _ = run!(sir_model3, sir_agent_step!, sir_model_step!, 2000; adata)

data1[(end-10):end, :]

11 rows × 3 columns

stepinfected_statusrecovered_status
Int64Int64Int64
1199085659
2199185857
3199285956
4199386055
5199486154
6199586154
7199686154
8199786352
9199886451
10199986550
11200086550

Now, we can plot the number of infected versus time

using CairoMakie
figure = Figure()
ax = figure[1, 1] = Axis(figure; ylabel = "Infected")
l1 = lines!(ax, data1[:, aggname(:status, infected)], color = :orange)
l2 = lines!(ax, data2[:, aggname(:status, infected)], color = :blue)
l3 = lines!(ax, data3[:, aggname(:status, infected)], color = :green)
figure[1, 2] =
    Legend(figure, [l1, l2, l3], ["r=$r1, beta=$β1", "r=$r2, beta=$β1", "r=$r1, beta=$β2"])
figure