# 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:

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!

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

step | infected_status | recovered_status | |
---|---|---|---|

Int64 | Int64 | Int64 | |

1 | 1990 | 856 | 59 |

2 | 1991 | 858 | 57 |

3 | 1992 | 859 | 56 |

4 | 1993 | 860 | 55 |

5 | 1994 | 861 | 54 |

6 | 1995 | 861 | 54 |

7 | 1996 | 861 | 54 |

8 | 1997 | 863 | 52 |

9 | 1998 | 864 | 51 |

10 | 1999 | 865 | 50 |

11 | 2000 | 865 | 50 |

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
```