# Lancement du groupe en épidémiologie et en santé publique du centre de recherches mathématiques

Ce vendredi, lancement des activités du groupe en épidémiologie et en santé publique du centre de recherches mathématiques de Montréal. D’autres évènements seront à venir…

# Talk at the Paris Machine Learning Meetup

This Wednesday, I will be presenting at the Paris Machine Learning Meetup our recent work on optimal control and COVID-19.

The slides can be downloaded from here.

# Automne 2020, à distance

C’est maintenant officiel, la session d’automne se fera (elle aussi) à distance, majoritairement. On en saura plus cet été, je pense…

# COVID19 pandemic control: balancing detection policy and lockdown intervention under ICU sustainability

After almost two months, with Romuald Elie, Chi (Tran Viet Chi) and Mathieu Laurière, we finally have a draft of the paper related to our recent work, entitled

Our model is an extention of the classical SIR model, but more realistic for the COVID-19 (pronounded “cider”)

We use scenarios to see the impact on various quantities

but also optimal control to see the best strategy, when it comes to lockdown, and testing. All comments are welcome…

# Testing for Covid-19 in the U.S.

For almost a month, on a daily basis, we are working with colleagues (Romuald, Chi and Mathieu) on modeling the dynamics of the recent pandemic. I learn of lot of things discussing with them, but we keep struggling with the tests. Paul, in Montréal, helped me a little bit, but I think we will still have to more to get a better understand. To but honest, we stuggle with two very simple questions

• how many people are tested on a daily basis ?

Recently, I discovered Modelling COVID-19 exit strategies for policy makers in the United Kingdom, which is very close to what we try to do… and in the document two interesting scenarios are discussed, with, for the first one, “1 million ‘reliable’ daily tests are deployed” (in the U.K.) and “5 million ‘useless’ daily tests are deployed”. There are about 65 millions unhabitants in the U.K. so we talk here about 1.5% people tested, on a daily basis, or 7.69% people ! It could make sense, but our question was, at some point, is that realistic ? where are we today with testing ? In the U.S. https://covidtracking.com/ collects interesting data, on a daily basis, per state.

url = "https://raw.githubusercontent.com/COVID19Tracking/covid-tracking-data/master/data/states_daily_4pm_et.csv" download.file(url,destfile="covid.csv") base = read.csv("covid.csv")

Unfortunately, there is no information about the population. That we can find on wikipedia. But in that table, the state is given by its full name (and the symbol in the previous dataset). So we new also to match the two datasets properly,

url="https://en.wikipedia.org/wiki/List_of_states_and_territories_of_the_United_States_by_population" download.file(url,destfile = "popUS.html") #pas contaminé 2/3 R=3 library(XML) tables=readHTMLTable("popUS.html") T=tables[[1]][3:54,c("V3","V4")] names(T)=c("state","pop") url="https://en.wikipedia.org/wiki/List_of_U.S._state_abbreviations" download.file(url,destfile = "nameUS.html") tables=readHTMLTable("nameUS.html") T2=tables[[1]][13:63,c(1,4)] names(T2)=c("state","symbol") T=merge(T,T2) T$population = as.numeric(gsub(",", "", T$pop, fixed = TRUE)) names(base)[2]="symbol" base = merge(base,T[,c("symbol","population")])

Now our dataset is fine… and we can get a function to plot the number of people tested in the U.S. (cumulated). Here, we distinguish between the positive and the negative,

