Virus Spread Simulation Revisited – Population Attributes and Healthcare Resources

Introduction

In my latest blog post (Simulating a Virus Spread – What you can do to help healthcare cope), I described the importance of social distancing in a pandemic in order to minimize the load on healthcare service, or as now is accepted as the concept of “flattening the curve”. I choose to revisit this post….from a hospital bed, being treated for….COVID-19.

I will make a quick recapitulation of the post, but to fully understand this blog post, I strongly urge you to read it through. Basically, I created a naive model describing the spread of a virus to monitor the proportion of individuals being susceptible to contract a virus, S(t), the proportion of individuals actually having been infected, I(t), and those that have recovered, R(t). I made it a point in the article that social distancing was a key ingredient to avoid flooding ICU:s with infected patients and thereby endangering patient safety. This is easily illustrated by the picture below:

spread0206 — The proportion of infected individuals in the left picture never reaches 20% of the entire population, while it is well above 50% in the right one. The difference is a direct consequence of a poor social distancing strategy. As you can see, the parameter beta is 0.2 in the left picture and 0.6 in the right one. Thus, the fewer contacts are, the less likely are infections.

As you can see, the three functions are dependent on time. As I pointed out, this model does not take into consideration several crucial factors to make it usable in real-world situations, such as the age distibution of the population, the influx of individuals or the age-dependance contact rate, $\beta$ . As was pointed out by S.Y.Del Valle, J.M.Hyman, H.W.Hethcote & S.G.Eubank in “Mixing patterns between age groups in social networks” (Social Networks, Volume 29, Issue 4, October 2007, pp. 539-554), the daily average number of people individuals come into contact with depends, among other things, on age, as can be seen in the graphic below:

Daily-average-number-of-contacts-per-person-in-age-group-j-The-average-number-of — Daily average number of contacts per person and age group. (ref. Article mentioned in the above paragraph)

What we then did in the previous post was to average out $\beta$ and determine S(t), I(t) and R(t) as the solutions of the system of differential equations (know as the SIR model and first expressed in this manner by by Kermack and McKendrick [Proc. R. Soc. A, 115, 772 (1927)]):

$\begin{aligned}\frac{\mathrm{d}S(t)}{\mathrm{d}t} &= -\beta \frac{S(t) \times I(t)}{N},\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t} &= \beta\frac{ S(t)\times I(t)}{N} - \gamma I(t),\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t} &= \gamma I(t). \end{aligned}$

Given some set of boundary conditions and where, as said above, $\beta$ is the contact rate and $1/\gamma$ is known as the recovery rate.

Intermediary States and Final State

In the simplistic model described above there where only three possibles states, namely the suceptible ones, the infected ones and the individuals that have recovered. It goes without saying that this is a rather unprecise image of reality. Indeed, there is a well known factor called “incubation period” that we haven’t taken into account in out. What actually happens is that an infected individual comes in contact with a susceptible one (healthy) who then contracts the virus. That person does not, however, immediately become infected, but instead enters an incubation phase. This phase ends with the emergence of the first symptoms. Now, every suceptible individual will be exposed in the same way as anybody would, hence, there is no variability of the exposed rate, call it $\lambda$ . This, however, does not mean that the impact of the virus on different individuals will be the same. Indeed, a multitude of records on the novel Corona have shown that the younger an infected individual is, the less likely he/she is to show symptoms.

So, taking this new Exposed state into consideration we need to rewrite our differential equations as

$\begin{aligned}\frac{\mathrm{d}S(t)}{\mathrm{d}t} &= -\beta \frac{S(t) \times I(t)}{N},\\ \frac{\mathrm{d}E(t)}{\mathrm{d}t}&= \beta \frac{S(t) \times I(t)}{N}-\lambda E(t),\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t} &= \lambda E(t) - \gamma I(t),\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t} &= \gamma I(t). \end{aligned}$

