Numerical simulations
The year is 2020. It’s a chilly December morning in San Diego, and you’re out of avocados for the fifth time that month. Ah well, guess it’s time for a grocery store run. You get into the car, and just before turning on the ignition, realize you’ve forgotten your mask on the dresser. With a slight groan, you make a dash for the warm comfort of your home. At least there’s some exercise involved right?
If you haven’t been living under a rock these past few months, you would be aware of the Covid-19 pandemic that has drastically changed the lives of millions around the globe. In addition to the economic, physical, and mental strains it has put on our lives, it has also become quite a political issue with many protesting against the numerous public health safety measures put in place. Wearing masks, social distancing, and proper, vigilante use of hand sanitizers and hand washing are some of the measures that have sparked debate in terms of their effectiveness in slowing the spread of this contagion.
With all these in mind it occurred to me; Could we simulate virus spread in a community in a way that allows us directly see the effects of these public health safety measures?
To answer that question, let’s have a brief understanding of simulations
Numerical simulations
A numerical simulation is essentially a computer program that mathematically models a physical system, allowing us to study the behavior of systems that are too complex to provide analytical solutions e.g non-linear systems. Some examples of these kinds of simulations include the finite difference method (FEM), finite element method, and agent-based Modeling. My focus in this article will be on agent-based modeling. Agent-based modeling is a type of simulation modeling technique with a large number of real-world applications such as traffic, customer flow management, stock market, operational risk, diffusion of innovation and adoption dynamics, etc.
In agent-based modeling, a system is modeled as a collection of autonomous decision-making entities (agents). These agents make decisions based on a set of rules and their individual situation and their repetitive interactions are the main features of this type of modeling. Monte Carlo methods are used to introduce randomness. Their ability to exhibit complex behavior patterns can provide valuable information about the dynamics of the real-world system it emulates.
Application to epidemiology
A particular area where the benefits of this type of modeling are evident is in simulating the spread of infectious disease. In a population of N individuals, an individual at any given time (t which is measured in days), can be in one of three categories:
- Infected and prone to infect others (infected) (I(t))
- Healthy but susceptible to infection (susceptible) (S(t))
- Recovered or removed (recovered) (R(t))
In other to simplify calculations, working with the fraction of the total population in each category (i.e divided by N) is used instead. Both sets of variables will give us the same information about the progress of the epidemic. A couple of assumptions are made in the modeling process;
- The only way an individual leaves the susceptible group is by becoming infected i.e The susceptible group isn’t influenced by external changes e.g births, immigration
- The rate of change of the number of susceptible individuals (dS/dt) is dependent on the original number of susceptible individuals (S(t))
- Under the assumption of a homogenously mixed population, each infected individual has a fixed number of contacts (b) per day with individuals in our system, and not all these contacts are susceptible. Therefore on average, each infected individual generates (b S(t)/N) new infected individuals per day
- A fixed fraction k of the infected group will recover during any given day.
We can illustrate how three categories of agents, in this case, will change over time in the diagram below:

The above assumptions are the basis of the classical SIR models of disease dynamics. This model is very straightforward and relies on the fact that when a susceptible and infected individual are in close proximity, the susceptible individuals contract the disease and become infected. Although reasonably predictive for infectious diseases, it is sometimes criticized for simplifying and unrealistic assumptions. Nevertheless, they are still useful in informing public health safety measures to mitigate and suppress outbreaks. For an in-depth look into the theory please visit [this link]
A simple python simulation
So can we make our own Simulation? The very first step involves defining the variables and rules governing our system. Specifically, we will be investigating the effects of following public safety guidelines in reducing the spread of COVID19. The basic guidelines are to
- Social distance (at least 6 feet apart)
- Wear a mask
- Wash your hands/hand sanitizer
Recently a vaccine was developed by Pfizer-BioNtech to combat the virus, so we will also investigate its effect on reducing the spread.
Simulation parameters
A square grid is used to represent the community containing the individuals (Agents). In our exam, a 100 x 100 square grid. We also assume:
- infection_radius (to reflect the 6 feet distance). which a
- recovery_threshold (to reflect recovering rate)
- death_threshold (to reflect death rate)
The idea is simple: At each iteration, the locations of the individuals in our system are shifted by a random amount shift. Then we check to see if an individual is within the radius of an infected individual (infection_radius). These individuals are then moved into the infected population. At every iteration, we can mirror (b S(t)/N) newly infected individuals per day, as well as k individual recoveries using Monte-Carlo methods. Samples the size of the infected population are drawn from a uniform distribution, and values that are higher and lower than our specified thresholds for recovery and deaths (recovery_threshold & death_threshold) are moved into the recovered and removed populations respectively. A uniform distribution is used here because all infected individuals are assumed to be equally likely to recover or die (a simplistic assumption I know). We also assume the recovered and dead population can not be infected again.
Public safety measures
We will treat public safety measures as entities we can introduce into our system. Specifically, we’ll consider the following variables:
- Mask
- SD
- Hygiene
- Vacc
Mask: If this entity is introduced, Mask is the factor by which the infection_radius drops i.e infection_radius X Mask if masks are used.
SD: This corresponds to the social distancing protocol and is the factor by which the shifts in individual locations at each iteration drops i.e shift/SD
Hygiene: The use of cleaning agents will reduce the transmissibility of the contagion. We can simulate this effect by randomly assigning an infectivity score from 1 to 0 to each infected individual, and only infected individuals with scores greater than this value can infect. It is an alternate way of saying a fraction of the susceptible population will not be infected.
Vacc: This corresponds to the use of a vaccine. If this variable is set in play, the assumption is the death rate will drop by a factor of Vacc i.e. death_threshold/Vacc.
Note: These values are relative values and not absolute. It’s usually best to run models relative to each other to quantify how changes in one parameter affect a system. You are free to use any values you like!
Implementation
Let’s first load in our libraries
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import style
style.use('dark_background')
from matplotlib import animation, rc
from IPython.display import HTML
rc('animation', html='html5')
We utilize the %matplotlib notebook command to allow our jupyter notebook to Livestream our plots dynamically, as opposed to statically plotting just a single figure. This is necessary in order to actually see our simulation. Now we define our parameters:
## Colors for the different states of the population
inf_col = 'red' # Infected
rem_col = 'lime' # Removed/reoovered
ded_col = 'blue' # Dead
sus_col = 'grey' # Susceptable
# Simulation parameters
infection_radius = 0.6 # Radius of infection (6 feet for covid-19)
recovery_threshold = (1-0.009975)# Assuming 99.75% recovery rate
death_threshold = 0.00021 # 2.1% death rate
stage_dimension = 100 (dimension of square grid)
# Initial population conditions
Pop = 800 # pop. size (per 100,000 of the population)
init = [Pop, # Total population
np.round(Pop*0.01), # infected (per 100 of the population)
0, # recovered population
0] # dead population
# Public safety variables
Mask_on=False
Mask=0.7
Hygiene_on = False
Hygiene = 0.9
SD_on=False
SD= 10
Vacc_on = False
Vacc = 2
We first define our colors to differentiate individuals/agents in the different population categories during the simulation. We then define the simulation parameters. The values I used were modeled based on values obtained at coronatracker at the time this article was written. We begin with a population of 800,0000 (scaled down to 800 for model purposes) and our assumption is 1% of the population has covid-19. The public safety variables are modeled to mimic strict policies i.e mask on, constant sanitization and disinfection, social distancing, and the use of vaccines.
Next, we define a class to model the system’s behavior: We begin by calculating how many individuals/agents belong in our susceptible population and then assigning colors to each point (individual) on our grid based on their population category. We then define the upper and lower bounds of the grid and generate random points within that grid from a uniform distribution. We then stream the data in a continuous loop and create the animation using the FuncAnimation method of the animation class. This method accepts a figure object (self.fig) upon which the plots are made, an initial function _initfunc (self.setup) which displays the locations of the individuals/agents before perturbation, and self.update which shows the population after perturbation and adjustment of population.
Now the most important bit, our data stream! So how does it look? well at each iteration we need to apply our random shifts to the population. If we have our social distancing flag on, this shift is reduced by a factor of SD. Then we isolate the infected population and calculate their distances to the susceptible population. We assign susceptible individuals whose distance from an infected individual is less than the infection_radius. Note that if the mask flag is on, this reduces the infection_radius by a factor Mask since the assumption is mask limits the spread of the disease. Next, we determine which infected individuals meet the criteria for being moved to either the recovered or removed population. If we apply the hygiene and vaccine flags, then only a certain fraction of infected individuals can infect to mimic fewer infections from using hand sanitizers and handwashing, and the death rate is decreased by a factor of Vacc to mimic the supposed drop in deaths due to vaccinations. To sum it all up, here is our simulation class defined below:
class Sim():
'''
The class Sim creates a simulation of the spread of COVID19 using user defined population and public health parameters.
'''
def __init__(self,init,Mask_on,Hygiene_on,SD_on,Vacc_on,gif_name):
self.N = init[0]
suc = self.N-init[1]-(init[2]+init[3])
self.cols = np.repeat([sus_col,inf_col,rem_col,ded_col],
[suc,init[1],init[2],init[3]])
self.low,self.high = -stage_dimension/2,stage_dimension/2
self.pos = np.random.uniform(self.low,self.high,size=(2,self.N))
self.stream = self.data_stream()
#print(self.removed)
self.fig,self.ax = plt.subplots();
self.ani = animation.FuncAnimation(self.fig, self.update, interval=200,
init_func=self.setup, blit=True, save_count=1)
## uncomment to save mp4 or gif. Increase save_counts for longer simulation file
#self.ani.save('simulation.mp4')
#self.ani.save(gif_name+'.gif', writer='imagemagick', fps=60)
HTML(self.ani.to_html5_video())
def setup(self):
'''
Function constructs frame based on the initial population state
'''
p,c= next(self.stream)
self.scat = self.ax.scatter(x=p[0,:],
y=p[1,:],
c=c.T)
return self.scat
def arg_within_radius(self,inf,susceptible,r):
'''
This function calculates distances between an infected and susceptible individuals
'''
if(Mask_on):
r=r*Mask
dist = np.sqrt(((inf-susceptible)**2).sum(axis=1))
return np.argwhere(dist<r).ravel()
def data_stream(self):
while True:
removed = 0
shift = 0.5*np.random.normal(0,1,size=(2,self.N))
if(SD_on):
shift=shift/SD
self.pos[0:2,:]+=shift
inf = np.argwhere(self.cols == inf_col).ravel()
inf_people = self.pos[:,inf].T
suc = np.argwhere(self.cols == sus_col).ravel()
suc_people = self.pos[:,suc].T
for i in inf_people:
if(Hygiene_on):
risk = np.random.rand(1)
if(risk > Hygiene):
infect = self.arg_within_radius(i,suc_people,infection_radius)
else:
pass
else:
infect = self.arg_within_radius(i,suc_people,infection_radius)
infect_index = suc[infect]
self.cols[infect_index] = inf_col
# adjust out of bound individuals
x1=np.where(self.pos[0:2,:]<self.low)
x2=np.where(self.pos[0:2,:]>self.high)
self.pos[0:2,:][x1],self.pos[0:2,:][x2] = self.low,self.high
r = np.random.uniform(0,1, size = inf.size)
rec = np.argwhere(r>recovery_threshold).T
recovered_idx = inf[rec]
self.cols[recovered_idx] = rem_col
d = np.random.uniform(0,1, size = (inf.size))
if (Vacc_on):
dead = np.argwhere(d<(death_threshold/Vacc)).T
else:
dead = np.argwhere(d<death_threshold).T
dead_idx = inf[dead]
self.cols[dead_idx] = ded_col
ded=np.argwhere(self.cols == ded_col).ravel()
ded_people = self.pos[:,dead].T
for dd in ded_people:
removed+=1
#self.removed+=removed
#print("An individual died")
yield self.pos,self.cols
def update(self,i):
'''
Function used in the animation function to plot population after pertubation and population assignment
'''
data,c= next(self.stream)
self.scat.set_offsets(data[0:2,:].T)
self.scat.set_sizes(np.zeros(self.N)+10)
self.scat.set_color(c.T)
return self.scat
Now with everything set, let’s run two simulations; one where we follow public health safety guidelines, and one where we don’t. Let’s begin with not following the public health guidelines.

