2b:null 2c:["$","title",null,{"children":"Interactive SEIR Model of NM | S. Scott Collis | Observable"}] 31:T5f2,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. The set of ordinary differential equations is given by`32: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 `33:Td0c,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 New Mexico, and the simulation runs for 8 months (240 days) by default. I've tuned the baseline parameters to match the NM 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 state * There is no direct modeling of the influence of age (everything is averaged over age) * There is no spatial dependancy. In other words, things are averaged over space and there are no net inflows or outflows. * 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`} as discrete functions of time. * The default model parameters do a decent job of matching the observed data (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. * It is also likely that mortality is not only a function of age but also time, as care and treatments improve. On 7/28 I added a mortality adjustment ${tex`M_a`} to model this resulting is much better data fits. ### 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.`34:T901,{ const data = vl.data({values:Urange,name:"our_data"}); 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 = 4000; //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 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 UScases = USdata .encode( vl.x().fieldT("date").axis({tickCount:"week",title:"Date"}), vl.y().fieldQ("dcases"), vl.color().value("green") ); const USdeaths = USdata .encode( vl.x().fieldT("date").axis({tickCount:"week",title:"Date"}), vl.y().fieldQ("ddeaths"), vl.color().value("red") ); 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) ); 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") return plot.render(); }35:Ta40,CrAndDeaths = { const ymin = 1; const ymax = 10e5; let yscale = {type:"log",base:10,domain:[ymin,ymax]}; if (ylinear) yscale = {}; 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({tooltip:true,filled:false}).data(usData) //.transform({filter:{field:"state",oneOf:["New Mexico"]}}) .encode( vl.x().fieldT("date"), vl.y().fieldQ("cases").scale(yscale), vl.color().value("green") ); const USdeaths = vl.markSquare({tooltip:true,filled:false}).data(usData) .transform({filter:{and:[{field:"state",oneOf:["New Mexico"]},{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(); }36:T4c0,tk = { if (scenario=="Baseline" || scenario=="Update 5/7/2020") { var phases = [toDay(startTime), toDay("2020-03-04"), toDay("2020-03-13"), toDay("2020-03-16"), toDay("2020-03-24"), toDay("2020-04-12")]; } else if (scenario=="Update 5/21/2020") { var phases = [toDay(startTime), toDay("2020-03-04"), toDay("2020-03-13"), toDay("2020-03-16"), toDay("2020-03-24"), toDay("2020-04-12"), toDay("2020-05-16")]; } else if (scenario=="Update 7/28/2020") { var phases = [toDay(startTime), toDay("2020-03-04"), toDay("2020-03-13"), toDay("2020-03-16"), toDay("2020-03-24"), toDay("2020-04-12"), toDay("2020-05-16"), toDay("2020-05-31")]; } else if (scenario=="Update 2/28/2021") { var phases = [toDay(startTime), toDay("2020-03-04"), toDay("2020-03-13"), toDay("2020-03-16"), toDay("2020-03-24"), toDay("2020-04-12"), toDay("2020-05-16"), toDay("2020-05-31"), toDay("2020-07-04"), toDay("2020-09-01"), toDay("2020-11-20")]; } return phases; }37:T1210,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.Ma = Ma; this.c = c_; 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; 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