drawing = function(st ="NY"){ sbase=base[base$symbol==st,c("date","positive","negative","population")] sbase$DATE = as.Date(as.character(sbase$date),"%Y%m%d") sbase=sbase[order(sbase$DATE),] par(mfrow=c(1,2)) plot(sbase$DATE,(sbase$positive+sbase$negative)/sbase$population,ylab="Proportion Test (/population of state)",type="l",xlab="",col="blue",lwd=3) lines(sbase$DATE,sbase$positive/sbase$population,col="red",lwd=2) legend("topleft",c("negative","positive"),lwd=2,col=c("blue","red"),bty="n") title(st) plot(sbase$DATE,sbase$positive/(sbase$positive+sbasenegative),ylab="Ratio of positive tests",ylim=c(0,1),type="l",xlab="",col="black",lwd=3) title(st)} Let us start with New York drawing("NY") As at now, 4% of the entiere population got tested… over 6 weeks…. The graph on the right is the proportion of people who tested positive… I won’t get back on that one here today, I keep it for our work. In New Jersey, we got about 2.5% of the entiere population tested, overall, drawing("NJ") Let us try a last one, Florida drawing("FL") As at today, it is 1.5% of the population, over 6 weeks. Overall, in the U.S. less than 0.1% people are tested, on a daily basis. Which is far from the 1.5% in the U.K. scenarios. Now, here come the second question, • what are we actually testing for ? On that one, my experience in biology is… very limited, and Paul helped me. He mentioned this morning a nice report, from a lab in UC Berkeley One of my question was for instance, if you get tested positive, and you do it again, can you test negative ? Or, in the context of our data, do we test different people ? are some people tested on a regular basis (perhaps every week) ? For instance, with antigen tests (Reverse Transcription Quantitative Polymerase Chain Reaction (RT-qPCR) – also called molecular or PCR – Polymerase Chain Reaction – test) we test if someone is infectious, while with antibody test (using serological immunoassays that detect viral-specific antibodies — Immunoglobin M (IgM) and G (IgG) — also called serology test), we test for immunity. Which is rather different… I have no idea what we have in our database, to be honest… and for the past six weeks, I have seen a lot of databases, and most of the time, I don’t know how to interpret, I don’t know what is measured… and it is scary. So, so far, we try to do some maths, to test dynamics by tuning parameters “the best we can” (and not estimate them). But if anyone has good references on testing, in the context of Covid-19 (for instance on specificity, sensitivity of all those tests) I would love to hear about it ! # “Family History” and Life Insurance Today and tomorrow, I will attend the Online International Conference in Actuarial science, data science and finance, organised by colleagues in Lyon. But I won’t be in Lyon, I will be at home, in Montréal… I will give a talk on wednesday afternoon, on a joint paper with Ewen Gallic and Olivier Cabrignac. Slides are available here, and if I can get a copy of the video, I will share it… # Telematics: a revolution ? Tomorrow morning, Laurence Barry will present some joint work at the Online International Conference in Actuarial science, data science and finance, organised by colleagues in Lyon. The paper, “Personalization as a Promise: Can Big Data Change the Practice of Insurance?” is online, as a Working Paper of the PARI chair, The purpose of this paper is to measure the impact of technologies from the Big Bang. data on thehe pricing ofhe products car insurance. The first part describes how the aggregated view buildsuit by statistics enables highlighting invisible regularities at the individual level. Despite a very granular segmentations in automobile insurance, the approach remained classificatory, hypothesizing the risk identity of individuals from the same class. The second part highlights the reversal of big data-induced perspective in the’analysis ofgiven ; awith theur volume and the new algorithms, the aggregate viewpoint is questioned.. The hypothesis of class homogeneity is becoming increasingly difficult to test. maintain, especially since predictive analysis boasts the ability to predict the rs results at the individual level. The third part is studying the’influence of telematics boxes able to import the new pinsurance aradigm automobile. However, a reading of the most recent research articles on a pricing automobile including this new monter that the epistemological leap, at least for now, has not taken place. # INF7100, à distance Suite au contexte de la covid-19, la session d’été commencera une semaine plus tard. Le cours INF7100 Initiation à la science des données et à l’intelligence artificielle – qui devait commencer mardi 28 avril commencera le 5 mai prochain. Les étudiants inscrits recevront un lien pour se connecter, pour une présentation globale. On le donnera conjointement avec Marie-Jean Meurs (du département d’informatique) et Jean-Hugues Roy (de l’école des médias). Je posterais des informations dans les semaines à venir sur les points que je vais aborder (essentiellement sur la science des données). # De la démarche scientifique en période de crise Dans une conférence donnée le 13 février 2020[i], intitulée contre la méthode, Didier Raoult affirmait « moi je n’ai jamais fait d’essais randomisés […] faire ça sur des maladies infectieuses, ça n’a pas de sens ». Cette vision était reprise dans une tribune plus détaillée, où face à « la méthode » (et « aux mathématiques »), Didier Raoult défendait (ce qu’il appelait) « la morale [et] l’humanisme » du serment d’Hippocrate. Comme il le rappelle, faire des groupes de contrôle, c’est « dire au malade qu’on va lui donner au hasard soit le médicament dont on sait qu’il marche, soit le médicament dont on ne sait pas s’il marche » (Raoult (2020a, 2020b)). Alors que cette méthode d’expériences randomisées est aujourd’hui saluée dans toutes les disciplines – comme le rappelle le prix Nobel d’économie attribué en 2019 à Esther Duflo, Michael Kremer et Abhijit Banerjee – comment un chercheur peut-il prendre une telle position, aujourd’hui ? Continue reading De la démarche scientifique en période de crise # Modeling Pandemics (3) In Statistical Inference in a Stochastic Epidemic SEIR Model with Control Intervention, a more complex model than the one we’ve seen yesterday was considered (and is called the SEIR model). Consider a population of size $N$, and assume that $S$ is the number of susceptible, $E$ the number of exposed, $I$ the number of infectious, and $R$ for the number recovered (or immune) individuals, \displaystyle{\begin{aligned}{\frac {dS}{dt}}&=-\beta {\frac {I}{N}}S\\[8pt]{\frac {dE}{dt}}&=\beta {\frac {I}{N}}S-aE\\[8pt]{\frac {dI}{dt}}&=aE-b I\\[8pt]{\frac {dR}{dt}}&=b I\end{aligned}}Between $S$ and $I$, the transition rate is $\beta I$, where $\beta$ is the average number of contacts per person per time, multiplied by the probability of disease transmission in a contact between a susceptible and an infectious subject. Between $I$ and $R$, the transition rate is $b$ (simply the rate of recovered or dead, that is, number of recovered or dead during a period of time divided by the total number of infected on that same period of time). And finally, the incubation period is a random variable with exponential distribution with parameter $a$, so that the average incubation period is $a^{-1}$. Probably more interesting, Understanding the dynamics of ebola epidemics suggested a more complex model, with susceptible people $S$, exposed $E$, Infectious, but either in community $I$, or in hospitals $H$, some people who died $F$ and finally those who either recover or are buried and therefore are no longer susceptible $R$. Thus, the following dynamic model is considered\displaystyle{\begin{aligned}{\frac {dS}{dt}}&=-(\beta_II+\beta_HH+\beta_FF)\frac{S}{N}\\[8pt]\frac {dE}{dt}&=(\beta_II+\beta_HH+\beta_FF)\frac{S}{N}-\alpha E\\[8pt]\frac {dI}{dt}&=\alpha E+\theta\gamma_H I-(1-\theta)(1-\delta)\gamma_RI-(1-\theta)\delta\gamma_FI\\[8pt]\frac {dH}{dt}&=\theta\gamma_HI-\delta\lambda_FH-(1-\delta)\lambda_RH\\[8pt]\frac {dF}{dt}&=(1-\theta)(1-\delta)\gamma_RI+\delta\lambda_FH-\nu F\\[8pt]\frac {dR}{dt}&=(1-\theta)(1-\delta)\gamma_RI+(1-\delta)\lambda_FH+\nu F\end{aligned}}In that model, parameters are $\alpha^{-1}$ is the (average) incubation period (7 days), $\gamma_H^{-1}$ the onset to hospitalization (5 days), $\gamma_F^{-1}$ the onset to death (9 days), $\gamma_R^{-1}$ the onset to “recovery” (10 days), $\lambda_F^{-1}$ the hospitalisation to death (4 days) while $\lambda_R^{-1}$ is the hospitalisation to recovery (5 days), $\eta^{-1}$ is the death to burial (2 days). Here, numbers are from Understanding the dynamics of ebola epidemics (in the context of ebola). The other parameters are $\beta_I$ the transmission rate in community (0.588), $\beta_H$ the transmission rate in hospital (0.794) and $\beta_F$ the transmission rate at funeral (7.653). Thus epsilon = 0.001 Z = c(S = 1-epsilon, E = epsilon, I=0,H=0,F=0,R=0) p=c(alpha=1/7*7, theta=0.81, delta=0.81, betai=0.588, betah=0.794, blambdaf=7.653,N=1, gammah=1/5*7, gammaf=1/9.6*7, gammar=1/10*7, lambdaf=1/4.6*7, lambdar=1/5*7, nu=1/2*7) If $\boldsymbol{Z}=(S,E,I,H,F,R)$, if we write $$\frac{\partial \boldsymbol{Z}}{\partial t} = SEIHFR(\boldsymbol{Z})$$where $SEIHFR$ is SEIHFR = function(t,Z,p){ S=Z[1]; E=Z[2]; I=Z[3]; H=Z[4]; F=Z[5]; R=Z[6] alpha=p["alpha"]; theta=p["theta"]; delta=p["delta"] betai=p["betai"]; betah=p["betah"]; gammah=p["gammah"] gammaf=p["gammaf"]; gammar=p["gammar"]; lambdaf=p["lambdaf"] lambdar=p["lambdar"]; nu=p["nu"]; blambdaf=p["blambdaf"] N=S+E+I+H+F+R dS=-(betai*I+betah*H+blambdaf*F)*S/N dE=(betai*I+betah*H+blambdaf*F)*S/N-alpha*E dI=alpha*E-theta*gammah*I-(1-theta)*(1-delta)*gammar*I-(1-theta)*delta*gammaf*I dH=theta*gammah*I-delta*lambdaf*H-(1-delta)*lambdaf*H dF=(1-theta)*(1-delta)*gammar*I+delta*lambdaf*H-nu*F dR=(1-theta)*(1-delta)*gammar*I+(1-delta)*lambdar*H+nu*F dZ=c(dS,dE,dI,dH,dF,dR) list(dZ)} We can solve it, or at least study the dynamics from some starting values library(deSolve) times = seq(0, 50, by = .1) resol = ode(y=Z, times=times, func=SEIHFR, parms=p) For instance, the proportion of people infected is the following plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="",col="red") lines(resol[,"time"],resol[,"H"],col="blue") # Modeling pandemics (2) When introducing the SIR model, in our initial post, we got an ordinary differential equation, but we did not really discuss stability, and periodicity. It has to do with the Jacobian matrix of the system. But first of all, we had three equations for three function, but actually$$\displaystyle{{\frac{dS}{dt}}+{\frac {dI}{dt}}+{\frac {dR}{dt}}=0}$$so it means that our problem is here simply in dimension 2. Hence\displaystyle {\begin{aligned}&X={\frac {dS}{dt}}=\mu(N-S)-{\frac {\beta IS}{N}},\\[6pt]&Y={\frac {dI}{dt}}={\frac {\beta IS}{N}}-(\mu+\gamma)I\end{aligned}}and therefore, the Jacobian of the system is$$\begin{pmatrix}\displaystyle{\frac{\partial X}{\partial S}}&\displaystyle{\frac{\partial X}{\partial I}}\\[9pt]\displaystyle{\frac{\partial Y}{\partial S}}&\displaystyle{\frac{\partial Y}{\partial I}}\end{pmatrix}=\begin{pmatrix}\displaystyle{-\mu-\beta\frac{I}{N}}&\displaystyle{-\beta\frac{S}{N}}\\[9pt]\displaystyle{\beta\frac{I}{N}}&\displaystyle{\beta\frac{S}{N}-(\mu+\gamma)}\end{pmatrix}$$We should evaluate the Jacobian at the equilibrium, i.e. $$S^\star=\frac{\gamma+\mu}{\beta}=\frac{1}{R_0}$$and$$I^\star=\frac{\mu(R_0-1)}{\beta}$$We should then look at eigenvalues of the matrix. Our very last example was times = seq(0, 100, by=.1) p = c(mu = 1/100, N = 1, beta = 50, gamma = 10) start_SIR = c(S=0.19, I=0.01, R = 0.8) resol = ode(y=start_SIR, t=times, func=SIR, p=p) plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="") We can compute values at the equilibrium mu=p["mu"]; beta=p["beta"]; gamma=p["gamma"] N=1 S = (gamma + mu)/beta I = mu * (beta/(gamma + mu) - 1)/beta and the Jacobian matrix J=matrix(c(-(mu + beta * I/N),-(beta * S/N), beta * I/N,beta * S/N - (mu + gamma)),2,2,byrow = TRUE) Now, if we look at the eigenvalues, eigen(J)values [1] -0.024975+0.6318831i -0.024975-0.6318831i

