Intro to Infection Modelling

Like everyone else, the COVID-19 pandemic has drastically changed how our group operates and interacts. However, we've been luckier than others in our operations since relatively few projects that are currently active require access to the university experimental facilities which are now closed. 

Everyone has a different way of coping with the anxiety and stress that an unprecedented event like this one produces. For many scientists, there is a certain comfort in information. Sometimes the act of defining, modeling, classifying can make you feel as if you have more control over an uncontrollable situation and, at least, you might gain insight that might be actionable. 

The next two articles on this blog, submitted by the Postdocs, discuss the pandemic. The first, below, tries to consider the modeling of these systems. The second tries to grapple with the things that scientists get "wrong," by laying bare how piecing together a real-world puzzle from incomplete information can sometimes lead to dramatic shifts in our understanding of how these systems operate. In times like these, the public often turns to scientists for advice and it's important to remember that none of us has a crystal ball.

Wherever you are, I hope you are as safe and secure as you can be and are managing the physical and mental toll. No matter your circumstances you aren't alone in this.

by Dr. Paul Godin

With the COVID-19 pandemic drastically impacting our day-to-day lives, I thought it worth examining how epidemiologist model infections and how social distancing can “flatten the curve”.  (Note, this idea was inspired by Peter Taylor, who provided the initial matlab code to run the simulations)

To do this we’ll start with the SEIR model (https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SEIR_model): Susceptible>Exposed>Infectious>Recovered/Dead. This model categorizes a population into 5 different categories:

1.     Susceptible: Persons who have not been exposed to the virus
and therefore aren’t carriers nor are they immune.

2.     Exposed: Persons who carry the virus but are not exhibiting symptoms. They can infect the susceptible population. They may also recover or move to exhibiting symptoms.

3.     Infectious: Persons who carry the virus and are exhibiting symptoms.
They can infect the susceptible population. They may also recover or die from the virus.

4.     Recovered: Persons who were either exposed or infected but have recovered and are now immune. They no longer carry the virus.

5.     Dead: Infectious persons who have died from the virus.

An additional assumption is that the total population is constant, i.e. no births, no deaths other than from the virus, and no migration. The total population in our model can be expressed as the sum of the 5 categories: P=S+E+R+I+R+D. Going forward we’ll normalize all our values relative to the total population.

Population modelling is a simple matter of setting up equations that relate how a population category changes with time, expressed as coupled differential equations. Let’s start with the exposed population; susceptible people become exposed after coming in contact with other exposed individuals or infected individuals. This can be expressed as:

where b_e and b_i are the rates at which an interaction between a susceptible person with exposed or infected results in an exposed person. But individuals can also move from exposed into either the recovered or infected groups, so we’ll add terms to the above equation to express those changes:

where c is the rate at which exposed individuals start expressing symptoms and are classified as infected, and k_e is the rate at which exposed individuals recover and develop an immunity. The first to terms depend on an interaction between to different groups, hence why the populations of those groups are included in those terms. How many individuals recover or start showing symptoms is a fraction of the existing exposed individuals, so those terms only depend on the current number of exposed individuals.

Let’s derive rates of change for the other groups. Infected individuals only come from the exposed population that develop symptoms. Infected individuals are lost when they either die or recover. Thus, the change of infected population with time is:

where k_i and m_i are the rates at which infected individuals recover and die respectively. The rate of the recovered and death populations are straight forward from the above equations:

To solve the above equations to express the populations as a function of time not possible to do analytically, so we’ll approximate the solution numerically. For a given population category, we can approximate the population at a later time as:

where P is a population category (S,E,I,R,D), t is time, and delta_t is a time-step. The smaller we make the time-step the more accurate our model becomes. These numerical solutions require an initial population value, and we can repeatedly apply the above equations at each time step to produce a full plot of the populations as a function of time. To start we’ll assume that there are no recovered, dead, or infected individuals in our total population. We’ll have a tiny fraction of individuals who are exposed, and the rest are susceptible (. Those of you paying attention may have noticed we’re missing a rate of change for the susceptible population, but we don’t actually need one, since we can express the susceptible population as everybody who isn’t in one of the other categories: S=1-E-I-R-D.

The last thing we need to determine the values of the rate constants in the above equations. This is where the expert epidemiologist would come in and talk about how they determine these numbers. But since I’m a physicist were just going to make up some initial values. To start lets use b_e=0.5, b_i=0.5, k_i=0.1, k_e=0.2, m_i=0.004, and c=0.2. With these numbers the rates of susceptible individuals becoming exposed is the same for interactions with exposed or infected individuals. Exposed individuals are more likely to recover than infected individuals. Exposed individuals are equally likely to recover as they are to develop full symptoms.

Let’s see what these model values look like:

With these values we see a sharp peak of infected/exposed individuals quickly into our simulation. Most people will ultimately become exposed and develop immunity, but a small significant fraction of the population did die. 

Now let’s quarantine infected individuals. To simulate this, we’re going to lower the infection rate b_i to 0.1, or to put it in words it’s less likely for known infected individuals to interact with the susceptible population. Let’s look at the model now:

With quarantines in place, the peak of infections occurs a few days later, and is more spread out (flattening the curve). Less individuals overall have been exposed and less have died. This is great! Quarantining has help!

Lastly, let’s add social distancing, we reduce b_e to 0.3, to simulate avoiding social interactions, but we’ll not drop it as much as we did b_i since people are still allowed to run errands and some people are still working essential services. Here’s what adding social distancing does to our model: 

Look at how flat that infection curve is! So, when you hear on the news experts saying we need social distancing this is what they are referring too. Of course, this is a simple model with quite a few assumptions in it. For example, there have been reported cases of recovered individuals testing positive again for COVID-19, so our assumption of complete immunity is not correct. Below is the code I used to create the plots using matlab (original source code from Peter Taylor at York University). Feel free to play around with the rate constants or add additional couplings between the groups if you want to try and make a more accurate model.

%Code to Run SEIR Model

%Time length and Time Step

t=0:0.05:500;

dt=0.05;

%Initial Population: Mostly susceptible, small fraction exposed, rest are 0

e = 0.00001; %initial exposed

in = 0.0; %initial infected

s = 1-e; %initial susceptible

r = 0; %initial recovered (immune)

d=0; %initial deaths

S(1)= s;IN(1)=in; R(1)= r; E(1)=e; %populations as a function of time

be=0.5; %exposed from exposed-susceptible interaction constant

bi=0.5; %exposed from infected-susceptible interaction constant

c=0.2; %infected from exposed constant

ke=0.2; %recover from exposed constant

ki=0.1; %recover from infected constant

mi=0.004; %death from infection constant

%calculate population fractions at each time-step

for j = 2:10001

%calculate current rates of change for given populations

dedt=bi*in*s +be*e*s-(c+ke)*e; %rate of change of exposed

dindt=c*e-(ki+mi)*in; %rate of change of infected

drdt=ki*in+ke*e; %rate of change of recovered

dddt=mi*in; %rate of change of deaths

%update populations for current time step

e=e+dedt*dt;

in = in + dindt*dt;

r = r + drdt*dt;

d = d+dddt*dt;

s = 1-e-in-r-d;

E(j)=e; IN(j)=in; R(j)=r; S(j)=s; D(j)= d;

end

%plot population projections

figure

plot(t,S,'blue')

xlim([0 400]); ylim([0 1]); xlabel('Time (Days)')

ylabel('Fraction of Population')

hold on; plot(t,E,'m');plot(t,IN,'r'); plot(t,R,'g');

plot(t,D,'black');

legend('Susceptible','Exposed','Infected','Recovered','Deceased')