import numpy as np, scipy as sp, pandas as pd
from matplotlib import pyplot as plt
from sklearn.metrics.pairwise import euclidean_distances
from matplotlib import animation
from IPython.display import HTML
from itertools import combinations
from collections import namedtuple
"""
Units: time is measured in 12 hour chunks, but there's no strict correspondence to units here.
"""
class EffectiveArea:
"""Tells us the boundaries of the area that people can live in.
"""
def __init__(self,):
self.xmin, self.xmax = 0, 20
self.ymin, self.ymax = 0, 20
class Person:
"""
In this version, after 14 days, you either recover or die
"""
def __init__(self,effectivearea,state='succeptible',distancing=False):
self.r = 0.15
self.days_infected = 0
self._typical_recovery_days = 14
self.mortality_rate = {True:0.05,False:0.02} # index is whether we're over capacity
self.distancing = distancing
if self.distancing:
self.m = 1000
else:
self.m = 1
self.state = state
if self.state == 'infected':
self.days_infected = np.random.randint(self._typical_recovery_days)
self.ea = effectivearea
self.x, self.y = np.random.uniform(self.ea.xmin,self.ea.xmax), np.random.uniform(self.ea.ymin,self.ea.ymax)
if self.distancing:
self.distance(force=True)
else:
self.undistance(force=True)
def distance(self,force=False):
if force or not self.distancing:
self.m = 1000
self.vx, self.vy = 0,0
self.distance = True
def undistance(self,force=False):
if force or self.distancing:
self.m = 1
self.vx, self.vy = np.random.normal(size=2) # Maxwell Boltzmann?
self.distancing = False
def move(self,dt,ea,overcapacity=False):
self.x, self.y = self.x + self.vx*dt, self.y + self.vy*dt
"""People don't recover instantly. Once it's been two weeks,
they start getting a chance to recover. We could also give
them a chance to die here.
"""
if self.state == 'infected':
if self.days_infected < self._typical_recovery_days:
self.days_infected += dt
else:
if np.random.uniform() < self.mortality_rate[overcapacity]:
self.state = 'dead'
self.m = 1000
self.vx, self.vy = 0,0
else:
self.days_infected = 0
self.state = 'recovered'
def collide(p1,p2):
"""This is the 2D elastic collision problem from your intro physics book.
"""
if p1.state == 'infected' and p2.state == 'succeptible':
p2.state = 'infected'
elif p2.state == 'infected' and p1.state == 'succeptible':
p1.state = 'infected'
m1, m2 = p1.m, p2.m
r1, r2 = np.array([p1.x,p1.y]), np.array([p2.x,p2.y])
v1, v2 = np.array([p1.vx,p1.vy]), np.array([p2.vx,p2.vy])
M = m1 + m2
d = np.linalg.norm(r1 - r2)**2
u1 = v1 - 2*m2 / M * np.dot(v1-v2, r1-r2) / d * (r1 - r2)
u2 = v2 - 2*m1 / M * np.dot(v2-v1, r2-r1) / d * (r2 - r1)
p1.vx,p1.vy = u1
p2.vx,p2.vy = u2
class Universe:
def __init__(self,
npeople, # How many people there are in the world
initial_infection_chance=0.1, # Initial fraction of population which is infected
distancing=0.0, # Fraction of population which practices distancing
hospital_beds_fraction = 1.0 # Better is 0.003 which is closer to reality
):
self.npeople = npeople
self.ea = EffectiveArea()
self.dt = 0.1
self.hospital_beds_fraction = hospital_beds_fraction
self.hospital_beds = npeople*hospital_beds_fraction
self.data = None # gets set in self.run
def _state():
if np.random.uniform() < initial_infection_chance:
return 'infected'
return 'succeptible'
def _distancing():
return np.random.uniform() < distancing
self.people = [Person(self.ea,_state(),_distancing()) for i in range(self.npeople)]
# self.color = {'succeptible':0.5,'infected':0.0,'recovered':0.7} # old color scheme
self.color = {'succeptible':'lightblue','infected':'red','recovered':'green','dead':'midnightblue'}
def _step(self):
"""iterate through one timestep
"""
points = list(zip([p.x for p in self.people],[p.y for p in self.people]))
dists = euclidean_distances(points,points)
close = dists < 2*self.people[0].r
close = close.tolist()
for (i,j) in combinations(range(self.npeople),2):
if close[i][j]: # a bit faster than numpy indexing once things get big.
collide(self.people[i],self.people[j])
# Are we over capacity?
ninfected = len([p for p in self.people if p.state == 'infected'])
# about 5% need hospitalization
overcapacity = 0.05 * ninfected > self.hospital_beds
for p in self.people:
p.move(self.dt,self.ea,overcapacity)
if p.x <= self.ea.xmin or p.x >= self.ea.xmax:
p.vx = -p.vx
if p.y <= self.ea.ymin or p.y >= self.ea.ymax:
p.vy = -p.vy
def run(self,steps,stop_distancing_at = None):
"""Run a simulation
Internal data looks like
```
x_coords[frame,particle_number]
```
"""
x_coords = np.zeros((steps,len(self.people)))
y_coords = np.zeros((steps,len(self.people)))
state = np.zeros((steps,len(self.people)),dtype='object')
# SIR model. `i` is a questionable variable name TBH
s,i,r,d = np.zeros(steps),np.zeros(steps),np.zeros(steps),np.zeros(steps)
def pop_count(people):
s,i,r,d = 0,0,0,0
for p in people:
if p.state == 'succeptible':
s += 1
elif p.state == 'infected':
i += 1
elif p.state == 'recovered':
r += 1
elif p.state == 'dead':
d += 1
return s,i,r,d
s[0],i[0],r[0],d[0] = pop_count(self.people)
x_coords[0] = [p.x for p in self.people]
y_coords[0] = [p.y for p in self.people]
state[0] = [p.state for p in self.people]
for step in range(1,steps):
if step == stop_distancing_at:
for p in self.people:
p.undistance()
self._step()
s[step],i[step],r[step],d[step] = pop_count(self.people)
x_coords[step] = [p.x for p in self.people]
y_coords[step] = [p.y for p in self.people]
state[step] = [p.state for p in self.people]
dtype = namedtuple('RunData',['x','y','state','s','i','r','d','steps'])
self.data = dtype(x_coords,y_coords,state,s,i,r,d,steps)
def draw(self,ax=None):
"""
A very simple method to draw the current state. Better graphing comes in the
animation functions.
"""
if ax is None:
fig,ax = plt.subplots(figsize=(5,5))
scat = ax.scatter([p.x for p in self.people],[p.y for p in self.people],
c = [self.color[p.state] for p in self.people],
marker='.')
return scat,
def getanim(u,show_capacity=True):
fig,ax = plt.subplots(figsize=(6,6))
ax.axis('off')
left, width = 0.1, 0.85
bottom, height = 0.1, 0.65
spacing = 0.02
rect_universe = [left, bottom, width, height]
rect_trend = [left, bottom + height + spacing, width, 0.2]
ax_universe = plt.axes(rect_universe)
ax_trend = plt.axes(rect_trend)
#ax_universe.axis('off') # I want the box around it.
ax_universe.axes.get_xaxis().set_ticks([])
ax_universe.axes.get_yaxis().set_ticks([])
ax_trend.axis('off')
s,i,r,d = u.data.s,u.data.i,u.data.r,u.data.d
line, = draw_stacked_plot(u,ax_trend,show_time=True,show_capacity=show_capacity)
#ax_trend.stackplot(range(len(s)), i, s, r, d, labels=['sick','healthy','recovered','dead'],
# colors=[u.color['infected'],u.color['succeptible'],u.color['recovered'],u.color['dead']])
scat, = u.draw(ax_universe)
def drawframe(i):
data = np.column_stack(([u.data.x[i],u.data.y[i]]))
scat.set_offsets(data)
colors = np.array([u.color[_] for _ in u.data.state[i]])
scat.set_color(colors)
#_s,_i,_r = np.zeros(len(u.data.s)),np.zeros(len(u.data.s)),np.zeros(len(u.data.s))
#_s[:i],_i[:i],_r[:i] = u.data.s[:i],u.data.i[:i],u.data.r[:i]
#ax_trend.stackplot(range(len(_s)), _s, _i, _r, labels=['s','i','r'],colors=['blue','green','yellow'])
#ax_trend.legend(loc='upper left')
line.set_data([i,i],[0,u.npeople])
return scat,line
anim = animation.FuncAnimation(fig, drawframe, frames=u.data.steps,
interval=20, blit=False, repeat=False)
return anim
def draw_stacked_plot(u,ax=None,show_time=False,show_capacity=True,print_description=False):
if ax is None:
fig,ax = plt.subplots()
ax.axis('off')
s,i,r,d = u.data.s,u.data.i,u.data.r,u.data.d
if print_description:
print(f'{100*d[-1]/u.npeople}% of people died in the figure below')
ax.stackplot(range(len(s)), i, s, r, d, labels=['sick','healthy','recovered','dead'],
colors=[u.color['infected'],u.color['succeptible'],u.color['recovered'],u.color['dead']])
ax.legend(loc='center left')
if show_capacity:
ax.hlines(y=u.hospital_beds*20,xmin=0,xmax=len(s),linestyle='dashed')
if show_time is True:
line, = ax.plot([0,0],[0,s[0]+i[0]+r[0]+d[0]],linestyle='dashed',color='black')
return line,
return
Social Distancing and Hospital Capacity for coronavirus
There was a great Washington Post article on social distancing. It comes with simulations to show you the effect of different measures, and I think it’s fantastic #SciCom. I want to talk it over with my kids, but I want to be ready for them to ask questions like “what if we did (some cool other idea) for social distancing?” I also thought this would be useful for classes like intro physics, data science, computational modeling, etc., especially if those classes can team up with humanities, arts, etc. So, I want to have some Python code to replicate it and play with.
I put the Python code on github so that it’s easy to grab and play with. I like having notebooks as standalone things, so all of the relevant code from when I wrote this up is also at the end of the notebook. I tried to make the code basic enough that it’s easy to add different models yourself; please do so!
After poking at this, my main thought is that we’d be best off flattening the curve for as long as society can stand, while pushing hard at capacity and cures.
I’m no epidemiologist, but I think the results below make sense. The b
This is all based on a classic model called the SIR model, in which people can either be “succeptible” to getting the disease (S), “infected” (I), or “recovered” (R). Recovered people are considered immune from the disease
There are lots of modifications you can make to that model, and you can google for them.
Let’s start with a model like the one in the Washington Post article: a town of about 200 people that starts with one or two sick people. How long does it take for everyone to get infected? How long for everyone to recover?
Here’s a first pass at it. Each person is represented by a dot. The blue dots are healthy (S for “succeptible”) people. The red dots are sick people (I), and the green dots are recovered (R) people.
We’ll make one change from the Washington Post model right from the beginning: if you get sick, there’s a chance you’ll die. Figuring out that chance is tricky, but it seems like 2% is a good estimate. So the dots will turn dark blue when people die.
(In case it’s not obvious in your browser, that thing above is a movie. If you hover over it, you should have a button to play it. All of the plots that look like it below are also movies.)
You probably notice two things:
Basically everyone gets infected.
That red shape is the curve people talk about when they say #FlattenTheCurve.
It doesn’t look like a lot of people died.
You can barely see the dark blue bit showing up at the top of the plot. If you look at the dots, you might not even have noticed the five dots that turned dark blue and stopped moving. But those five dots represent 2.5% of the population. If we scale this model up to the US population of 330 million, 2.5% is 8.5 million people dead!
Luckily, this model probably doesn’t scale up perfectly to the whole US. But that’s a worst case scenario of “what if we did nothing” and it’s terrible.
Social Distancing
So what can we do? One idea is social distancing or social isolation. That’s the idea that you just stay at home and don’t interact with people. How much can that help? Let’s follow the Washington Post article and see what happens if 25% of people isolate.
That made almost no perceptible difference!
A bit about models like this: you should look at the general trends. Here, you can see that the shape of the red curve is about the same. The model uses a lot of randomness, so running it 100 times will give 100 different results. Looking for exact values of small numbers (e.g. the exact number of deaths) isn’t going to be as informative. For better estimates of that, you want to run this with a larger population, and you want to run it for longer. We’ll stick to qualitative results.
What if 50% of people distance?
Oof. Again, it’s about the same, though you can see some difference in that some people are not actually getting infected!
What if 75% of people distance?
This has got to be better, right?
That’s definitely better, but not perfect. What if we had really extreme distancing?
What if 90% of people distance?
That flattens the curve for a lot longer, and gives people a lot longer to deal with things.
But we’re also clearly not done with the pandemic, as there are still infected people in the population. What if we run that for longer?
So, if we could do an essentially perfect job of social distancing until the disease was eradicated, we win. That’s a tough sell, though, because I’m reading estimates that the time required for complete eradication could be much longer than 3-4 months. Like 12-18 months. And we haven’t even thought about hospital capacity.
What about hospital capacity?
What if we run out of hospital beds? In that case, obviously more people are going to die. The estimates of mortality rates when we have hospital capacity seem to range quite a bit, but 2% seems like a good estimate (I’m not linking to any sources here because things keep changing. Feel encouraged to download the code and re-run things yourself if you want a different estimate!).
The overall mortality rate without hospitals seems more like 5%. That’s not a factor of 10 or anything, but it’s quite a bit. And it becomes a bigger difference if hospitals get better at treating patients.
What’s our capacity? Googling around, it looks like we have a hospital capacity of around 3 beds per 1000 people in the US.
Bigger population: in order to see this, we’re going to want a bigger population! Both in terms of having a large enough population to see the capacity matter, and so differences of 1-2 people don’t make such a difference. That means things will take longer to compute. Drawing the awesome animation actually takes a long time, so I’m just going to show the charts you saw on top of those figures, pre-drawn, with no animations for the rest of this. Hopefully that makes sense at this point, since you’ve seen several animations.
So here’s the nobody-does-anything plot, original size. The dashed line at the bottom is the hospital capacity. In this case, here are the exact numbers: * with 3 beds per 1000 people, that’s 0.6 beds for our 200 people. * We expect 5% of sick people to need beds. That means that, out of 20 sick people, one of them is expected to need a bed. So, if we’re just talking about corona virus, we need one bed for every 20 sick people. One bed means we can handle 20 sick people. Two beds means we can handle 40 sick people. 0.6 beds means we can handle 20*0.6 = 12 sick people. So the dashed line in the chart below is at 12 people
NOTE if you want to know the percentage of people who died in the above simulation, you can see it in text right below the
Pythoncode/right above the figure. This one says “4% of people died in the figure below”. You’ll see a similar bit of text on all of the next figures, but it’s kind of hard to see, and I’m not great at text formatting.And now with 1000 people
You can see the awful trend there: the number of deaths starts going up dramatically once everyone is sick and the hospitals are over capacity. This only becomes more pronounced when hospitals get better at treating people!
Let’s say we start a bit later, with 5% of the population infected. That’s probably not a terrible guess. Given our lack of testing, it could even be higher than that.
That’s bad.
Does distancing help?
What if we have moderate social distance of 50%?
And what if we have extreme distancing of 90%?
And, just as a reminder, 5% of the US population is 16 million people. If pure social distancing saves 1%, that’s about 3 million people we would save. As a side note, I think the numbers from my simulations are higher than the ones I’m seeing because I’m modeling lots of worst-case scenarios.
What if we build more beds
We could certainly do something like mobilize the Army Corps of Engineers to build a bunch of hospitals. That’s what other countries did a while ago. What if we increased our hospital capacity by a factor of 5? I’m guessing that’s reasonable because we don’t need great multi-purpose hospitals. We need temporary stuff for one particular purpose. Maybe this includes the navy hospitals? Maybe it includes repurposed dorm rooms? I’m hoping. But I’m not a civil engineer; I could be wrong here.
So that’s way better. Not perfect, but way better.
What if people stop distancing?
With 5x the capacity
The above assumes we keep distancing until basically everyone’s recovered. What if we build these beds, and people stop distancing when it looks like the curve is going down?
So that’s bad, it’s still an improvement ofer not distancing, but it’s bad.
With 2x the capacity
What if we don’t get a 5x extra capacity? What if we do something (still seemingly awesome), like double our capacity instead?
Here’s with no distancing
Here’s with extreme distancing
And here’s with extreme distancing that stops early
Again, now we’ve flattened the curve, and we need to run for longer. You can’t quite tell, because I’m not showing time on the x-axis, but here’s running that for twice as long.
In the first of those simulations, we happened to stop distancing after we were actually under capacity. That seems to make a noticeable distance. But either way, it’s not great.
Conclusions
Can I make any conclusions from this? Even though this isn’t a quantitatively correct model, I think a few things are probably true:
Seems like we’d be best off flattening the curve for as long as society can stand, while pushing hard at capacity and cures.
Running it yourself
If you want to run it yourself, here’s the code. Instead of
up top, just run the cell below.
If you’re playing around with the code, here are some things to note: