State Model

Let’s model this as a simple markov chain model - where you progress from healthy to infected to shedding (\(p_1\)) to symptomatic (\(p_2\)) to sick (\(p_3\)) to dead (\(p_4\)) and equally you can go from sheeding, symptomatic or sick to healthy (immune) (\(p_5\),\(p_6\),\(p_7\)). We can then let the number of shedding and symptomatic cases drive the number of infected (by some \(p_0\) value). We can let the cycle time be days and simply simulate.

Here we’re making the assumptions that sick (hospitalised) don’t drive infections and also that recovery directly from infected isn’t possible.

We can calculate the actual \(R_0\) as the number of infected from people that are shedding or symptomatic: \[ R_0 = \sum_{i=0}^{\infty}(1-p_2)^i p_0 + \sum_{i=0}^{\infty}\sum_{j=0}^{\infty}(1-p_2)^i (1-p_3)^j p_0 p_2 = \frac{p_0}{p_2}+\frac{p_0}{p_3} \]

This assume we’re ignoring recovery from shedding or symptomatic but we can easily augment to include this: \[ R_0 = \sum_{i=0}^{\infty}(1-p_2-p_5)^i p_0 + \sum_{i=0}^{\infty}\sum_{j=0}^{\infty}(1-p_2-p_5)^i (1-p_3-p_6)^j p_0 p_2 = \frac{p_0}{p_2+p_5}+\frac{p_0 p_2}{(p_2+p_5)(p_3+p_6)} \]

The transition probabilities can be modelled so these depend on the time the individual has been in a given state. Using an Erlang distribution seems like a natural choice for this. The probability density functions for various shape, \(k\) and rate, \(\lambda\) values are shown below.

One can easily formulate an agent based model which simply simulates the markov chain model with such distributions.

simulate <- function(N=1e6,R0=2.75,
                     k1=10,l1=4,
                     k2=3,l2=1.5,
                     k3=14,l3=2,
                     k4=4,l4=1,
                     k5=5,l5=0.5,w5=.85,
                     k6=10,l6=1,w6=.9,
                     k7=6,l7=1,w7=.95,
                     day.socdist=0,R0.socdist=0.5,
                     n.days=365,n.initial=10,random.seed=2) {
    set.seed(random.seed)
    x <- rep(0,N) # healthy
    t <- rep(1,N) # time in state
    x[1:n.initial] <- 1 # initial infected population
    t[1:n.initial] <- sample(7,n.initial,TRUE) # spread initial infection time
    n.infected <- rep(NA,n.days); 
    n.shedding <- rep(NA,length(n.infected)); 
    n.symptomatic <- rep(NA,length(n.infected)); 
    n.sick <- rep(NA,length(n.infected)); 
    n.dead <- rep(NA,length(n.infected)); 
    n.immune <- rep(NA,length(n.infected));
    tot.shedding <- 0
    tot.symptomatic <- 0
    tot.sick <- 0
    t.infection <- rep(NA,N)
    t.shedding <- rep(NA,N)
    t.symptomatic <- rep(NA,N)
    t.sick <- rep(NA,N)
    t.dead <- rep(NA,N)
    t.immune <- rep(NA,N)

    p0 <- R0/(sum(sapply(0:50,function(i) (1-(1-w5)*pgamma(i,k2,l2)-
                                           w5*pgamma(i,k5,l5))^i))+
              sum(sapply(0:50,function(i) (1-w5)*pgamma(i,k2,l2)*
                                          (1-(1-w5)*pgamma(i,k2,l2)-
                                           w5*pgamma(i,k5,l5))^i))*
              sum(sapply(0:50,function(j) (1-(1-w6)*pgamma(j,k3,l3)-
                                           w6*pgamma(j,k6,l6))^j)))
    
    for (day in 1:length(n.infected)) {
        if (day==day.socdist) p0 <- R0.socdist/(R0/p0)

        n.infected[day] <- length(i1<-which(x==1))
        n.shedding[day] <- length(i2<-which(x==2))
        n.symptomatic[day] <- length(i3<-which(x==3))
        n.sick[day] <- length(i4<-which(x==4))
        n.dead[day] <- sum(x==5)
        n.immune[day] <- sum(x==6)
        if (n.infected[day]+n.shedding[day]+n.symptomatic[day]+n.sick[day]==0)
            break
        ii.chtot <- c()
        if (length(i4)>0) {
            p1 <- (1-w7)*pgamma(1:max(t[i4]),k4,l4)
            p2 <- w7*pgamma(1:length(p1),k7,l7)
            for (k in 1:length(p1)) {
                iii <- i4[t[i4]==k]
                z <- sample(c(0,1,2),length(iii),
                            TRUE,c(1-p1[k]-p2[k],p1[k],p2[k]))
                x[j1<-iii[z==1]] <- 5 # dead
                x[j2<-iii[z==2]] <- 6 # immune
                t.dead[j1] <- day
                t.immune[j2] <- day
                ii.chtot <- c(ii.chtot,j1,j2)
            }
        }
        if (length(i3)>0) {
            p1 <- (1-w6)*pgamma(1:max(t[i3]),k3,l3)
            p2 <- w6*pgamma(1:length(p1),k6,l6)
            for (k in 1:length(p1)) {
                iii <- i3[t[i3]==k]
                z <- sample(c(0,1,2),length(iii),
                            TRUE,c(1-p1[k]-p2[k],p1[k],p2[k]))
                x[j1<-iii[z==1]] <- 4 # sick
                x[j2<-iii[z==2]] <- 6 # immune
                t.sick[j1] <- day
                t.immune[j2] <- day
                ii.chtot <- c(ii.chtot,j1,j2)
                tot.sick <- tot.sick+length(j1)
            }
        }
        if (length(i2)>0) {
            p1 <- (1-w5)*pgamma(1:max(t[i2]),k2,l2)
            p2 <- w5*pgamma(1:length(p1),k5,l5)
            for (k in 1:length(p1)) {
                iii <- i2[t[i2]==k]
                z <- sample(c(0,1,2),length(iii),
                            TRUE,c(1-p1[k]-p2[k],p1[k],p2[k]))
                x[j1<-iii[z==1]] <- 3 # symptomatic
                x[j2<-iii[z==2]] <- 6 # immune
                ii.chtot <- c(ii.chtot,j1,j2)
                t.symptomatic[j1] <- day
                t.immune[j2] <- day
                tot.symptomatic <- tot.symptomatic+length(j1)
            }
        }
        if (length(i1)>0) {
            p <- pgamma(1:max(t[i1]),k1,l1)
            for (k in 1:length(p)) {
                iii <- i1[t[i1]==k]
                z <- sample(c(0,1),length(iii),TRUE,c(1-p[k],p[k]))
                x[j<-iii[z==1]] <- 2 # shedding
                ii.chtot <- c(ii.chtot,j)
                t.shedding[j] <- day
                tot.shedding <- tot.shedding+length(j)
            }
        }
        i <- sample(N,trunc(sum(x==2 | x==3)*p0),TRUE)
        ii.chtot <- c(ii.chtot,i0<-i[msk<-x[i]==0])
        t.infection[i0] <- day
        x[i][msk] <- 1 # infecting uninfected
        t <- t+1
        t[ii.chtot] <- 1        
    }
    return(list(n.infected=n.infected[1:day],
                n.shedding=n.shedding[1:day],
                n.symptomatic=n.symptomatic[1:day],
                n.sick=n.sick[1:day],
                n.dead=n.dead[1:day],
                n.immune=n.immune[1:day],
                t.infection=t.infection,
                t.shedding=t.shedding,
                t.symptomatic=t.symptomatic,
                t.sick=t.sick,
                t.dead=t.dead,
                t.immune=t.immune,
                tot.infected=sum(x>0),
                tot.shedding=tot.shedding,
                tot.symptomatic=tot.symptomatic,
                tot.sick=tot.sick,
                tot.dead=sum(x==5),N=N))
}

Epidemiological and clinical observations

Overall we can from epidemiological and clinical observations attain that:

Based on the above we adopt \(R_0=2.75\), \(p_1 \sim E(10,4)\), \(p_2 \sim E(3,1.5)\), \(p_3 \sim E(14,2)\), and \(p_4 \sim E(4,1)\), where \(E\) denotes the Erlang distribution. We empirically consider the recovery distributions and fix these as \(p_5 \sim E(5,0.5)\), \(p_6 \sim E(10,1)\), and \(p_7 \sim E(6,1)\) with weights of 85%, 90% and 95% respectively.

Simulations

When we look at overall population dynamics the total population size becomes important, but for computational efficiencies we leave this at \(10^6\) unless otherwise stated.

s1 <- simulate(n.days=150)
simplot(s1,"Herd immunity approach")

Which provides a motality rate from those who become sick of 10.48 %, and of those symptomatic of 2.27 %, and a percentage of symptomatics becoming sick of 21.65 %, and a percentage of asymptomatics of 40.33 %. The overall population mortality rate is 1.34 %.

We can also investigate the distribution of progression from infected to symptomatic, sick or dead.

And we see this is reasonably consistent with the epidemiological and clinical observations.

Let’s compare with known data from the UK (circles) and Italy (triangles) where we adjust manually the time of observations to coincide with the simulations (\(-7\) for the UK and \(+8\) for Italy) based on day 1 in the observational data being 22 January 2020.

Which also seems reasonably consistent (at least until one starts to intervene by measures such as social distancing).

Interventions

Social distancing can be done at various efficiencies which we can simply model by adjusting \(R_0\) at the day, \(t\) when measures take effect (\(R_t\)). Based on the assumption that we need to see a reasonable number of fatalities before starting social distancing it seems reasonable to assume we start this around day 50 in practice (125 dead and 2094 sick). We can then investigate the effectiveness on different social distancing measures.

If we have perfect social distancing \(R_0=0\):

This leads to a mortality rate of 2.21 % for those becoming symptomatic, and of 0.55 % (5484 dead) for the overall population.

Assuming we have imperfect social distancing with \(R_0=0.25\):

This leads to a mortality rate of 2.27 % for those becoming symptomatic, and of 0.68 % (6811 dead) for the overall population.

As opposed to a more relaxed social distancing with \(R_0=0.5\):

Which leads to a mortality rate of 2.26 % for those becoming symptomatic, and of 0.83 % (8270 dead) for the overall population.

And let’s contrast this to a similar social distancing initiated earlier at day 40 (15 dead and 221 sick):

Which leads to a mortality rate of 2.25 % for those becoming symptomatic, and of 0.25 % (2544 dead) for the overall population.

Let’s generate an overall simulation showing the impact of time of social distancing with a given R_0. Note that we do this simulation with a reduced population size of \(10^5\) which implies that the timings may be optimistic (too short).

Considering the UK

Based on the above let’s now make a large simulation with a total population of 60MM (similar to the UK population size). We experienced lockdown on 23 March which seems to correspond approximately to day 55 in the simulations above (387 dead and 6002 sick (hospitalised)).

Initially we consider an \(R_0\) of 0.5:

Which leads to a mortality rate of 2.24 % for those becoming symptomatic, and of 0.15 % (89678 dead) for the overall population.

Zooming in on the first 100 days we have:

Suggesting that we will see the number of people hospitalised peaking about 2 weeks (day 71 ) after initiating social distancing but with a long tail.

Which we can contrast with a more effective one, using \(R_0=0.25\)

Which leads to a mortality rate of 2.23 % for those becoming symptomatic, and of 0.06 % (37594 dead) for the overall population. Also, in this case we have resolved the epidemic (no new infections, nor anyone able to infect) by day 189 , and the number of people sick peaks on day 70 .

Zooming in on the first 100 days we have:

And let’s now consider perfect social distancing

Which leads to a mortality rate of 2.25 % for those becoming symptomatic, and of 0.04 % (24992 dead) for the overall population. In this case we have resolved the epidemic by day 90 , and the number of people sick peaks on day 70 .

Zooming in on the first 90 days we have:

Of course, if we don’t practice social distancing at all, we will have infections spreading until we reach herd immunity with the consequence of predicted fatality rates as shown above of 1.34 % (this would not be notably affected by the population size):

Which, as expected, leads to a mortality rate of 2.23 % for those becoming symptomatic, and of 1.32 % (791036 dead) for the overall population. In this case we have resolved the epidemic by day 149 , and the number of people sick peaks on day 88 .

