From 6ca605bc53cf62d204c50cc3d22033f319976128 Mon Sep 17 00:00:00 2001 From: Gimenez Date: Mon, 30 Mar 2020 15:12:29 +0200 Subject: [PATCH] Add files via upload --- sir_covid19.html | 94 +++++++++++++++++++++++++++++++++++++----------- 1 file changed, 73 insertions(+), 21 deletions(-) diff --git a/sir_covid19.html b/sir_covid19.html index 0ac8390..7b8a455 100644 --- a/sir_covid19.html +++ b/sir_covid19.html @@ -369,13 +369,13 @@

SIR covid-19

Motivation

Last week, the Montpellier Evolution and Ecology seminar was broadcasted online. The organisers had the awesome idea of inviting experts in epidemiology and evolution to talk about the covid-19. Mircea Sofonea talked about estimating epidemiological parameters of the virus spread, Nicolas Bierne talked about the genomic epidemiology of SARS-CoV-2 and Sylvain Gandon talked about epidemiology and evolution of SARS-CoV2.

-

The talks were recorded, and the videos can be found on Youtube here.

-

I am no expert in genomics, and Mircea, Samuel Alizon and their colleagues have made their work available in several reports that include R codes. Here, I will focus my attention on Sylvain’s work with Sébastien Lion and Troy Day. I am gonna reproduce their analyses in R. In brief, Sylvain Gandon and colleagues use a SIR model to illustrate the effect of different strategies of social distancing on the epidemic. Sylvain’s talk starts at the 49th minute here and lasts approx. 25 minutes. Watch it, you will not regret it. If you’d like to know more about epidemiological models and R, I strongly recommend Epidemics: Models and Data in R by Ottar N. Bjornstad which comes with R codes and Shiny apps.

-

Below I follow Sylvain’s talk step by step and the comments are basically transcripts of his words.

+

The talks were recorded, and the videos can be found on Youtube here along with the questions/answers there.

+

I am no expert in genomics, and Mircea, Samuel Alizon and their colleagues have made their work available in several reports that include R codes. Here, I will focus my attention on Sylvain’s work with Sébastien Lion and Troy Day. I am gonna reproduce their analyses in R. In brief, Sylvain Gandon and colleagues use a SIR model to illustrate the effect of different strategies of social distancing on the epidemic. Sylvain’s talk starts at the 49th minute here and lasts approx. 25 minutes. Watch it, you will not regret it. If you’d like to know more about epidemiological models and R, I strongly recommend Epidemics: Models and Data in R by Ottar N. Bjornstad which comes with R codes and Shiny apps. A quick introduction to SIR models in relation to covid-19 can be found there.

+

Below I follow Sylvain’s talk step by step and the comments are basically transcripts of his words. GitHub does not support inline formulas, which explains why they look funny (latex style) below; you might want to have a look the html file or the pdf file or re-run the analyses with the Rmd file in RStudio.

SIR model

-

We start with a classic SIR model omitting birth and death. We have three compartments \(S\), \(I\) and \(R\), for the susceptible, infected and recovered stocks, with \(N = S + I + R\). Total population size is constant through time. The dynamic of the system is described by a system of ordinary differential equations (ODEs): \[\frac{dS}{dt}=-\beta I S/N, \frac{dI}{dt}=\beta I S/N - \gamma I, \frac{dR}{dt}=-\gamma I\] Let us have a look to the dynamic of this system. We load the deSolve package for solving ODEs and the tidyverse package for data manipulation and visualisation, and set defaults for data visualisation:

+

We start with a classic SIR model omitting birth and death. We have three compartments \(S\), \(I\) and \(R\), for the susceptible, infected and recovered stocks, with \(N = S + I + R\). Total population size is constant through time. The dynamic of the system is described by a system of ordinary differential equations (ODEs): \[\frac{dS}{dt}=-\beta I S/N, \frac{dI}{dt}=\beta I S/N - \gamma I, \frac{dR}{dt}=\gamma I\] Let us have a look to the dynamic of this system. We load the deSolve package for solving ODEs and the tidyverse package for data manipulation and visualisation, and set defaults for data visualisation:

library(deSolve)
 library(tidyverse)
 theme_set(theme_light(base_size=16))
@@ -427,7 +427,7 @@

SIR model

y = 'Prevalence') + scale_colour_viridis_d(name = NULL, breaks = c('susceptible','infected')) -

+

The unimodal infected curve shows that the epidemic occurs between 50 and 100 days. Then, because recovery means immunity for life in this model, the number of susceptible individuals at risk of getting infected decreases through time, and the epidemic eventually stops because there is no more susceptible host to infect.

@@ -457,7 +457,7 @@

No social distancing

y = 'Prevalence') + scale_colour_viridis_d(name = NULL, breaks = c('susceptible','infected')) -

+

We may use a prevalence on log scale which shows an exponential growth in the early stage of the epidemic:

init       <- c(S = 1-1e-5, I = 1e-5, R = 0)
 parameters <- c(mu = 0, N = 1, beta = 0.2, gamma = 0.1)
@@ -480,7 +480,7 @@ 
No social distancing

   scale_colour_viridis_d(name = NULL,
                          breaks = c('susceptible','infected')) +
   coord_trans(y = "log10")
-

+

If we increase \(R_0\), the epidemic is faster:

# R0 = 2
 init       <- c(S = 1-1e-5, I = 1e-5, R = 0)
@@ -581,7 +581,7 @@ 
No social distancing

 
 library(patchwork)
 linearRknot2 + logRknot2 + linearRknot35 + logRknot35
-

+

On the right, we see that the slope of the infected increases when we go from \(R_0=2\) (top right) to \(R_0=3.5\) (bottom right). Also on the left, we see that the drop in the susceptible is faster when we go from \(R_0=2\) (top left) to \(R_0=3.5\) (bottom left). Flattening the curve is decreasing the \(R_0\) and going from the bottom left to the top left panel.

@@ -639,7 +639,7 @@

With permanent social distancing

y = 'Prevalence') + scale_colour_viridis_d(name = NULL, breaks = c('susceptible','infected')) -

+

Now we introduce social distancing after 30 days and have increasing social distancing with \(c_{\beta} = 0.3, 0.6, 0.9\). We get:

## c_beta = 0.3
 
@@ -824,8 +824,35 @@ 
With permanent social distancing

 (linearcb03 | logcb03) / 
 (linearcb06 | logcb06) /
 (linearcb09 | logcb09)
-

+

When you do social distancing (the grey shaded area), you flatten the curve of the epidemic (top left to middle left to bottom left), and eventually, you eradicate the disease (bottom right). The thing with permanent social distancing is that it is unfeasible in the long term (socially, economically, etc..).

+

We can dig a bit deeper and see the effect of social distancing on the total epidemic size. Calculating total epidemic size boils down to finding out the equilibrium of the ODEs system. To do so, we will use numerical integration and the rootSolve package:

+
library(rootSolve)
+
+init <- c(S = 1-1e-5, I = 1e-5, R = 0)
+lgrid <- 100
+grid <- seq(0.01, 0.7, length = lgrid)
+res <- rep(NA,lgrid)
+index <- 1
+for (i in grid){
+  parameters <- c(N = 1, beta = 0.25, gamma = 0.1, cbeta = i, cgamma = 0, tSD = 30)
+  times      <- c(0, 300)
+  out <- runsteady(y = init, times = times, func = sirSD, parms = parameters)  
+  res[index] <- out$y[3]
+  index <- index + 1
+}
+

We may visualize the result:

+
res %>% 
+  as_tibble() %>%
+  add_column(cbeta = grid) %>%
+  ggplot(aes(x = cbeta, y = value)) + 
+  geom_line(lwd=1.5) + 
+  labs(title = '',
+       x = 'Intensity of social distancing',
+       y = 'Total epidemic size') + 
+  scale_y_continuous(breaks = seq(0, 1, by = 0.1))
+

+

If there is no social distancing, at the end of the epidemic, almost \(90\%\) individuals get infected. In contrast, if we apply a social distancing of \(c_{\beta} = 0.6\) or higher, then the number of infected individuals is almost null.

What if stop social distancing after some time, or do social distancing repeatedly over several periods of time?

@@ -863,7 +890,7 @@

Transitory social distancing

}

With \(R_0=2.5\) and \(c_{\beta} = 0.8\), we get:

init       <- c(S = 1-1e-5, I = 1e-5, R = 0)
-parameters <- c(mu = 0, N = 1, beta = 0.25, gamma = 0.1, cbeta = 0.8, cgamma = 0, tSDmin = 30, tSDmax = 180)
+parameters <- c(N = 1, beta = 0.25, gamma = 0.1, cbeta = 0.8, cgamma = 0, tSDmin = 30, tSDmax = 180)
 times      <- seq(0, 400, by = 1)
 rects <- data.frame(xstart = 30,
                 xend = 180,
@@ -890,9 +917,10 @@ 
Transitory social distancing

        y = 'Prevalence') + 
   scale_colour_viridis_d(name = NULL,
                          breaks = c('susceptible','infected'))
-

And on the log scale:

+

+

And on the log scale:

