# Epidemic modelingLessons for the real world

In this section we'll use our models as a lens through which to examine some important questions for real-world epidemic planning. It's worth repeating the caveat that such model-based conclusions should always be held skeptically and weighed rigorously against broader real-world considerations. However, this exercise can nevertheless be useful, because it can help build intuition about the basic mechanisms underlying important phenomena which are present in the models and known by experts to be present in real-world epidemics as well.

What difference does it make if we start with a handful of infectious individuals rather than just one? The difference is profound and has significant implications for real-world infectious diseases.

Exercise
Modify the code for the function SIR_simulation below to set the number of initially infectious individuals to 10 rather than 1. What happens to the resulting histogram?

Use your observation to draw a conclusion about the strategy of temporary reduction of in the context of this model.

using Plots
import Base.Iterators: countfrom

function SIR_simulation(population_size, infection_probability)
statuses = fill("susceptible", population_size, 1)
statuses[1, 1] = "infectious"
transmissions = []
for t in countfrom(2)
n_infectious = sum(statuses[:, t-1] .== "infectious")
if n_infectious == 0
break
end
statuses = [statuses fill("susceptible", population_size)]
for k in 1:population_size
if statuses[k, t-1] == "recovered" || statuses[k, t-1] == "infectious"
statuses[k, t] = "recovered"
end
end
for j in 1:population_size
if statuses[j, t-1] == "infectious"
for k in 1:population_size
if statuses[k, t-1] == "susceptible"
if rand() < infection_probability
push!(transmissions, [(j, t-1),(k, t)])
statuses[k, t] = "infectious"
end
end
end
end
end
end
statuses, transmissions
end

population_size = 20
infection_probability = 2 / population_size
statuses, transmissions = SIR_simulation(population_size, infection_probability)
statuses

n = population_size = 1000
p = infection_probability = 2/population_size
n_recovered(n, p) = sum(SIR_simulation(n, p)[1][:, end] .== "recovered")
histogram([n_recovered(n, p) for _ in 1:5_000], xlims = (0, 1000),
nbins = 100, label = "", xlabel = "final number recovered",
ylabel = "number of runs")

Solution. What we see is that even with 10 initially infectious individuals, the disease spreads to a substantial percentage of the population with very high probability:

using Plots
import Base.Iterators: countfrom

function SIR_simulation(population_size,
infection_probability,
n_initially_infected)
statuses = fill("susceptible", population_size, 1)
statuses[1:n_initially_infected, 1] .= "infectious"
transmissions = []
for t in countfrom(2)
n_infectious = sum(statuses[:, t-1] .== "infectious")
if n_infectious == 0
break
end
statuses = [statuses fill("susceptible", population_size)]
for k in 1:population_size
if statuses[k, t-1] == "recovered" || statuses[k, t-1] == "infectious"
statuses[k, t] = "recovered"
end
end
for j in 1:population_size
if statuses[j, t-1] == "infectious"
for k in 1:population_size
if statuses[k, t-1] == "susceptible"
if rand() < infection_probability
push!(transmissions, [(j, t-1),(k, t)])
statuses[k, t] = "infectious"
end
end
end
end
end
end
statuses, transmissions
end

n = population_size = 1000
p = infection_probability = 2/population_size
n_recovered(n, p) = sum(SIR_simulation(n, p, 10)[1][:, end] .== "recovered")
histogram([n_recovered(n, p) for _ in 1:5_000], xlims = (0, 1000),
nbins = 100, label = "", xlabel = "final number recovered",
ylabel = "number of runs")

The final proportion infected depends on the value but not on the number of individuals who infectious initially.

While we derived this observation in the context of an SIR model, it does apply to real-life outbreaks as well: unless the disease is , matters far more than the number of individuals who are infectious to start.

When you heart the argument that an infectious disease should be ignored because the number of cases is still a small proportion of the population, it's analogous to arguing that an oven fire should be ignored until it at least engulfs the kitchen: flammability matters more than the size of the initial fire, unless the initial fire is nonexistent.

This observation also helps us see why a temporary reduction in does not prevent a spike in cases later once goes back up. Unless the conditions for a lower value are sustained long enough to eliminate the virus altogether, the post-reduction phase is essentially a fresh run of the model with large number of susceptible individuals, a high value, and enough initial infections to ensure widespread infection. To avoid a population-level outbreak, we need either a long-term reduction in or widespread immunity achieved through a vaccine.

## Containment vs Mitigation

The SIR model may be made more realistic in a wide variety of ways. One of the most glaring omissions is the lack of in the model: in the real world, the virus spreads by the kind of physical contact which is much more likely between people who live close to one another. Let's see what happens when we incorporate spatial relationships into the model.

Let's arrange our individuals into a square grid, and we'll let each person transmit to its four neighbors:

using Plots

function infection_plot(statuses)
p = plot(ratio = 1, legend = false, axis = false, grid = false)
for (status, color) in (("susceptible", "gray"),
("infectious", "tomato"),
("recovered", "darkgreen"))
pts = [Tuple(i) for i in CartesianIndices(statuses) if statuses[i] == status]
if length(pts) > 0
scatter!(p, pts, color = color, markersize = 6)
end
end
p
end

"Return the four nodes adjacent to (i, j)"
function neighbors(i, j, sidelength, barrier = Set())
filter(k -> all(1 .≤ k .≤ sidelength) && k ∉ barrier, [(i+1, j), (i-1, j), (i, j+1), (i, j-1)])
end

function spatial_simulation(sidelength, n_timesteps, R₀, barrier = Set())
statuses = fill("susceptible", sidelength, sidelength, n_timesteps)
statuses[1:3, 1:3, 1] .= "infectious"
for t in 2:n_timesteps
for i in 1:sidelength
for j in 1:sidelength
if statuses[i, j, t - 1] == "infectious"
nbs = neighbors(i, j, sidelength, barrier)
for (k,l) in nbs
if statuses[k, l, t-1] == "susceptible"
if rand() < R₀/length(nbs)
statuses[k, l, t] = "infectious"
end
end
end
statuses[i, j, t] = "recovered"
elseif statuses[i, j, t - 1] == "recovered"
statuses[i, j, t] = "recovered"
end
end
end
end
statuses
end

sidelength = 30
n_timesteps = 150
barrier_height = 25
barrier = Set([(10, k) for k in 1:barrier_height])
R₀ = 2.2
statuses = spatial_simulation(30, 150, R₀, barrier)
t = 150
infection_plot(statuses[:, :, t])
scatter!(collect(barrier), color = :black, markershape = :square)

Varying , we obtain the following animation:

Exercise
Experiment with other variations on this model to explore the following idea: barriers don't seem to help much unless they manage to completely contain the virus, because community spread quickly becomes the dominant mode of transmission.

This exercise is more open-ended than the others in this course. You might want to do this exploration in a Jupyter notebook.


Bruno