and depending on which initial conditions (or boundary conditions) we chose, the shapes of the S(t), I(t), E(t) and R(t) curves will look different. Let us start by choosing an adequate value for the incubation period. According the official estimates, the incubation period for the novel coronavirus COVID-19 is 2-14 days. To make life a little easier we will assume the mean incubation period length is 8 days. For the time being we will not make any difference between age groups, so looking at the contact rate $\beta$ , we’ll just assume that it averages around 16 people.

To make the exercise a little more realistic, let’s assume that we are looking at Stockholm and its population. Stockholm’s population in the first quarter of 2020 was measured to be 829 417 inhabitants. Let’s assume that the entire city had initially 50 infected individuals. This is not an unreasonable assumption given what we know happened at the beginning of the COVID-19 pandemic in Sweden as many returned to Stockholm from Italian ski-resorts. Given these assumptions and the observations above we can initiate our model in the following way

N = 829417 # Stockholm's population 2020
D = 6.0 
gamma = 1.0 / D
delta = 1.0 / 8
R_0 = 50.0 #The number of people and infected individual can infect 
beta = R_0 * gamma 
S0, E0, I0, R0 = N-1, 1, 0, 0 # Boundary conditions

Our initial system of differential equations that describe th SEIR-model can be expressed in Python in the following way:

def deriv(y, t, N, beta, gamma, delta):
S, E, I, R = y
dSdt = -beta * S * I / N
dEdt = beta * S * I / N - delta * E
dIdt = delta * E - gamma * I
dRdt = gamma * I
return dSdt, dEdt, dIdt, dRdt

and solved using

t = np.linspace(0, 100, 100)
y0 = S0, E0, I0, R0 
ret = odeint(deriv, y0, t, args=(N, beta, gamma, delta))
S, E, I, R = ret.T

where ‘odeint’ is the integration of ordinary differential equations. Plotting the resulting curves gives us the following picture:

SEIRModelSTHLM_no_dead — In this scenario we assumed that there were initially 50 infected individuals in Stockholm from day 1. Note that given the assumptions we made, it took barely 17 days before more than half the population became exposed and that y day 21 around a quarter of the cities population is infected.

The scenario described above shows clearly that very few infected individuals are needed for a medical situation to quickly become uncontrolled. Even if it is true that we have made no differentiation between age groups or contact rate and so forth, it barely took 20 days since the introduction of 50 infected individuals among some 830 000 (that is 0,006%) to reach some 250 000 infected individuals. This is without any possible doubt a situation that any healthcare system could not deal with. Of course, all infected do not require medical care, but still.

In the above model, we disregarded the fact that people actually can die of a virus infection such as COVID-19. So, let’s tweak our model a little to take this outcome into consideration. In the previous model, people Susceptible to be contract the virus would first be Exposed after having encountered an infected individual. When their first symptoms appear they transition from being exposed to being Infected and finally Recover from the virus. But in reality, the virus can, in some cases, turn out to be fatal. So there are actually to separate branches in the chain of events, being Susceptible -> Exposed -> Infected -> Recover (or Death). In mathematical terms this implies the following differential equations:

$\begin{aligned}\frac{\mathrm{d}S(t)}{\mathrm{d}t} &= -\beta \frac{S(t) \times I(t)}{N},\\ \frac{\mathrm{d}E(t)}{\mathrm{d}t}&= \beta \frac{S(t) \times I(t)}{N}-\lambda E(t),\\ \frac{\mathrm{d}I(t)}{\mathrm{d}t} &= \lambda E(t) - (1-\alpha) \gamma I(t)-\alpha \rho I(t),\\ \frac{\mathrm{d}R(t)}{\mathrm{d}t} &=(1-\alpha) \gamma I(t),\\ \frac{\mathrm{d}D(t)}{\mathrm{d}t} &=\alpha \rho I(t). \end{aligned}$