init       <- c(S = 1-1e-5, I = 1e-5, R = 0)
-parameters <- c(mu = 0, N = 1, beta = 0.25, gamma = 0.1, cbeta = 0.8, cgamma = 0, tSDmin = 30, tSDmax = 180)
+parameters <- c(N = 1, beta = 0.25, gamma = 0.1, cbeta = 0.8, cgamma = 0, tSDmin = 30, tSDmax = 180)
 times      <- seq(0, 400, by = 1)
 rects <- data.frame(xstart = 30,
                 xend = 180,
@@ -920,10 +948,10 @@ 
Transitory social distancing

   scale_colour_viridis_d(name = NULL,
                          breaks = c('susceptible','infected')) +
   coord_trans(y = "log10")
-

+

We see that social distancing does have an effect by inverting the slope of infected individuals from 30 days onward; however, when we stop social distancing at 180 days then there are still a lot of susceptible individuals at risk of being infected, and there is a second wave occuring. Now if we decrease the intensity of social distancing, say with \(c_{\beta} = 0.4\), then there will not be a second wave:

init       <- c(S = 1-1e-5, I = 1e-5, R = 0)
-parameters <- c(mu = 0, N = 1, beta = 0.25, gamma = 0.1, cbeta = 0.4, cgamma = 0, tSDmin = 30, tSDmax = 180)
+parameters <- c(N = 1, beta = 0.25, gamma = 0.1, cbeta = 0.4, cgamma = 0, tSDmin = 30, tSDmax = 180)
 times      <- seq(0, 400, by = 1)
 rects <- data.frame(xstart = 30,
                 xend = 180,
@@ -951,10 +979,10 @@ 
Transitory social distancing

   scale_colour_viridis_d(name = NULL,
                          breaks = c('susceptible','infected')) +
   coord_trans(y = "log10")
-

+

This is because during the social distancing period (the grey shaded area in the figure below), the intensity is strong enough for a lot of immunity to build up (herd immunity), and the density of susceptible individuals decreases which means less opportunity for transmission:

init       <- c(S = 1-1e-5, I = 1e-5, R = 0)
-parameters <- c(mu = 0, N = 1, beta = 0.25, gamma = 0.1, cbeta = 0.4, cgamma = 0, tSDmin = 30, tSDmax = 180)
+parameters <- c(N = 1, beta = 0.25, gamma = 0.1, cbeta = 0.4, cgamma = 0, tSDmin = 30, tSDmax = 180)
 times      <- seq(0, 400, by = 1)
 rects <- data.frame(xstart = 30,
                 xend = 180,
@@ -981,7 +1009,31 @@ 
Transitory social distancing

        y = 'Prevalence (log)') + 
   scale_colour_viridis_d(name = NULL,
                          breaks = c('susceptible','infected'))
-

+

+

Let us look at the effect of transitory social distancing on the total epidemic size:

+
init <- c(S = 1-1e-5, I = 1e-5, R = 0)
+lgrid <- 100
+grid <- seq(0.01, 0.7, length = lgrid)
+res <- rep(NA,lgrid)
+index <- 1
+for (i in grid){
+  parameters <- c(N = 1, beta = 0.25, gamma = 0.1, cbeta = i, cgamma = 0, tSDmin = 30, tSDmax = 180)
+  times      <- c(0, 350)
+  out <- runsteady(y = init, times = times, func = sirSDtransitory, parms = parameters)  
+  res[index] <- out$y[3]
+  index <- index + 1
+}
+

We may visualize the result:

+
res %>% 
+  as_tibble() %>%
+  add_column(cbeta = grid) %>%
+  ggplot(aes(x = cbeta, y = value)) + 
+  geom_line(lwd=1) + 
+  labs(title = '',
+       x = 'Intensity of social distancing',
+       y = 'Total epidemic size') +
+  ylim(0,1)
+

Even though this optimal social distancing intensity does the job, it is still unsatisfying because at the peak of the epidemic, there is still a large proportion of the population which is infected.

@@ -1054,7 +1106,7 @@

Social distancing in multiple bouts

y = 'Prevalence') + scale_colour_viridis_d(name = NULL, breaks = c('susceptible','infected')) -

+

And on the log scale:

init <- c(S = 1-1e-5, I = 1e-5, R = 0)
 parameters <- c(mu = 0, N = 1, beta = 0.25, gamma = 0.1, cbeta = 0.9, cgamma = 0, 
@@ -1084,11 +1136,11 @@ 
Social distancing in multiple bouts

   labs(title = 'Multiple bouts social distancing',
        subtitle = expression(paste('Log scale',~c[beta],'=0.9')),
        x = 'Time (days)',
-       y = 'Prevalence (log)') + 
+       y = 'Prevalence') + 
   scale_colour_viridis_d(name = NULL,
                          breaks = c('susceptible','infected')) +
   coord_trans(y = "log10")
-

+

We see that during the periods of social distancing (grey shaded areas), we decrease the size of the epidemic, then in between periods of social distancing, the epidemic might come back, or you might succeed in eradicating the epidemic. Also, the maximum level of prevalence is lower than in permanent or transitory social distancing, which might help in decreasing the number of patients on hospitals and buying some time to develop a vaccine.