Finally, it may be worthwhile considering what happens if we practice social distancing which is still effective, but too poor, e.g. with \(R_0=0.75\):

Which leads to a mortality rate of 2.23 % for those becoming symptomatic, and of 0.56 % (336549 dead) for the overall population, and a much later peak for number of people sick on day 106 .

Zooming in on the first 100 days we have:

And one can observe the hump in the number of people sick.

Update in Response to Reviewer Feedback

The simulations provided above included data from Italy and the UK to March 28th when the initial simlations and the model above were performed. Below we’ve updated the simulations based on data available as of 31 August.

Total cases and deaths in the UK were obtained from https://coronavirus.data.gov.uk/. This data was norminally provided based on a start data of Jan-22. The UK lockdown began on 16 March 2020 when Matt Hancock told the House of Commons that all unnecessary social contact should cease. It has however not until 23 March that Boris Johnson told the country that people “must” stay at home and certain businesses must close.

Below the differences between the actual lockdown from a model perspective and window in which it actually happened is considered.

as.Date("2020-03-16")-(as.Date("2020-01-22")+6)
## Time difference of 48 days
as.Date("2020-03-23")-(as.Date("2020-01-22")+6)
## Time difference of 55 days
s2a <- simulate(N=6e7,day.socdist=50,R0.socdist=0.5,n.days=200,random.seed=3)
new.sick <- c()
for (i in 1:max(s2a$t.sick,na.rm=TRUE)) {
    new.sick <- c(new.sick,sum(s2a$t.sick==i,na.rm=TRUE))
}
plot(1:length(new.sick),cumsum(new.sick),log="y",ylim=c(1,1000000),type="l",col=2,xlab="Days",ylab="Cumulative Observations",
     main=expression("Relaxed social distancing at day 50 in the UK, R"[0]*"=0.5"))
lines((1:length(s2a$n.dead)),s2a$n.dead,col=1)
points((1:length(totalCases))-6,totalCases,col=2)
points((1:length(totalDead))-6,totalDead,col=1,pch=4)
legend("bottomright",fill=1:2,legend=c("Dead","Sick/Cases"),cex=0.5)

s2b <- simulate(N=6e7,day.socdist=48,R0.socdist=0.6,n.days=200,random.seed=3)
new.sick <- c()
for (i in 1:max(s2b$t.sick,na.rm=TRUE)) {
    new.sick <- c(new.sick,sum(s2b$t.sick==i,na.rm=TRUE))
}
plot(1:length(new.sick),cumsum(new.sick),log="y",ylim=c(1,1000000),type="l",col=2,xlab="Days",ylab="Cumulative Observations",
     main=expression("Relaxed social distancing at day 48 in the UK, R"[0]*"=0.6"))
lines((1:length(s2b$n.dead)),s2b$n.dead,col=1)
points((1:length(totalCases))-6,totalCases,col=2)
points((1:length(totalDead))-6,totalDead,col=1,pch=4)
legend("bottomright",fill=1:2,legend=c("Dead","Sick/Cases"),cex=0.5)

s2c <- simulate(N=6e7,day.socdist=55,R0.socdist=0.3,n.days=200,random.seed=3)
new.sick <- c()
for (i in 1:max(s2c$t.sick,na.rm=TRUE)) {
    new.sick <- c(new.sick,sum(s2c$t.sick==i,na.rm=TRUE))
}
plot(1:length(new.sick),cumsum(new.sick),log="y",ylim=c(1,1000000),type="l",col=2,xlab="Days",ylab="Cumulative Observations",
     main=expression("Relaxed social distancing at day 55 in the UK, R"[0]*"=0.3"))
lines((1:length(s2c$n.dead)),s2c$n.dead,col=1)
points((1:length(totalCases))-6,totalCases,col=2)
points((1:length(totalDead))-6,totalDead,col=1,pch=4)
legend("bottomright",fill=1:2,legend=c("Dead","Sick/Cases"),cex=0.5)

Which can be summarised in a single plot