where $\alpha$ is the death rate and $\rho$ is the counterpart of the rate of recovery $\gamma$ for those that die. I will call it the non-survival rate. It is calculated in the following manner: If in average, individuals survive 10 days from the time they are infected until they die, then $\rho = 1/10$ .

Let us revisit the scenario above, only to introduce death as a second possible outcome of being infected. I am going to assume that the virus has a death mortality rate of 5% (disregarding age group specific mortality) and that those who die do so in 7 days.

N = 829417 # Stockholm's population 2020
D = 6.0 
gamma = 1.0 / D
delta = 1.0 / 8
R_0 = 50.0
beta = R_0 * gamma 
alpha = 0.05
rho = 1/7
S0, E0, I0, R0, D0 = N-1, 1, 0, 0, 0

The system of differential equations are given by:

def deriv(y, t, N, beta, gamma, delta, alpha, rho):
S, E, I, R, D = y
dSdt = -beta * S * I / N
dEdt = beta * S * I / N - delta * E
dIdt = delta * E - (1 - alpha) * gamma * I - alpha * rho * I
dRdt = (1 - alpha) * gamma * I
dDdt = alpha * rho * I
return dSdt, dEdt, dIdt, dRdt, dDdt

Solving the problem is again by integration:

t = np.linspace(0, 100, 100) # Grid of time points (in days)
y0 = S0, E0, I0, R0, D0 # Initial conditions vector
ret = odeint(deriv, y0, t, args=(N, beta, gamma, delta, alpha, rho))
S, E, I, R, D = ret.T

Plotting the result:

SEIRModelSTHLM_dead — In this scenario we assumed that there were initially 50 infected individuals in Stockholm from day 1. Note that given the assumptions we made, it took barely 17 days before more than half the population became exposed and that day 21, around a quarter of the cities population is infected. This scenario is different from the previous one in that not all individuals do recover from the virus, some die. In this case the probability of dying is 0.1.

It is noteworthy thing is that the virus spreads extremely quickly. Indeed, by day 30 there are no more susceptible individuals.

Indeed, the entire population is either in the infected group, the recovered and imune group or dead. In other words, herd immunity is reached very quickly. However, as you can see in the picture, this come at the price of approximately 50 000 deaths. This is something that no-one is prepared to pay, and rightly so. I am not going to enter any discussions on the different strategies that countries around the world have adopted because my opinion on the matter is irrelevant. I am neither an epidemiologist nor a virologist and cannot bring any new insight to the table.

At the end of my previous blog and at the beginning of this one, argued that many of the parameters I used for different models were averages for the entire population rather than groups specific values. I started for instance with the notion of contact rate, denoted $\beta$ and pointed out that individuals at different ages have different types of environment they interact with.

In the graph to the right, one can observe that children below the age of 5 interact with an average of 12 other individuals, which would be parents, siblings, grand-parents and some occasional nanny or babysitter. Past that age, theit contact increase due to school and in their twenties either studies or work implies a larger contact net. As Time passes by and we approach retirement our network of people shrinks.

It is an important observation in the context of pandemics explains why it is mostly relevant to minimize social contacts for individuals for ages between 20 and 60, since they are the one contributing most to the average contact rate. I will not deal with contact rate with a age group approach and will instead use this approach to investigate population distribution into age groups and age group specific mortality $\alpha$ . However, if we remember that R0 is the number of people an infected individual can infect, it should come as no surprise that R0 must be age dependent, at least according to the reasons I mentioned about how our human contact depend on age. It should therefor make sense to derive some function that emulates this. If you have read the previous blog on this subject, you should remember that I made the argument that reducing contacts (that is influencing the parameter $\beta$ , and thereby R0 since R0 = $\beta/ \gamma$ ) would have an enormous impact on the flattening of the curve. But, as I also pointed out, this approach is naive, because there is absolutely no way to implement this in pointwise manner. It is more likely to be a continuous process. Now, what we actually want to achieve is diminishing contacts between individual. The strategy one adopts determine how quickly this happens. The function to the left is a logistic function and is actually a rather good candidate for the type of function we need to consider. Of course, in this example it shows the opposite of what we would like to do. Indeed, the function displayed here is:

$f(x) = \frac{1}{1+\exp(-k(x-x_{0}))}$

What we need is a function that slopes the other way, that is that shows the number of social contacts BEFORE an intervention, R0_usual, and fewer AFTER, R0_intervened. How quickly this change occurs depends entirely on the strategy adopted and the length of time the strategy is put into action. This last observation ought to translate by how the slope of the logistic function looks like, i.e. the value of the parameter k. Going from no social distancing rules at day D to total lockdown day D+1 implies a very steep slope, i.e. a high value of k. We thus set our function to be:

$R0(t) = \frac{R0_{usual} - R0_{intervened}}{1+\exp(-k(t-t_{0}))} + R0_{intervened}$

Note that is $R0_{usual} = R0_{intervened}$ , then $R0(t) = R0_{intervened} = R0_{usual}$ , which means that nothing has been done. Let’s now see how what we have done previously can be transformed by using $R0(t)$ instead of the constant R0-value. We’ll conduct a few experiments with different parameter settings:

Early lockdown and extreme social distancing

In this experiment, presence of the virus is detected almost immediately after its arrival in Stockholm and the authorities choose to lockdown the city from have reduced social contacts from 50 to 2 in less than 10 days:

N = 829417 # Stockholm's population 2020
D = 6.0 
gamma = 1.0 / D
delta = 1.0 / 8
R0_usual, k, t0, R0_intervened = 25.0, 0.5, 10, 2

def logistic_R_0(t):
    return (R0_usual-R0_intervened) / (1 + np.exp(-k*(-t+t0))) + R0_intervened

def beta(t):
    return logistic_R_0(t) * gamma

alpha = 0.05
rho = 1/7
S0, E0, I0, R0, D0 = N-1, 1, 0, 0, 0

As you can see from the plot below, the number of fatalities is rather low compared to strategies involving no restrictions. The downside is the length of time under which such a strategy has to be withheld and the consequences that follow. Again, this article is not a forum for disscussing these issues and I am confident that our authorities do engage in these disscussions. Note, that herd immunity is reached already after 170 days. Can the restrictions be withdrawn by then?

EXPERIMENT1R0f — In this strategy, lockdown is done almost immediately after the discovery of the virus and extreme social distancing is held under the 400 days of the experiment.

Early lockdown, gradual adoption of extreme social distancing

In this experiment, we do exactly what was done in the previous experiment but transition from normal social contacts to extreme social distancing in a smoother fasion.

N = 829417 # Stockholm's population 2020
D = 6.0 
gamma = 1.0 / D
delta = 1.0 / 8
R0_usual, k, t0, R0_intervened = 25.0, 0.05, 10, 2
def logistic_R_0(t):
return (R0_usual-R0_intervened) / (1 + np.exp(-k*(-t+t0))) + R0_intervened
def beta(t):
return logistic_R_0(t) * gamma
alpha = 0.05
rho = 1/7
S0, E0, I0, R0, D0 = N-1, 1, 0, 0, 0

We solve the system of differential equation under the following initial conditions

t = np.linspace(0, 100, 100) 
y0 = S0, E0, I0, R0, D0 
ret = odeint(deriv, y0, t, args=(N, beta, gamma, delta, alpha, rho))
S, E, I, R, D = ret.T

and plot in

plotseird(t, S, E, I, R, D, L)

The slope, as you can see was changed from 0.5 to 0.05. This means that instead of going from 25 social contacts to 2 in 10 days, we allowed people to change their habits in 100 days. The result, in fatalities is a slight increase in the number of deaths but a huge increase in the number of infected individual under a very short period of time. Since we have no data on the individuals having been infected, it could be the difference a lot of people having light symptoms and people flooding ICU:s all over the city.