Now, how about all measures in place?

All though these are snippets from the simulation before reaching steady-state, we notice that there’s only a single death in our scenario where all public health safety measures were followed, as opposed to two deaths in the initial scenario lacking any safety measures. In fact, if you run both simulations till steady-state, you’ll find the second scenario with safety measures had fewer deaths and a slower rate of spread of the contagion than the first scenario. If you don’t believe me, I invite you to run the code yourself.
Insights
As you would have expected, following the safety guidelines appears to have slowed the spread of the contagion and resulted in fewer deaths. But there are some caveats. This model is inherently simplistic; we not only used very simple assumptions and mathematical formulas, we assumed conditions were fixed for all individuals but in reality, there are other factors that must be taken into account; underlying conditions, length of exposure, etc. If we wanted a more complex model we would need to account for those. On the other hand, by keeping the conditions as simple as possible we can get some insights on a broad scale and by running different simulations where some parameters are varied, gain some understanding as to which variables have more of an impact. Overall simulation appears to be a very useful tool in making decisions. So with all that being said, the question of the day is:
Mask on or Mask off?

This project was mainly inspired by these resources ([[[here](https://www.linkedin.com/in/ibinabo-bestmann/)](https://github.com/Ibinaldo?tab=repositories)](https://www.coursera.org/learn/modeling-simulation-natural-processes), and here). The code will be provided here in the not so distant future, and my Linkedin profile can be found here. Thank you for reading! (and hopefully clapping!)
Ibinabo Bestmann