This suggest that irrespective of the increased testing capacity in UK we still observe that the number of people sick rather than symptomatic track with the observed cumulative case numbers which may indicate a too relaxed view as what constitute sick vs symptomatic. We also observe some evidence that the actual lockdown in practice was initiated ahead of 23 March and closer to the initial messaging on initiating social distancing from Matt Hancock as the best post hoc fit seems to coinside with a lockdown around day 50 (corresponding to 18 March) as opposed to day 55 used above. We also adjust day 0 in the simulations to be 28 January instead of 29 January, as we now, with the increased set of data, can better estimate this.

Considering South Africa

As a comparator we consider South Africa which has a population size similar to that of the UK and is a low and middle income country (LMIC). South Africa followed a response similar to that of the UK with lockdown announced on 23 March 2020.

s3 <- simulate(N=6e7,day.socdist=50,R0.socdist=.85,n.days=220,w5=.95,w6=0.95,w7=0.98)

All in all, we observe that the increase in probabilities of the transition into the immune state together with the increase in \(R_0\) compared to the UK is sufficient to explain the much lower mortality rates and the much wider duration, and later peak, of reported cases compared to that observed in the UK. This simulation also suggest that about \(1/3\) of the population has been exposed to SARS-CoV-2, a much higher proportion that observed in the UK simulations.

Session Information

This document was genered using R markdown and the following session:

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 4.0.2 (2020-06-22)
##  os       macOS Catalina 10.15.6      
##  system   x86_64, darwin17.0          
##  ui       X11                         
##  language (EN)                        
##  collate  en_GB.UTF-8                 
##  ctype    en_GB.UTF-8                 
##  tz       Europe/London               
##  date     2020-09-06                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package     * version date       lib source        
##  assertthat    0.2.1   2019-03-21 [1] CRAN (R 4.0.0)
##  backports     1.1.8   2020-06-17 [1] CRAN (R 4.0.2)
##  callr         3.4.3   2020-03-28 [1] CRAN (R 4.0.0)
##  cli           2.0.2   2020-02-28 [1] CRAN (R 4.0.0)
##  codetools     0.2-16  2018-12-24 [1] CRAN (R 4.0.2)
##  crayon        1.3.4   2017-09-16 [1] CRAN (R 4.0.0)
##  desc          1.2.0   2018-05-01 [1] CRAN (R 4.0.0)
##  devtools      2.3.1   2020-07-21 [1] CRAN (R 4.0.2)
##  digest        0.6.25  2020-02-23 [1] CRAN (R 4.0.0)
##  ellipsis      0.3.1   2020-05-15 [1] CRAN (R 4.0.0)
##  evaluate      0.14    2019-05-28 [1] CRAN (R 4.0.0)
##  fansi         0.4.1   2020-01-08 [1] CRAN (R 4.0.0)
##  fs            1.4.2   2020-06-30 [1] CRAN (R 4.0.2)
##  glue          1.4.1   2020-05-13 [1] CRAN (R 4.0.0)
##  htmltools     0.5.0   2020-06-16 [1] CRAN (R 4.0.2)
##  knitr         1.29    2020-06-23 [1] CRAN (R 4.0.2)
##  magrittr      1.5     2014-11-22 [1] CRAN (R 4.0.0)
##  memoise       1.1.0   2017-04-21 [1] CRAN (R 4.0.0)
##  pkgbuild      1.1.0   2020-07-13 [1] CRAN (R 4.0.2)
##  pkgload       1.1.0   2020-05-29 [1] CRAN (R 4.0.0)
##  prettyunits   1.1.1   2020-01-24 [1] CRAN (R 4.0.0)
##  processx      3.4.3   2020-07-05 [1] CRAN (R 4.0.2)
##  ps            1.3.3   2020-05-08 [1] CRAN (R 4.0.0)
##  R6            2.4.1   2019-11-12 [1] CRAN (R 4.0.0)
##  remotes       2.2.0   2020-07-21 [1] CRAN (R 4.0.2)
##  rlang         0.4.7   2020-07-09 [1] CRAN (R 4.0.2)
##  rmarkdown     2.3     2020-06-18 [1] CRAN (R 4.0.2)
##  rprojroot     1.3-2   2018-01-03 [1] CRAN (R 4.0.0)
##  sessioninfo   1.1.1   2018-11-05 [1] CRAN (R 4.0.0)
##  stringi       1.4.6   2020-02-17 [1] CRAN (R 4.0.0)
##  stringr       1.4.0   2019-02-10 [1] CRAN (R 4.0.0)
##  testthat      2.3.2   2020-03-02 [1] CRAN (R 4.0.0)
##  usethis       1.6.1   2020-04-29 [1] CRAN (R 4.0.0)
##  withr         2.2.0   2020-04-20 [1] CRAN (R 4.0.0)
##  xfun          0.16    2020-07-24 [1] CRAN (R 4.0.2)
##  yaml          2.2.1   2020-02-01 [1] CRAN (R 4.0.0)
## 
## [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library