In this strategy, lockdown is done almost immediately after the discovery of the virus and extreme social distancing is gradually installed and held under 100 days.

Moderate approaches and moderate social distancing? Other attractive strategies

In these expertiments, the main goal has be total confinment or put in a more “romantic term”, Isolation, by useing an extreme measure restrictin social contacts with others to a bare minimum of two. It has worked, must I would be very surprised is psychologists world advice anyone to actually go through with such a strategy for en extended legth of time. Urge you as a reader to attemp the design of own strategies that you imagine would work. Think clearly about your agendas. Do you want to save lives? If so, which ones? Currently ill people or future illnesses? Do you want heard immunity? If so, at was cost?

An important aspect which I would love to devoted an entire post sbout are the social, moral and ethical aspects of the entire problem.It is no longer something that ponder upon anymore, namely, The problem of Pain. You do no need to have ant spiritual backround for this question to have relevance in this context. Many, Among which C. S. Lewis, Have spend a life time trying to answer this question and I must say, it even ha it’s place here. I’ll leave you to that excercise.

I will finalize this already incredibly length post by looking att how to deal with age groups.

Age group dependence of death rate

No one can have missed the fact that most individual that die of COVID-19 in IUC rooms are predomidantly elderly people. There are many reasons for this: with age often comes poorer overall health, but also an increased probability for comorbidities. It is therefore resonable to think that looking att age-dependent mortality $\alpha(t) = \alpha_{t}$ , where $t$ is counted in years. Doing us will give us a better picture of how the virus impacts society. Until now, our model has assumed that everyone in the population had the same risk of dying, when all the empiric data we’ve seen so far tell us this is not true. Since all the examples above related to Stockholm, I chose to look att the official statistics for the city on SCB (Statistics Sweden or in Swedish Statistiska Centralbyrån), collected and organized the following data:

TABLESWE — Data collected from SCB about the population distributiion in Stockholm in 10-year classes. The bottom row in the table gives the percentage of the entire population is represented by the age group.

At the same time, mortality in COVID-19 for different age groups has been estimated to by different institutions

MORTSLITYCOVID — Age group mortality rate of Covid 19. Note that these numbers have been gathered from WHO om march 1. The number may, since then, have been revised.

Let’s try to apply this to our model. Until now, in every experiment we conducted, we gave an average mortality for the entire population

alpha = 0.05

What we now need to do is to replace this by 1) Giving information about the population distribution over age groups as well as 2) the age group specific mortality. Unfortunately, the population distribution is not given with the same age group disctinctions as mortality, but for this excercise, we’ll do some approximations. One very important thing to notice and that has been extensively discussed is the impact that a large amount of infected individuals would have on healthcare services ability to cope. Too many infected in an old population with severe risk factors imply a great number of deaths. So, to take this fact into account, I have introduced a factor, IMPACT, that will influence mortality if there are many infected individuals. This factors “efficiency” will be greater, the older the population is and the more individuals are infected.

agegroup_specific_mortality = {"0-29": 0.17, "30-39": 0.22, "40-49": 0.20, "50-59": 0.14, "60-69": 0.13, "70-79": 0.11 , "80-89":0.025, "90.":0.005}
pop_size_by_agegroup = {"0-29": 0.0, "30-39": 0.003, "40-49": 0.004, "50-59": 0.01, "60-69": 0.035, "70-79": 0.125 , "80-89":0.197, "90.":0.227}
IMPACT = 1
alphaPop = sum(agegroup_specific_mortality[i] * pop_size_by_agegroup[i] for i in list(agegroup_specific_mortality.keys()))

Lets, see how this all plays out. This is in the setting of a forced social distancing strategy with a rather smooth transistion (from 25 contact to 2 in a period of 100 days).

