45:["$","title",null,{"children":"Interactive SEIR Model of the US | S. Scott Collis | Observable"}] 4d:T66a,md`Solves a stochastic ${tex`\textsf{SEIR}`} compartment model for disease progression based on the formulation of [Yeo, Yeo & Yeo]( https://doi.org/10.1101/2020.03.24.20043026). A nice feature of this formulation is that it tracks the number of infected individuals in the community ${tex`I_c`} separately from those who have been identified through testing ${tex`I_i`} and are therfore *issolated* and no longer contributing to new exposures. This allows one to see the advantage of testing to reduce the number of infectious individuals who can expose the susceptable population. The Wikipedia entry on [compartmental models](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology) also provides helpful background. This is my first use of JavaScript, Observable notebooks, and Vega-Lite plotting so I'm sure that there are improvements to be made. For those that are old-school, I have also implemented a version of this model in FORTRAN and it is available at https://github.com/sscollis/sanseir, although it lacks some of the features of the JavaScript version available here. In particular, the JavaScript/Observable version presented here differs by * Using an updated mortality rate for the US, * Using incubation and infection periods and statistics from a Johns Hopkins University study, and * Tuning of ${tex`R_o`} and the testing fractions to notionally coorespond to national changes while also matching to US national data. * A recent update to allow for time varying mortality to model the changing demographics of the infected population. The set of ordinary differential equations is given by`4e:T541,tex`\frac{dS}{dt} = -\tilde\beta I_c \frac{S}{N} \\[1ex] \frac{dE}{dt} = \tilde\beta I_c \frac{S}{N} - \alpha E \\[1ex] \frac{dI_i}{dt} = \tilde\alpha\rho E - \tilde\gamma I_i \\[1ex] \frac{dI_c}{dt} = \tilde\alpha(1-\rho) E - \tilde\gamma I_c \\[1ex] \frac{dR}{dt} = \tilde\gamma(1-c) (I_i + I_c) \\[1ex] \frac{dD}{dt} = \tilde\gamma c (I_i + I_c)` // // There is no need to track deaths separately in h versus c right now unless // the mortality is different in those populations (which it is). // To do so in a meaningful way would likely require that we track // I_h = in hospital (the mortality is thought to be 1.4%) // I_i = issolated in the community (would not contribute to new exposures) // I_c = in the community // Then // D_h = 0.014 * I_h // D_c = 0.0066 * (I_h+I_i+I_c) - 0.014*I_h // Therefore // c_c = D_c/(I_i+I_c) = (0.0066 * (I_h+I_i+I_c) - 0.014*I_h)/(I_i+I_c) // // A key here is that increasing testing should increase I_i but not really // I_h as the number of people in a hospital isn't really dependent on testing // but instead on the rate of severe disease. // // \frac{dR_h}{dt} = \tilde\gamma(1-c) I_h \\[1ex] // \frac{dR_c}{dt} = \tilde\gamma(1-c) I_c \\[1ex] // \frac{dD_h}{dt} = \tilde\gamma c I_h \\[1ex] // \frac{dD_c}{dt} = \tilde\gamma c I_c `4f:Tef0,md`Here, ${tex`S`} is the number of susceptable, ${tex`E`}, exposed, ${tex`I`}, infectious, ${tex`R`} recovered, and ${tex`D`} deceased individuals. The key parameters are: Parameter | Description ----------------------|----------------------------------------- ${tex`\tilde\beta`} | contact rate ${tex`\tilde\beta = R_o \tilde\gamma`} ${tex`\tilde\alpha`} | incubation rate ${tex`\rho`} | fraction of exposed that are tested ${tex`\tilde\gamma`} | infectious rate ${tex`c`} | mortality ${tex`I_o`} | Initial number of infectious individuals in the community where ${tex`R_o`} is the reproduction rate. We use the notation ${tex`\tilde\alpha = \alpha^{-1}`} to denote the rate versus the period. The incubation and infectious periods are optionally sampled from a probability distribution (with the minimum allowable value to be 1 day). The following default values are calibrated to the current outbreak data from the United States. The start of the simulation defaults to three days before the first reported case in the US, 1/21/2020, and the simulation runs for 8 months (240 days) by default. I've tuned the baseline parameters to match the US data with parameters that are reasonably within ranges found in the literature. Information regarding mortality is available in a [recent article](// https://www.usnews.com/news/health-news/articles/2020-03-30/odds-of-hospitalization-death-with-covid-19-rise-steadily-with-age-study). The mean infectious and incubation times in the baseline directly use the distributions and values from the recent [Johns Hopkins University study by Lauer et al.](https://github.com/HopkinsIDD/ncov_incubation). ### Notes: * This model does not account for new infected individuals traveling into the country * There is no direct modeling of the influence of age (everything is averaged over age) * 7/28/2020: However, I have added a time varying mortality adjustment which can, indirectly, model the effect of changes in average age (and therefore mortality) of the infected population. * There is no spatial dependancy. In other words, things are averaged over space and there are no net inflows or outflows. This is a very significant approximation! Sophisticated epidemiological models would capture the different mitigations, behaviors, and impacts at the state or county level. * It is assumed that individuals in ${tex`I_i`} are perfectly issolated (typically not true in practice) * The results are extremely sensitive to ${tex`R_o`} and yet this is the primary parameter to model things like social distancing, travel restrictions (mixing), population density, etc. * The different phases of the outbreak are modeled by changes in ${tex`R_o`} and testing, ${tex`\rho`}, and as of the 7/28/2020 update, the mortality adjustment ${tex`M_a`} as discrete functions of time. * The default model parameters do a decent job of matching the observed US date (at least through 4/23/2020). However, this does not imply that the model will successfully predict the future! * I have built in a time offset as recoveries and deaths lag infections and this is required to match the data and a delay of 7 days seems to match the data adequately. * It is also likely that mortality is not only a function of age but also time, as care and treatments improve. Adding this could allow a better fit to the data (this was done on 7/28/2020) ### Recommended Usage: While you are welcome to play with all the parameters in this interactive notebook, initialy I suggest simply trying out the different ***Phases of the Outbreak Scenarios*** which depict the influence of changing social-distancing restrictions as well as the possible impact of increased testing to identify and issolate individuals who are infectious.`50:T6dd,md`Using this radio dial, you can select between different scenarios. Note that these are just for illustration and are not intended to predict specific outcomes. * "Baseline" matches the initial disease outbreak with continued social-distancing * "5/9/2020" a modification of the Baseline that improves match to data with slight increase in ${tex`R_o`} * "6/18/2020" recalibration to add in change in mitigations in mid May * **"7/28/2020"** recalibration to add in change due to Memorial day events and new mitigataions in early July and mortality rate adjustment to reflect reduced mortality due to change in demographics of infected to younger population * "Two-waves" (hypothetical) releases social-distancing on May first and then reinstates on May 15 * "Three-waves" (hypothetical) brings another wave starting on June 1 releasing June 15, and * "Three-waves+testing" (hypothetical) is the same as the prior, but ramping up testing to issolate sick individuals. #### Remarks on 7/28/2020 Calibration Note that in my initial attempt at a **7/28/2020** recalibration, the deaths were systemically overpredicted which was likely due to the heavy dependance of mortality on the age of infected individuals. Events around Memorial day 2020 appeared to disproportionately infect younger individuals and the model did not have an age dependent mortality rate. So, I introduced a mortality adjustment that scales the default mortality as a function of time to model this change in infected demographics and possibly also changes due to new treatments. As time progressses, I find that the mortality rate after the decrease around Memorial Day is rising again, perhaps due to younger individuals now infecting more vulnerable populations.`51:T951,areaPlot = { //const xAxis = vl.x().fieldQ("Time").scale({domain:[0,xscale]}); //const xAxis = vl.x().fieldT("date") //.scale({domain:[toDate(0),toDate(xscale)]});//.axis({format:"%B",title:"Date"}); const data = vl.data({values:Urange,name:"our_data"}); const usdata = vl.markLine({clip:true,tooltip:true,filled:false,interpolate:"natural", point:{opacity:1,shape:"square",filled:false,fill:'white'}}) .data(usData).transform([{calculate:"(datum.Min+datum.Max)/2","as":"Avg"}]); const xScale = { domain:[ toDateTimeObject(toDay(startTime)), toDateTimeObject(xscale) ] }; const xAxis = vl.x().fieldT("date").axis({tickCount:"week",title:"Date"}).scale(xScale); const ymin = 1; const ymax = 450000; //const yScale = { type: "log", base: 10, domain : [ymin, ymax] } ; const yScale = { domain : [ymin, ymax] } ; // {domain:{data:"our_data",field:"Max"}}; // {domain:[0,P+1]}; const color = vl.color().fieldN("Compartment").sort(vl.max("Min").order("descending")); const min = vl.markLine({clip:true,tooltip:true}) .data(Urange) .transform({filter:{and:[{field:"Compartment",oneOf:toPlot1.Compartments}]}}) .encode( color, xAxis, vl.y().fieldQ("Min").scale(yScale), ); const max = vl.markLine({clip:true,tooltip:true}) .data(Urange) .transform({filter:{and:[{field:"Compartment",oneOf:toPlot1.Compartments}]}}) .encode( color, xAxis, vl.y().fieldQ("Max").scale(yScale), ); const area = vl.markArea({clip:true,tooltip:true}) .data(Urange) .transform({filter:{and:[{field:"Compartment",oneOf:toPlot1.Compartments}]}}) .encode( color, xAxis, vl.y().fieldQ("Max").scale(yScale), vl.y2().fieldQ("Min"), vl.opacity().value(0.25) ); const UScases = usdata.encode( vl.x().fieldT("date"), vl.y().fieldQ("dcases"), vl.color().value("green") ); const USdeaths = usdata.encode( vl.x().fieldT("date"), vl.y().fieldQ("ddeaths"), vl.color().value("red") ); var layers = [min,max,area] if (toPlot1.Data.find( (name) => name === "cases" )) layers.push(UScases); if (toPlot1.Data.find( (name) => name === "deaths" )) layers.push(USdeaths); const plot = vl .layer(layers) .width(width-160) .title("Compartments from SEIR Model") .config({axisX:{ticks:true}}) return plot.render(); }52:Ta51,CrAndDeaths = { const ymin = 1; const ymax = 50e6; let yscale = {type:"log",base:10,domain:[ymin,ymax]}; if (ylinear) yscale = {}; //const xAxis = vl.x().fieldQ("Time").scale({domain:[2,xscale]});//.axis({format:"%B",title:"Date"}); const xScale = { domain:[ toDateTimeObject(toDay(startTime)), toDateTimeObject(xscale) ] }; const xAxis = vl.x().fieldT("date").axis({tickCount:"week",title:"Date"}).scale(xScale); const yTitle = {title:"Individuals"} const color = vl.color().fieldN("Compartment").sort(["Cr","D","C","R","Rh","Rc","Dh","Dc"]); //.transform([{filter:{field:"Compartment",oneOf:["D","Cr"]}}, // {calculate:"(datum.Min+datum.Max)/2","as":"Avg"}]) const trans = [{calculate:"(datum.Min+datum.Max)/2","as":"Avg"}, {filter:{and:[{field:"Compartment",oneOf:toPlot.Compartments}, {field:"Min",gt:0.1},{field:"Max",gt:0.1} ]}}]; const min = vl.markLine({clip:true,tooltip:true}) .data(Urange).transform(trans) .encode( vl.y().fieldQ("Min").scale(yscale), xAxis,color ); const max = vl.markLine({clip:true,tooltip:true}) .data(Urange).transform(trans) .encode( vl.y().fieldQ("Max").scale(yscale).axis(yTitle), xAxis,color ); const avg = vl.markLine({clip:true,tooltip:true}) .data(Urange).transform(trans) .encode( vl.y().fieldQ("Avg").scale(yscale), vl.opacity().value(0.5),xAxis,color ); const area = vl.markArea({clip:true,tooltip:true}) .data(Urange).transform(trans) .encode( vl.y().fieldQ("Max").scale(yscale), vl.y2().fieldQ("Min"), vl.opacity().value(0.15),xAxis,color ); const phase = vl.markLine({clip:true,tooltip:true}).encode( vl.y().fieldQ("Phase").axis({tickCount:5}).scale(yscale), ); const UScases = vl.markCircle({clip:true,tooltip:true,filled:false}).data(usData) .encode( vl.x().fieldT("date"), vl.y().fieldQ("cases").scale(yscale), vl.color().value("green") ); const USdeaths = vl.markSquare({clip:true,tooltip:true,filled:false}).data(usData) .transform({filter:{field:"deaths",gt:0}}) .encode( vl.x().fieldT("date"), vl.y().fieldQ("deaths").scale(yscale), vl.color().value("red") ); // // Select which US data to show // var layers=[min,max,area,avg]; if (toPlot.Data.find( (name) => name === "deaths" )) layers.push(USdeaths); if (toPlot.Data.find( (name) => name === "cases" )) layers.push(UScases); return vl .layer(layers) .width(width-175) .title("Cummulative Compartments") .config({axisX:{ticks:true}}) .render() //.toJSON(); }53:T681,tk = { if (scenario=="Baseline" || scenario=="5/9/2020") { var phases = [0, toDay("2020-02-15"), toDay("2020-03-07"), toDay("2020-03-20"), toDay("2020-03-28"), toDay("2020-04-12")]; } else if (scenario=="6/18/2020") { var phases = [0, toDay("2020-02-15"), toDay("2020-03-07"), toDay("2020-03-20"), toDay("2020-03-28"), toDay("2020-04-12"), toDay("2020-05-16")]; } else if (scenario=="7/28/2020") { var phases = [0, toDay("2020-02-15"), toDay("2020-03-07"), toDay("2020-03-20"), toDay("2020-03-28"), toDay("2020-04-12"), toDay("2020-05-16"), toDay("2020-05-31"), toDay("2020-07-13")]; } else if (scenario=="2/28/2021") { var phases = [0, toDay("2020-02-15"), toDay("2020-03-07"), toDay("2020-03-20"), toDay("2020-03-28"), toDay("2020-04-12"), toDay("2020-05-16"), toDay("2020-05-31"), toDay("2020-07-13"), toDay("2020-09-20"), toDay("2021-01-01")]; // // The remainder are hypothetical // } else if (scenario=="Two-waves") { var phases = [0, toDay("2020-02-15"), toDay("2020-03-07"), toDay("2020-03-20"), toDay("2020-03-28"), toDay("2020-04-12"), toDay("2020-05-01"), toDay("2020-05-15")]; } else if (scenario=="Three-waves" || scenario=="Three-waves+testing") { var phases = [0, toDay("2020-02-15"), toDay("2020-03-07"), toDay("2020-03-20"), toDay("2020-03-28"), toDay("2020-04-12"), toDay("2020-05-01"), toDay("2020-05-15"), toDay("2020-06-01"), toDay("2020-06-15")]; } return phases; }54:T455,rho = { if (scenario=="Baseline" || scenario=="5/9/2020") { var phases = [0.0, 0.0, 0.01, 0.1, 0.1, 0.1]; } else if (scenario=="6/18/2020") { var phases = [0.0, 0.0, 0.01, 0.1, 0.1, 0.1, 0.15, 0.20]; } else if (scenario=="7/28/2020") { var phases = [0.0, 0.0, 0.01, 0.1, 0.1, 0.1, 0.15, 0.20, 0.20]; } else if (scenario=="2/28/2021") { var phases = [0.0, 0.0, 0.01, 0.1, 0.1, 0.1, 0.15, 0.20, 0.20, 0.2, 0.2]; // // The remaining are hypothetical // } else if (scenario=="Two-waves") { var phases = [0.0, 0.0, 0.01, 0.1, 0.1, 0.1, 0.15, 0.15]; } else if (scenario=="Three-waves") { var phases = [0.0, 0.0, 0.01, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]; } else if (scenario=="Three-waves+testing") { var phases = [0.0, 0.0, 0.01, 0.1, 0.1, 0.1, 0.15, 0.15, 0.2, 0.2]; } return phases; }55:T1269,class SEIR { constructor(P_=P,Io_=Io,alpha_=ALPHA,gamma_=GAMMA,c_=c, to_=0,tf_=tf,nt_=250,Fa_=Fa, stochastic_=stochastic,samples_=samples, k_=k,n_=n,distribution_=distribution) { this.P = P_; this.Io = Io_; this.Ro = Ro; this.tk = tk this.alpha = alpha_; this.gamma = gamma_; this._rho = rho; this.c = c_; this.Ma = Ma; this.startTime = new Date(startTime); this.to = to_; this.tf = tf_; this.nt = nt_; this.dt = (this.tf-this.to)/this.nt; this.Fa = Fa_; this.phase = 0; this.distribution = distribution_; this.n = n; if (this.distribution==="Erlang") { this.alphaDist = new Erlang(this.alpha/k,k); this.gammaDist = new Erlang(this.gamma/k,k); } else if (this.distribution==="Normal") { this.alphaDist = new Normal(this.alpha,k); this.gammaDist = new Normal(this.gamma,k); } else if (this.distribution==="LogNormal") { this.alphaDist = new LogNormal(Math.log(this.alpha)-Math.log(k)**2/2,Math.log(k)); this.gammaDist = new LogNormal(Math.log(this.gamma)-Math.log(k)**2/2,Math.log(k)); } else if (this.distribution==="Gamma") { this.alphaDist = new Gamma(this.alpha/k,k); this.gammaDist = new Gamma(this.gamma/k,k); } else if (this.distribution==="JHU") { //this.alphaDist = new Gamma(0.95,5.81); this.alphaDist = new LogNormal(1.62,0.42); this.gammaDist = new Gamma(4.620,0.41); } this.stochastic = stochastic; this.alphai; // current sample this.gammai; // current sample this.samples = samples_; this.O = []; } get beta() { return this.Ro[this.phase]/this.gammai; } get rho() { return this._rho[this.phase]; } get scaled_c() { return this.Ma[this.phase]*this.c; } get setIC() { var IC = [this.P-this.Io, 0, 0, this.Io, 0, 0, 0, 0]; return IC; } setPhase(t) { let p=0; for(let k=0; kthis.tk[k]) p=k; } this.phase=p; return p; } rhs(U,t,dUdt) { let N = 0 for (let i=0; i<6; i++) N += U[i]; const Ni = 1.0/N; // inverse of total effective population const S = U[0]; // susceptable population const E = U[1]; // exposed population const Ih = U[2]; // Infected population in hospital (issolation) const Ic = U[3]; // Infected population in community const Rh = U[4]; // Recovered population in hospital (issolation) const Rc = U[5]; // Recovered population in community this.setPhase(t); // set the phase of the outbreak const alphat = 1.0/this.alphai; // incubation rate const gammat = 1.0/this.gammai; // infectious rate const one = 1.0; const zero = 0.0; const c = this.scaled_c; // set the current value of mortality dUdt[0] = -this.beta*Ic*S*Ni; dUdt[1] = this.beta*Ic*S*Ni - alphat*E; dUdt[2] = alphat*this.rho*E - gammat*Ih; dUdt[3] = alphat*(one-this.rho)*E - gammat*Ic; dUdt[4] = gammat*(one-c)*Ih; dUdt[5] = gammat*(one-c)*Ic; dUdt[6] = gammat*c*Ih; dUdt[7] = gammat*c*Ic; } stochastic_solve() { var U = [[0.,0.,0.,0.,0.,0.,0.,0.]]; var t = [0]; var ph = [0]; U[0] = this.setIC; let rhs = this.rhs.bind(this); for (let i=0; i