Day, Michael. 2020. “Covid-19: identifying and isolating asymptomatic people helped eliminate virus in Italian village.” BMJ, March, 1–1.

Lai, Chih-Cheng, Yen Hung Liu, Cheng-Yi Wang, Ya-Hui Wang, Shun-Chung Hsueh, Muh-Yen Yen, Wen-Chien Ko, and Po-Ren Hsueh. 2020. “Asymptomatic carrier state, acute respiratory disease, and pneumonia due to severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2): Facts and myths.” Journal of Microbiology, Immunology and Infection, March, 1–9.

Lauer, Stephen A, Kyra H Grantz, Qifang Bi, Forrest K Jones, Qulu Zheng, Hannah R Meredith, Andrew S Azman, Nicholas G Reich, and Justin Lessler. 2020. “The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application.” Annals of Internal Medicine, March, 1–7.

Li, Long quan, Tian Huang, Yong qing Wang, Zheng ping Wang, Yuan Liang, Tao bi Huang, Hui yun Zhang, Wei ming Sun, and Yu ping Wang. 2020. “2019 novel coronavirus patients’ clinical characteristics, discharge rate and fatality rate of meta-analysis.” Journal of Medical Virology, March, jmv.25757–7.

Linton, Natalie M, Tetsuro Kobayashi, Yichi Yang, Katsuma Hayashi, Andrei R Akhmetzhanov, Sung-mok Jung, Baoyin Yuan, Ryo Kinoshita, and Hiroshi Nishiura. 2020. “Incubation Period and Other Epidemiological Characteristics of 2019 Novel Coronavirus Infections with Right Truncation: A Statistical Analysis of Publicly Available Case Data.” Journal of Clinical Medicine 9 (2): 538–9.

Mizumoto, Kenji, Katsushi Kagaya, Alexander Zarebski, and Gerardo Chowell. 2020. “Estimating the asymptomatic proportion of coronavirus disease 2019 (COVID-19) cases on board the Diamond Princess cruise ship, Yokohama, Japan, 2020.” Euro Surveillance : Bulletin Europeen Sur Les Maladies Transmissibles = European Communicable Disease Bulletin 25 (10): 454.

Nishiura, Hiroshi, Tetsuro Kobayashi, Ayako Suzuki, Sung-mok Jung, Katsuma Hayashi, Ryo Kinoshita, Yichi Yang, et al. 2020. “Estimation of the asymptomatic ratio of novel coronavirus infections (COVID-19).” International Journal of Infectious Diseases, March, 1–9.

Rodriguez-Morales, Alfonso J, Jaime A Cardona-Ospina, Estefanía Gutiérrez-Ocampo, Rhuvi Villamizar-Peña, Yeimer Holguin-Rivera, Juan Pablo Escalera-Antezana, Lucia Elena Alvarado-Arnez, et al. 2020. Travel Medicine and Infectious Disease, March, 101623.

Singhal, Tanu. 2020. “A Review of Coronavirus Disease-2019 (COVID-19).” Indian Journal of Pediatrics 87 (4): 281–86.

Wang, Yixuan, Yuyi Wang, Yan Chen, and Qingsong Qin. 2020. “Unique epidemiological and clinical features of the emerging 2019 novel coronavirus pneumonia (COVID-19) implicate special control measures.” Journal of Medical Virology, March, jmv.25748–28.

Zhao, Shi, Qianyin Lin, Jinjun Ran, Salihu S Musa, Guangpu Yang, Weiming Wang, Yijun Lou, et al. 2020. “Preliminary estimation of the basic reproduction number of novel coronavirus (2019-nCoV) in China, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak.” International Journal of Infectious Diseases 92 (March): 214–17.