or more precisely $2\pi/b$ where $a\pm ib$ are the conjuguate eigenvalues

2 * pi/(Im(eigen(J)$values[1])) [1] 9.943588 we have a damping period of 10 time lengths (10 days, or 10 weeks), which is more or less what we’ve seen above, The graph above was obtained using p = c(mu = 1/100, N = 1, beta = 50, gamma = 10) start_SIR = c(S=0.19, I=0.01, R = 0.8) resol = ode(y=start_SIR, t=times, func=SIR, p=p) plot(resol[1:1e5,"time"],resol[1:1e5,"I"],type="l",xlab="time",ylab="",lwd=3,col="red") yi=resol[,"I"] dyi=diff(yi) i=which((dyi[2:length(dyi)]*dyi[1:(length(dyi)-1)])&lt;0) t=resol[i,"time"] arrows(t[2],.008,t[4],.008,length=.1,code=3) If we look carefully. at the begining, the duration is (much) longer than 10 (about 13)… but it does converge towards 9.94 plot(diff(t[seq(2,40,by=2)]),type="b") abline(h=2 * pi/(Im(eigen(J)$values[1]))

So here, theoretically, every 10 weeks (assuming that our time length is a week), we should observe an outbreak, smaller than the previous one. In practice, initially it is every 13 or 12 weeks, but the time to wait between outbreaks decreases (until it reaches 10 weeks).

# Modeling pandemics (1)

The most popular model to model epidemics is the so-called SIR model – or Kermack-McKendrick. Consider a population of size $N$, and assume that $S$ is the number of susceptible, $I$ the number of infectious, and $R$ for the number recovered (or immune) individuals, \displaystyle {\begin{aligned}&{\frac {dS}{dt}}=-{\frac {\beta IS}{N}},\\[6pt]&{\frac {dI}{dt}}={\frac {\beta IS}{N}}-\gamma I,\\[6pt]&{\frac {dR}{dt}}=\gamma I,\end{aligned}}so that $$\displaystyle{{\frac{dS}{dt}}+{\frac {dI}{dt}}+{\frac {dR}{dt}}=0}$$which implies that $S+I+R=N$. In order to be more realistic, consider some (constant) birth rate $\mu$, so that the model becomes\displaystyle {\begin{aligned}&{\frac {dS}{dt}}=\mu(N-S)-{\frac {\beta IS}{N}},\\[6pt]&{\frac {dI}{dt}}={\frac {\beta IS}{N}}-(\gamma+\mu) I,\\[6pt]&{\frac {dR}{dt}}=\gamma I-\mu R,\end{aligned}}Note, in this model, that people get sick (infected) but they do not die, they recover. So here, we can model chickenpox, for instance, not SARS.

The dynamics of the infectious class depends on the following ratio:$$\displaystyle{R_{0}={\frac {\beta }{\gamma +\mu}}}$$ which is the so-called basic reproduction number (or reproductive ratio). The effective reproductive ratio is $R_0S/N$, and the turnover of the epidemic happens exactly when $R_0S/N=1$, or when the fraction of remaining susceptibles is $R_0^{-1}$. As shown in Directly transmitted infectious diseases:Control by vaccination, if $S/N the disease (the number of people infected) will start to decrease.

Want to see it  ? Start with

mu = 0 beta = 2 gamma = 1/2

for the parameters. Here,  $R_0=4$. We also need starting values

epsilon = .001 N = 1 S = 1-epsilon I = epsilon R = 0

Then use the ordinary differential equation solver, in R. The idea is to say that $\boldsymbol{Z}=(S,I,R)$ and we have the gradient $$\frac{\partial \boldsymbol{Z}}{\partial t} = SIR(\boldsymbol{Z})$$where $SIR$ is function of the various parameters. Hence, set

p = c(mu = 0, N = 1, beta = 2, gamma = 1/2) start_SIR = c(S = 1-epsilon, I = epsilon, R = 0)

The we must define the time, and the function that returns the gradient,

times = seq(0, 10, by = .1) SIR = function(t,Z,p){ S=Z[1]; I=Z[2]; R=Z[3]; N=S+I+R mu=p["mu"]; beta=p["beta"]; gamma=p["gamma"] dS=mu*(N-S)-beta*S*I/N dI=beta*S*I/N-(mu+gamma)*I dR=gamma*I-mu*R dZ=c(dS,dI,dR) return(list(dZ))}

To solve this problem use

library(deSolve) resol = ode(y=start_SIR, times=times, func=SIR, parms=p)

We can visualize the dynamics below

par(mfrow=c(1,2)) t=resol[,"time"] plot(t,resol[,"S"],type="l",xlab="time",ylab="") lines(t,resol[,"I"],col="red") lines(t,resol[,"R"],col="blue") plot(t,t*0+1,type="l",xlab="time",ylab="",ylim=0:1) polygon(c(t,rev(t)),c(resol[,"R"],rep(0,nrow(resol))),col="blue") polygon(c(t,rev(t)),c(resol[,"R"]+resol[,"I"],rev(resol[,"R"])),col="red")

We can actually also visualize the effective reproductive number is $R_0S/N$, where

R0=p["beta"]/(p["gamma"]+p["mu"])

The effective reproductive number is on the left, and as we mentioned above, when we reach 1, we actually reach the maximum of the infected,

plot(t,resol[,"S"]*R0,type="l",xlab="time",ylab="") abline(h=1,lty=2,col="red") abline(v=max(t[resol[,"S"]*R0&gt;=1]),col="darkgreen") points(max(t[resol[,"S"]*R0&gt;=1]),1,pch=19) plot(t,resol[,"S"],type="l",xlab="time",ylab="",col="grey") lines(t,resol[,"I"],col="red",lwd=3) lines(t,resol[,"R"],col="light blue") abline(v=max(t[resol[,"S"]*R0&gt;=1]),col="darkgreen") points(max(t[resol[,"S"]*R0&gt;=1]),max(resol[,"I"]),pch=19)

And when adding a $\mu$ parameter, we can obtain some interesting dynamics on the number of infected,

times = seq(0, 100, by=.1) p = c(mu = 1/100, N = 1, beta = 50, gamma = 10) start_SIR = c(S=0.19, I=0.01, R = 0.8) resol = ode(y=start_SIR, t=times, func=SIR, p=p) plot(resol[,"time"],resol[,"I"],type="l",xlab="time",ylab="")