N = 829417 # Stockholm's population 2020
D = 6.0
gamma = 1.0 / D
delta = 1.0 / 8
R0_usual, k, t0, R0_intervened = 25.0, 0.05, 10, 2

def logistic_R_0(t):
return (R0_usual-R0_intervened) / (1 + np.exp(-k*(-t+t0))) + R0_intervened

def beta(t):
return logistic_R_0(t) * gamma

agegroup_specific_mortality = {"0-29": 0.17, "30-39": 0.22, "40-49": 0.20, "50-59": 0.14, "60-69": 0.13, "70-79": 0.11 , "80-89":0.025, "90.":0.005}
pop_size_by_agegroup = {"0-29": 0.0, "30-39": 0.003, "40-49": 0.004, "50-59": 0.01, "60-69": 0.035, "70-79": 0.125 , "80-89":0.197, "90.":0.227}
IMPACT = 1
alphaPop = sum(agegroup_specific_mortality[i] * pop_size_by_agegroup[i] for i in list(agegroup_specific_mortality.keys()))
rho = 1/9 
S0, E0, I0, R0, D0 = N-1, 1, 0, 0, 0

t = np.linspace(0, 100, 100) 
y0 = S0, E0, I0, R0, D0

ret = odeint(deriv, y0, t, args=(N, beta, gamma, delta, alpha_opt, rho))
S, E, I, R, D = ret.T
R0_over_time = [logistic_R_0(i) for i in range(len(t))] 
Alpha = [IMPACT * I[i]/N + alpha_opt for i in range(len(t))]

The result is over 50 000 dead and the total number of infected reaches almost 200, which is not a very successful strategy. The overall mortality rate is over 23% on day 60.

Let’s redo this experiment and change the strategy. We’ll introduce an early lockdown with rapid restrict measures and assume a low impact on healthcare services.

N = 829417 # Stockholm's population 2020
D = 6.0
gamma = 1.0 / D
delta = 1.0 / 8
R0_usual, k, t0, R0_intervened = 25.0, 1, 10, 2

def logistic_R_0(t):
return (R0_usual-R0_intervened) / (1 + np.exp(-k*(-t+t0))) + R0_intervened

def beta(t):
return logistic_R_0(t) * gamma

agegroup_specific_mortality = {"0-29": 0.17, "30-39": 0.22, "40-49": 0.20, "50-59": 0.14, "60-69": 0.13, "70-79": 0.11 , "80-89":0.025, "90.":0.005}
pop_size_by_agegroup = {"0-29": 0.0, "30-39": 0.003, "40-49": 0.004, "50-59": 0.01, "60-69": 0.035, "70-79": 0.125 , "80-89":0.197, "90.":0.227}
IMPACT = 0.01
alphaPop = sum(agegroup_specific_mortality[i] * pop_size_by_agegroup[i] for i in list(agegroup_specific_mortality.keys()))
rho = 1/9 
S0, E0, I0, R0, D0 = N-1, 1, 0, 0, 0

t = np.linspace(0, 400, 100) 
y0 = S0, E0, I0, R0, D0

ret = odeint(deriv, y0, t, args=(N, beta, gamma, delta, alphaPop, rho))
S, E, I, R, D = ret.T
R0TimeLine = [logistic_R_0(i) for i in range(len(t))] 
Alpha = [IMPACT * I[i]/N + alphaPop for i in range(len(t))]

The result is around 15 000 deaths and a maximal death rate of 2,79%. This is due to the fact that healthcare services were given the chance to cope with the influx of patients.

There are many experiments that couls be rum and I’ll let you play around with this at you leisure. As for me, time to rest. I am writting this post from a hospital bed in Stockholm being treated for…..yes, COVID-19. Data Science never stops!

Introduction

Intermediary States and Final State

Early lockdown and extreme social distancing

Moderate approaches and moderate social distancing? Other attractive strategies

Age group dependence of death rate

Dela det här:

Like this:

Related

Leave a Reply Cancel reply