As I mentioned in my teaching blog, the population dynamics practical class in our Ecology subject examined the dynamics of COVID-19 this year, rather than doing an experiment with ciliates. For part of the prac, we simulated a basic SIR model (Susceptible, Infectious, Recovered) to illustrate the dynamics of an epidemic.

The model won’t represent the dynamics of infections of the coronavirus perfectly (e.g., the model assumes that infected individuals are immediately infectious, they remain equally infectious until recovery, and infected individuals are not quarantined), but it gives a reasonable sense of what might be expected under particular scenarios.

The situation in North Dakota, where rates of infection are at the highest per capita levels recorded for any state in the US got me thinking about what might be in prospect there. So let’s step through some modelling and see…

Currently (19 October) around 32000 people have been recorded as being infected (in a state of 762000 people, around 4%) and ~26000 have recovered, putting the recorded number of infected people at roughly 6000. Given the expanding nature of the epidemic, there are likely many more people currently infected (we’ll get to that later).

The *R* value is key to any SIR model. This parameter is the expected number of people that will be infected by each infectious person over the course of that person’s infection in the circumstance where the vast majority of the population is susceptible. It is possible to roughly estimate this value from rates at which new infections accumulate. I’ve estimated for North Dakota that it is roughly 1.7 at the moment. This amount averages over people who are identified as being infectious and are isolated and those who continue to spread the disease in the community.

If we assume that people remain infectious for 20 days on average, then the recovery rate parameter for the SIR model will be 0.05 (=1/20). Further, the daily infection rate will be *R*/20 (=0.085 if R = 1.7).

Those parameters then let us simulate SIR dynamics. We can start with 26000 recovered individuals. We know there are at least 6000 people currently infected, but the true number will be greater than that. If we assume the true number of infected people is currently 9000, then the simulated number of new infections per day matches the observed data, so we’ll run with that.

After 100 more days (roughly the end of January 2021), the SIR model predicts that the number of people infected per day will have passed its peak of 3500-4000 people. The number of people with infections would also be near its peak of 70000 people, which is roughly 10 times the current number. Bear in mind that North Dakota’s ICU units are already close to capacity – they would be completely swamped by a 10 fold increase in cases.

After 100 days, more than 250000 people are projected to have been infected (~35% of the population), and with an infection fatality rate of 1% we would expect 2500 to have died, with the number of deaths still increasing rapidly at that point.

These outcomes at day 100 are alarmingly bad, and the model is only projecting further misery from that point. Now, the projections might be overly pessimistic. Notably, mandates for restrictions might be brought in before these outcomes are achieved. However, North Dakota’s governor seems reluctant to mandate any restrictions.

In the absence of mandates, individuals are likely to change their behaviour. Currently, less than 5% of the population has been recorded as infected with the coronavirus. With only a small minority of those experiencing extremely poor health outcomes, people might not realise the challenge that lies ahead. However, as the proportion of people infected increases and the impacts of the disease on friends, family and acquaintances become more apparent, people might choose to take more precautions of their own. Changed behaviour, whether mandated or voluntary, will reduce the *R* value below current levels. That will reduce the number of people infected at anyone time and the rate at which cases accumulate.

Consider reducing R from 1.7 to 1.5. This small change (approximately 10% fewer people being infected by each infected person) highlights why small reductions in measures are important. For example, masks cannot prevent infections completely, but even small reductions are useful. Or consider if everyone reduced the number of interactions with other people in high transmission environments by even 10%, then the reduction in *R* would generate much better public health outcomes. When it comes to this pandemic, the 1-percenters really add up to being very useful.

With a stroke of luck, only a fraction of the population will be susceptible to infection. That is speculation at this stage, but we can use the model to answer what would happen in only 400000 people rather than 782000 were susceptible. That would reduce the expected number of deaths by the end of January to 1800 (instead of 2500) although deaths would still be increasing rapidly at that time. The number of infected people at any one time would peak at 35000 (instead of 70000). That is a better, rather hopeful scenario, but still bad.

If you want to run your own simulations, you are welcome to use the following R code. I’d be interested if you find it useful (or not):

# Simulate the SIR model in discrete time steps # Written by Michael McCarthy to approximate dynamics in North Dakota TotalN <- 762000 # Pop size of ND is estimated as 762062 Infected <- 9000 # The initial number of people infected Recovered <- 26000 # How many people have recovered so far - these are assumed # to be no longer infectious. Note, "recovered" people could be # anything between completely well to 100% dead. The main point # is that they are no longer susceptible or infectious. Susceptible <- TotalN - Infected - Recovered MaxT <- 200 # how many times steps for the simulation (days) R0 <- 1.7 # potential new infections per infected person over # the course of the infection gamma <- 0.05 # recovery rate per day of infected people beta <- R0 * gamma # beta = potential new infections per infected person per day for(t in 1:MaxT){ # for each time step (day) in the simulation NewlyI <- beta * Infected[t] * Susceptible[t] / TotalN # calculate the number of new infections NewlyR <- gamma * Infected[t] # calculate the number of new recoveries nextS <- Susceptible[t] - NewlyI # number of susceptible in the next time step (subtract new infections) nextI <- Infected[t] + NewlyI - NewlyR # number of infected in the next time step # (add new infections, subtract recoveries) nextR <- Recovered[t] + NewlyR # number of recovered in the next time step (add new recoveries) # join new amounts to the previous values, building vectors of values Susceptible <- c(Susceptible, nextS) Infected <- c(Infected, nextI) Recovered <- c(Recovered, nextR) } # Now Susceptible, Infected and Recovered will be vectors # with (t+1) values, each being the population size in each year. # Now plot the results of the simulation… # first set up the graphics to plot with smallish margins par(mfrow=c(1,1), mar=c(4,4,2,2)) # generate some values for the time axis (0, 1, 2, …, MaxT) Time <- 0:MaxT plot(Time, Susceptible, type="l", col="black", ylim=c(0, TotalN), ylab="Number") # first plot the susceptible population in black lines(Time, Infected, type="l", col="red") # add a line with the Infected number in red lines(Time, Recovered, type="l", col="blue") # add a line with the Recovered number in blue # Add some labels to the plot to help identify the results text(x=0, y=Susceptible[1]/1.1, labels="Susceptible", pos=4, col="black") text(x=MaxT, y=Recovered[MaxT]/1.1, labels="Recovered", pos=2, col="blue") text(x=MaxT, y=TotalN/20, labels="Infected", pos=2, col="red") plot(Time, Infected, type="l", col="red") # plot the number of people who are newly infected on each day NewPerDay <- Susceptible[1:(length(Susceptible)-1)] - Susceptible[2:length(Susceptible)] plot(ylab="New Cases per Day", xlab="Time", 1:MaxT, NewPerDay, type="l", col="orange")

Oh, and North Dakota is one of my favourite Lyle Lovett songs: