# 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

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

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.