Identifying the measurements required to estimate rates of COVID-19 transmission, infection, and detection, using variational data assimilation
IIdentifying the measurements required to estimate rates of COVID-19transmission, infection, and detection, using variational data assimilation
Eve Armstrong ∗ , Manuela Runge , and Jaline Gerardin Department of Physics, New York Institute of Technology, New York, NY 10023, USA Department of Astrophysics, American Museum of Natural History, New York, NY 10024, USA Department of Preventive Medicine, Northwestern University, Chicago IL, 60611
August 4, 2020
Abstract
We demonstrate the ability of statistical data assimilation to identify the measurements required for accuratestate and parameter estimation in an epidemiological model for the novel coronavirus disease COVID-19. Ourcontext is an effort to inform policy regarding social behavior, to mitigate strain on hospital capacity. Themodel unknowns are taken to be: the time-varying transmission rate, the fraction of exposed cases that requirehospitalization, and the time-varying detection probabilities of new asymptomatic and symptomatic cases. Insimulations, we obtain accurate estimates of undetected (that is, unmeasured) infectious populations, by measuringthe detected cases together with the recovered and dead - and without assumed knowledge of the detection rates.Given a noiseless measurement of the recovered population, excellent estimates of all quantities are obtained usinga temporal baseline of 101 days, with the exception of the time-varying transmission rate at times prior to theimplementation of social distancing. With low noise added to the recovered population, accurate state estimatesrequire a lengthening of the temporal baseline of measurements. Estimates of all parameters are sensitive to thecontamination, highlighting the need for accurate and uniform methods of reporting. The aim of this paper isto exemplify the power of SDA to determine what properties of measurements will yield estimates of unknownparameters to a desired precision, in a model with the complexity required to capture important features of theCOVID-19 pandemic.
I. INTRODUCTION
The coronavirus disease 2019 (COVID-19) is burdeninghealth care systems worldwide, threatening physical andpsychological health, and devastating the global economy.Individual countries and states are tasked with balancingpopulation-level mitigation measures with maintainingeconomic activity. Mathematical modeling has been usedto aid policymakers’ plans for hospital capacity needs,and to understand the minimum criteria for effectivecontact tracing [1]. Both state-level decision-making andaccurate modeling benefit from quality surveillance data.Insufficient insufficient testing capacity, however, espe-cially at the beginning of the epidemic in the UnitedStates, and other data reporting issues have meant thatsurveillance data on COVID-19 is biased and incom-plete [2–4]. Models intended to guide intervention policymust be able to handle imperfect data.Within this context, we seek a means to quantify what data must be recorded in order to estimate specific un-known quantities in an epidemiological model of COVID-19 transmission. These unknown quantities are: i) thetransmission rate, ii) the fraction of the exposed popula-tion that acquires symptoms sufficiently severe to requirehospitalization, and iii) time-varying detection probabil-ities of asymptomatic and symptomatic cases. In thispaper, we demonstrate the ability of statistical data as-similation (SDA) to quantify the accuracy to which theseparameters can be estimated, given certain properties ofthe data including noise level.SDA is an inverse formulation [5]: a machine learningapproach designed to optimally combine a model withdata. Invented for numerical weather prediction [6–11],and more recently applied to biological neuron mod-els [12–18], SDA offers a systematic means to identifythe measurements required to estimate unknown modelparameters to a desired precision.Data assimilation has been presented as a means ∗ [email protected] a r X i v : . [ q - b i o . P E ] J u l or general epidemiological forecasting [19], and onework has examined variational data assimilation specif-ically - the method we employ in this paper - for esti-mating parameters in epidemiological models [20]. Re-lated Bayesian frameworks for estimating unknown prop-erties of epidemiological models have also been ex-plored [21,22]. To date, there have been two employmentsof SDA for COVID-19 specifically. Ref [23] used a simpleSIR (susceptible/infected/recovered) model, and Ref [24]expanded the SIR model to include a compartment ofpatients in treatment.Two features of our work distinguish this paper asnovel. First, we expand the model in terms of the num-ber of compartments. The aim here is to capture keyfeatures of COVID-19 such that the model structure is rel-evant for questions from policymakers on containing thepandemic. These features are: i) asymptomatic, presymp-tomatic, and symptomatic populations, ii) undetectedand detected cases, and iii) two hospitalized populations:those who do and do not require critical care. For ourmotivations for these choices, see Model . Second, weemploy SDA for the specific purpose of examining thesensitivity of estimates of time-varying parameters tovarious properties of the measurements, including thedegree of noise (or error) added. Moreover, we aim todemonstrate the power and versatility of the SDA tech-nique to explore what is required of measurements tocomplete a model with a dimension sufficiently high tocapture the policy-relevant complexities of COVID-19transmission and containment - an examination that hasnot previously been done.To this end, we sought to estimate the parame-ters noted above, using simulated data representing ametropolitan-area population loosely based on New YorkCity. We examined the sensitivity of estimations to: i)the subpopulations that were sampled, ii) the temporalbaseline of sampling, and iii) uncertainty in the sampling.Results using simulated data are threefold. First, rea-sonable estimations of time-varying detection probabili-ties require the reporting of new detected cases (asymp-tomatic and symptomatic), dead, and recovered. Second,given noiseless measurements, a temporal baseline of101 days is sufficient for the SDA procedure to capturethe general trends in the evolution of the model popula-tions, the detection probabilities, and the time-varyingtransmission rate following the implementation of socialdistancing. Importantly, the information contained in themeasured detected populations propagates successfully tothe estimation of the numbers of severe undetected cases. Third, the state evolution - and importantly the popu-lations requiring inpatient care - tolerates low ( ∼ fivepercent) noise, given a doubling of the temporal baselineof measurements; the parameter estimates do not toleratethis contamination.Finally, we discuss necessary modifications prior totesting with real data, including lowering the sensitivityof parameter estimates to noise in data. II. MODEL
The model is written in 22 state variables, each repre-senting a subpopulation of people; the total populationis conserved. Figure 1 shows a schematic of the modelstructure. Each member of a Population S that becomesExposed ( E ) ultimately reaches either the Recovered ( R )or Dead ( D ) state. Absent additive noise, the model is deter-ministic.
Five variables correspond to measured quantitiesin the inference experiments.As noted, the model is written with the aim to in-form policy on social behavior and contact tracing so asto avoid exceeding hospital capacity. To this end, themodel resolves asymptomatic-versus-symptomatic cases,undetected-versus-detected cases, and the two tiers ofhospital needs: the general (inpatient, non-intensive careunit (ICU)) H versus the critical care (ICU) C populations. Figure 1:
Schematic of the model . Each rectangle representsa population. Note the distinction of asymptomatic cases, un-detected cases, and the two tiers of hospitalized care: H and C . The aim of including this degree of resolution is to informpolicy on social behavior so as to minimize strain on hospitalcapacity. The red ovals indicate the variables that correspondto measured quantities in the inference experiments. H and C populations because hos-pital inpatient and ICU bed capacities are the key healthsystem metrics that we aim to avoid straining. Any policythat we consider must include predictions on inpatientand ICU bed needs. Preparing for those needs is a keyresponse if or when the epidemic grows uncontrolled.For details of the model, including the reaction equa-tions and descriptions of all state variables and parame-ters, see Appendix A . III. METHOD A. General inference formulation
SDA is an inference procedure, or a type of machinelearning, in which a model dynamical system is assumedto underlie any measured quantities. This model F canbe written as a set of D ordinary differential equationsthat evolve in some parameterization t as:d x a ( t ) d t = F a ( x ( t ) , p ( t )) ; a =
1, 2, . . . , D ,where the components x a of the vector x are the modelstate variables, and unknown parameters to be estimatedare contained in p ( t ) . A subset L of the D state vari-ables is associated with measured quantities. One seeksto estimate the p unknown parameters and the evolu-tion of all state variables that is consistent with the L measurements.A prerequisite for estimation using real data is thedesign of simulated experiments, wherein the true valuesof parameters are known. In addition to providing a consistency check, simulated experiments offer the op-portunity to ascertain which and how few experimentalmeasurements, in principle, are necessary and sufficientto complete a model.B. Optimization framework
SDA can be formulated as an optimization, wherein a costfunction is extremized. We take this approach, and writethe cost function in two terms: 1) one term representingthe difference between state estimate and measurement(measurement error), and 2) a term representing modelerror. It will be shown below in this Section that treat-ing the model error as finite offers a means to identifywhether a solution has been found within a particularregion of parameter space. This is a non-trivial problem,as any nonlinear model will render the cost function non-convex. We search the surface of the cost function via thevariational method, and we employ a method of anneal-ing to identify a lowest minumum - a procedure that hasbeen referred to loosely in the literature as variationalannealing (VA).The cost function A used in this paper is written as: A ( x ( n ) , p ) = J ∑ j = L ∑ l = R lm ( y l ( n ) − x l ( n )) + N − ∑ n = D ∑ a = R af ( x a ( n + ) − f a ( x ( n ) , p ( n ))) .(1)One seeks the path X = x ( ) , ..., x ( N ) , p ( ) , ... p ( N ) instate space on which A attains a minimum value . Notethat Equation 1 is shorthand; for the full form, see Ap-pendix A of Ref [18]. For a derivation - beginning with thephysical Action of a particle in state space - see Ref [25].The first squared term of Equation 1 governs thetransfer of information from measurements y l to modelstates x l . The summation on j runs over all discretizedtimepoints J at which measurements are made, whichmay be a subset of all integrated model timepoints. Thesummation on l is taken over all L measured quantities.The second squared term of Equation 1 incorpo-rates the model evolution of all D state variables x a .The term f a ( x ( n )) is defined, for discretization, as: [ F a ( x ( n )) + F a ( x ( n + ))] . The outer sum on n is takenover all discretized timepoints of the model equations ofmotion. The sum on a is taken over all D state variables. It may interest the reader that one can derive this cost function by considering the classical physical Action on a path in a state space,where the path of lowest Action corresponds to the correct solution [25] m and R f are inverse covariance matrices for themeasurement and model errors, respectively. We takeeach matrix to be diagonal and treat them as relativeweighting terms, whose utility will be described belowin this Section.The procedure searches a ( D ( N + ) + p ( N + )) -dimensional state space, where D is the number of statevariables, N is the number of discretized steps, and p isthe number of unknown parameters. To perform simu-lated experiments, the equations of motion are integratedforward to yield simulated data, and the VA procedureis challenged to infer the parameters and the evolutionof all state variables - measured and unmeasured - thatgenerated the simulated data.This specific formulation has been tested with chaoticmodels [26–29], and used to estimate parameters in mod-els of biological neurons [13, 14, 16, 18, 30, 31], as well asastrophysical scenarios [32].C. Annealing to identify a solution on a non-convexcost function surface
Our model is nonlinear, and thus the cost function sur-face is non-convex. For this reason, we iterate - or anneal- in terms of the ratio of model and measurement error,with the aim to gradually freeze out a lowest minimum.This procedure was introduced in Ref [33], and has sincebeen used in combination with variational optimizationon nonlinear models in Refs [11, 18, 30, 32] above. Theannealing works as follows.We first define the coefficient of measurement error R m to be 1.0, and write the coefficient of model error R f as: R f = R f ,0 α β , where R f ,0 is a small number nearzero, α is a small number greater than 1.0, and β is ini-tialized at zero. Parameter β is our annealing parameter.For the case in which β =
0, relatively free from modelconstraints the cost function surface is smooth and there exists one minimum of the variational problem that isconsistent with the measurements. We obtain an estimateof that minimum. Then we increase the weight of themodel term slightly, via an integer increment in β , andrecalculate the cost. We do this recursively, toward thedeterministic limit of R f (cid:29) R m . The aim is to remainsufficiently near to the lowest minimum to not becometrapped in a local minimum as the surface becomes re-solved. We will show in Results that a plot of the cost as afunction of β reveals whether a solution has been foundthat is consistent with both measurements and model. IV. THE EXPERIMENTS A. Simulated experiments
We based our simulated locality loosely on New YorkCity, with a population of 9 million. For simplicity, weassume a closed population. Simulations ran from aninitial time t of four days prior to 2020 March 1, the dateof the first reported COVID-19 case in New York City [34].At time t , there existed one detected symptomatic casewithin the population of 9 million. In addition to thatone initial detected case, we took as our initial conditionson the populations to be: 50 undetected asymptomatics,10 undetected mild symptomatics, and one undetectedsevere symptomatic.We chose five quantities as unknown parameters to beestimated (Table 1): 1) the time-varying transmission rate K i (t); 2) the detection probability of mild symptomaticcases d Sym ( t ) , 3) the detection probability of severe symp-tomatic cases d Sys ( t ) , 4) the fraction of cases that becomesymptomatic f sympt , and 5) the fraction of symptomaticcases that become severe enough to require hospitaliza-tion f severe . Here we summarize the key features that wesought to capture in modeling these parameters; for theirmathematical formulatons, see Appendix B . Table 1:
Unknown parameters to be estimated. K i , d Sym , and d Sys are taken to be time-varying. Parameters f sympt and f severe are constant numbers, as they are assumed to reflect an intrinsic property of the disease. The detection probability ofasymptomatic cases is taken to be known and zero. Parameter DescriptionK i ( t ) Time-varying transmission rate d Sym ( t ) Time-varying detection probability of mild symptomatics d Sys ( t ) Time-varying detection probability of symptomatics requiring hospitalization f sympt Fraction of positive cases that produce symptoms f severe Fraction of symptomatics that are severe K i (often referred to as the ef-fective contact rate) in a given population for a giveninfectious disease is measured in effective contacts perunit time. This may be expressed as the total contact ratemultiplied by the risk of infection, given contact betweenan infectious and a susceptible individual. The contactrate, in turn, can be impacted by amendments to socialbehavior .As a first step in applying SDA to a high-dimensionalepidemiological model, we chose to condense the signif-icance of K i into a relatively simple mathematical form.We assumed that K i was constant prior to the imple-mentation of a social-distancing mandate, which theneffected a rapid transition of K i to a lower constant value.Specifically, we modeled K i as a smooth approximationto a Heaviside function that begins its decline on March22, the date that the stay-at-home order took effect inNew York City [38]: 25 days after time t . For furthersimplicity, we took K i to reflect a single implementationof a social distancing protocol, and adherence to thatprotocol throughout the remaining temporal baseline ofestimation.Detection rates impact the sizes of the subpopulationsentering hospitals, and their values are highly uncer-tain [3, 4]. Thus we took these quantities to be unknown,and - as detection methods will evolve - time-varying.We also optimistically assumed that the methods willimprove, and thus we described them as increasing func-tions of time. We used smoothly-varying forms, the firstlinear and the second quadratic, to preclude symmetriesin the model equations. Meanwhile, we took the detec-tion probability for asymptomatic cases ( d As ) to be knownand zero, a reasonable reflection of the current state oftesting.Finally, we assigned as unknowns the fraction of casesthat become symptomatic ( f sympt ) and fraction of symp-tomatic cases that become sufficiently severe to requirehospitalization ( f severe ), as these fractions possess highuncertainties (Refs [39] and [40], respectively). As they re-flect an intrinsic property of the disease, we took them tobe constants. All other model parameters were taken tobe known and constant ( Appendix A ); however, the valuesof many other model parameters also possess significantuncertainties given the reported data, including, for ex-ample, the fraction of those hospitalized that require ICUcare. Future VA experiments can treat these quantities asunknowns as well.
Figure 2:
Schematic of the four simulated experiments.
The simulated experiments are summarized in theschematic of Figure 2. They were designed to probe theeffects upon estimations of three considerations: a) thenumber of measured subpopulations, b) the temporalbaseline of measurements, and c) contamination of mea-surements by noise. To this end, we designed a “base”experiment sufficient to yield an excellent solution, andthen four variations on this experiment.The base experiment ( i in Figure 2) possesses the fol-lowing features: a) five measured populations: detectedasymptomatic As det , detected mild symptomatic Sym det ,detected severe symptomatic
Sys det , Recovered R , andDead D ; b) a temporal baseline of 101 days, beginningon 2020 February 26; c) no noise in measurements.The three variations on this basic experiment ( ii through iv in Figure 2), incorporate the following in-dependent changes. In Experiment ii , the R population isnot measured - an example designed to reflect the currentsituation in some localities (e.g. Refs [3, 4]).Experiment iii includes a ∼ five percent noise level(for the form of additive noise, see Appendix C ) in thesimulated R data, and Experiment iv includes that noiselevel in addition to a doubled temporal baseline.For each experiment, twenty independent calculationswere initiated in parallel searches, each with a randomly-generated set of initial conditions on state variable andparameter values. For technical details of all experimentaldesigns and implementation, see Appendix C . V. RESULT A. General findings
The salient results for the simulated experiments i through iv are as follows: The reproduction number R , in the simplest SIR form, can be written as the effective contact rate divided by the recovery rate. Inpractice, R is a challenge to infer [22, 35–37]. K i (t) at times prior to the onset ofsocial distancing;ii (absent a measurement of Population R ): Poor esti-mate of all quantities;iii ( ∼
5% additive noise in R ): Poor estimates of allquantities;iv ( ∼
5% additive noise in R , with a doubled baselineof 201 days): Estimates of state evolution are robustto noise, while parameter estimates are sensitive tonoise.Figures of the estimated time evolution of state vari-ables and time-varying parameters are shown in theirrespective subsections, and the estimates of the staticparameters are listed in Table 2.B. Base Experiment i
The base experiment that employed five noiseless mea-sured populations over 101 days yielded an excellentsolution in terms of model evolution and parameter es-timates. Prior to examining the solution, we shall firstshow the cost function versus the annealing parameter β ,as this distribution can serve as a tool for assessing thesignificance of a solution. Table 2:
Estimates of static parameters f sympt and f severe overall simulated experiments. For Experiments i and iv , the re-ported numbers are taken from the annealing iteration witha value of parameter β of 32 and 40, respectively: once thedeterministic limit has been reached (see text). For Experiment ii , an attempt was made to retrieve parameter estimates at β =
2; that is: before the solution grows unstable exponen-tially (see Figure 5). See specific subsections for details of eachexperiment.
Experiment f sympt (true: 0.6) f severe (true: 0.07) Mean Variance Mean Variancei × − × − ii – iii – iv Cost function plotted at each annealing step β forthe base experiment i , for twenty paths in state space, where β scales the rigidity of the imposed model constraint. At low β the procedure endeavours to fit the measured variables tothe simulated measurements. As β increases, the cost increasesuntil it approaches a plateau (around β = Figure 3 shows the evolution of the cost throughoutannealing, for the ten distinct independent paths thatwere initiated; the x-axis shows the value of AnnealingParameter β , or: the increasing rigidity of the modelconstraint. At the start of iterations, the cost function ismainly fitting the measurements to data, and its value be-gins to climb as the model penalty is gradually imposed.If the procedure finds a solution that is consistent notonly with the measurements, but also with the model,then the cost will plateau. In Figure 4, we see this hap-pen, around β =
30, with some scatter across paths. Thereported estimates in this Subsection are taken at a valueof β of 32: on the plateau. The significance of this plateauwill become clearer upon examining the contrasting caseof Experiment ii .We now examine the state and parameter estimatesfor the base experiment i . For all experiments, eachsolution shown is representative of the solution for alltwenty paths. Figure 4 shows an excellent estimate ofall state variables during the temporal window in whichthe measured variables were sampled. For consistency inillustrating the time evolution of all state variables, weuse the state estimates for the Recovered ( R ) and Dead( D ) populations, which are cumulative, rather than followstandard epidemiological practice of showing incident R or D . The time-varying parameters are also estimatedwell, excepting K i (t) at times prior to its steep decline.We noted no improvement in this estimate for K i (t), fol-lowing a tenfold increase in the temporal resolution ofmeasurements (not shown). The procedure does appear6 igure 4: Estimates of the state - measured and unmeasured - variables, and the time-varying parameters K i , d Sym , and d Sys , for the base experiment i . Excellent estimates are obtained of all states and parameters, except early values of K i prior tothe implementation of social distancing; see text. Results are taken at a value for annealing parameter β of 32. to recognize that a fast transition in the value of K i oc-curred at early times, and that that value was previouslyhigher. It will be important to investigate the reason forthis failure in the estimation of K i at early times, to ruleout numerical issues involved with the quickly-changingderivative .C. Experiment ii: no measurement of R
Figure 5 shows the cost as a function of annealing forthe case with no measurement of Recovered Population R . Without examining the estimates, we know from theCost( β ) plot that no solution has been found that is con-sistent with both measurements and model: no plateauis reached. Rather, as the model constraint strengthens,the cost increases exponentially.Indeed, Figure 6 shows the estimation, taken at β = β , we conclude that there exists insufficient information in As det , Sym det , Sys det , and D alone to corral the proce-dure into a region of state-and-parameter space in whicha model solution is possible. We repeated this experi-ment with a doubled baseline of 201 days, and noted noimprovement (not shown). Figure 5:
Cost versus β for Experiment ii : R is not measured. As β increases, the cost increases indefinitely, indicating thatno solution has been found that is consistent with both mea-surements and model dynamics. As noted in
Experiments , we chose K i to reflect a rapid adherence to social distancing at Day 25 following time t , which then remained inplace through to Day 101. For the form of K i , see Appendix B .) igure 6: Estimates for Experiment ii : without a measurement of Population R . This result is taken at β =
2, prior to theexponential runaway in the cost. Estimates of unmeasured states and time-varying parameters are poor.Figure 7:
Estimates for Experiment iv : low noise added to Population R and with a doubled temporal baseline of 201 days. The noise added to R propagates to some unmeasured States ( S , E , As , and As det ), but the overall evolution is captured well.The noise precludes an estimate of the time-varying parameters (not shown). Results are reported using a value for β of 40. Experiments iii and iv: low noise added
In Experiment iii , the low noise added to R yielded apoor state and parameter estimate (not shown). Witha doubled temporal baseline of measurements (Experi-ment iv ), however, the state estimate became robust to thecontamination. Figure 7 shows this estimate. While the ∼ five percent noise added to Population R propagatesto the unmeasured States S , E , and P , the general stateevolution is still captured well. Importantly, the popula-tions entering the hospital are well estimated. Note thatsome low state estimates (e.g. As ) are not perfectly offsetby high estimates (e.g. Sym ). The addition of noise inthese numbers - by definition - breaks the conservationof the population. Finally, the parameter estimates forExperiment iv do not survive the added contamination(not shown). VI. CONCLUSION
We have endeavoured to illustrate the potential of SDAto systematically identify the specific measurements, tem-poral baseline of measurements, and degree of measure-ment accuracy, required to estimate unknown modelparameters in a high-dimensional model designed to ex-amine the complex problems that COVID-19 presents tohospitals. In light of our assumed knowledge of somemodel parameters, we restrict our conclusions to gen-eral comments. We emphasize that estimation of the fullmodel state requires measurements of the detected casesbut not the undetected, provided that the recovered anddead are also measured. The state evolution is tolerantto low noise in these measurements, while the parameterestimates are not.The ultimate aim of SDA is to test the validity of model estimation using real data, via prediction . In ad-vance of that step, we are performing a detailed study ofthe model’s sensitivity to contamination in the measur-able populations As det , Sym det , Sys det , R , and D . Concur-rently we are examining means to render the parameterestimation less sensitive to noise, via various additionalequality constraints in the cost function, and looseningthe assumption of Gaussian-distributed noise. In partic-ular, we shall require that the time-varying parametersbe smoothly-varying. It will be important to examine thestability of the SDA procedure over a range of choicesfor parameter values and initial numbers for the infectedpopulations.This procedure can be expanded in many directions.Examples include: 1) defining additional model parame-ters as unknowns to be estimated, including the fractionof patients hospitalized, the fraction who enter criticalcare, and the various timescales governing the reactionequations; 2) imposing various constraints regarding theunknown time-varying quantities, particularly transmis-sion rate K i (t), and identifying which forms permit asolution consistent with measurements; 3) examiningmodel sensitivity to the initial numbers within each pop-ulation; 4) examining model sensitivity to the temporalfrequency of data sampling. Moreover, it is our hopethat the exercises described in this paper can guide theapplication of SDA to a host of complicated questionssurrounding COVID-19. VII. ACKNOWLEDGEMENTS
Thank you to Patrick Clay from the University of Michi-gan for discussions on inferring exposure rates givensocial distancing protocols. 9 ppendix A: Details of the model
Table 3:
State variables of the COVID-19 transmission model.
The “detected” qualifier signifies that the population has beentested and is positive for COVID-19.
Variable DescriptionS
Susceptible E Exposed As det Asymptomatic, detected As Asymptomatic, undetected
Sym det
Symptomatic mild, detected
Sym
Symptomatic mild, undetected
Sys det
Symptomatic severe, detected
Sys
Symptomatic severe, undetected H det Hospitalized and will recover, detected H det Hospitalized and will go to critical care and recover, detected H det Hospitalized and will go to critical care and die, detected H Hospitalized and will recover, undetected H Hospitalized and will go to critical care and recover, undetected H Hospitalized and will go to critical care and die, undetected C det In critical care and will recover, detected C det In critical care and will die, detected C In critical care and will recover, undetected C In critical care and will die, undetected R Recovered D Dead
Reaction equations
The blue notation specified by overbrackets denotes the correspondence of specific terms to the reactions betweenthe populations depicted in Figure 1.d S d t = − → E (cid:122) (cid:125)(cid:124) (cid:123) K i · S · [ in f ectious + ( in f ectious det × reduced )] N • in f ectious = As + P + Sym + Sys + H + H + H + C + C • in f ectious det = As det + Sym det + Sys det d E d t = S → (cid:122) (cid:125)(cid:124) (cid:123) K i · S · [ in f ectious + ( in f ectious det × reduced )] / N − → As det (cid:122) (cid:125)(cid:124) (cid:123) − f sympt t in f ection · E · d As − → As (cid:122) (cid:125)(cid:124) (cid:123) − f sympt t in f ection · E · ( − d As ) − → P (cid:122) (cid:125)(cid:124) (cid:123) f sympt t in f ection · E d As det d t = E → (cid:122) (cid:125)(cid:124) (cid:123) − f sympt t in f ection · E · d As − → R (cid:122) (cid:125)(cid:124) (cid:123) t R , a · As det d As d t = E → (cid:122) (cid:125)(cid:124) (cid:123) − f sympt t in f ection · E · ( − d As ) − → R (cid:122) (cid:125)(cid:124) (cid:123) t R , a · As P d t = E → (cid:122) (cid:125)(cid:124) (cid:123) f sympt t in f ection · E − → Sym det (cid:122) (cid:125)(cid:124) (cid:123) − f severe t sympt · P · d Sym − → Sym (cid:122) (cid:125)(cid:124) (cid:123) − f severe t sympt · P · ( − d Sym ) − → Sys det (cid:122) (cid:125)(cid:124) (cid:123) f severe t sympt · P · d Sys − → Sys (cid:122) (cid:125)(cid:124) (cid:123) f severe t sympt · P · ( − d Sys ) d Sym det d t = P → (cid:122) (cid:125)(cid:124) (cid:123) − f severe t sympt · P · d Sym − → R (cid:122) (cid:125)(cid:124) (cid:123) t R , m · Sym det d Sym d t = P → (cid:122) (cid:125)(cid:124) (cid:123) − f severe t sympt · P · ( − d Sym ) − → R (cid:122) (cid:125)(cid:124) (cid:123) t R , m · Sym d Sys det d t = P → (cid:122) (cid:125)(cid:124) (cid:123) f severe t sympt · P · d Sys − → H det (cid:122) (cid:125)(cid:124) (cid:123) f H t H · Sys det − → H det (cid:122) (cid:125)(cid:124) (cid:123) f C t H · Sys det − → H det (cid:122) (cid:125)(cid:124) (cid:123) f D t H · Sys det d Sys d t = P → (cid:122) (cid:125)(cid:124) (cid:123) f severe t sympt · P · ( − d Sys ) − → H (cid:122) (cid:125)(cid:124) (cid:123) f H t H · Sys − → H (cid:122) (cid:125)(cid:124) (cid:123) f C t H · Sys − → H (cid:122) (cid:125)(cid:124) (cid:123) f D t H · Sys d H det d t = Sys det → (cid:122) (cid:125)(cid:124) (cid:123) f H t H · Sys det − → R (cid:122) (cid:125)(cid:124) (cid:123) t R , h · H det d H det d t = Sys det → (cid:122) (cid:125)(cid:124) (cid:123) f C t H · Sys det − → C det (cid:122) (cid:125)(cid:124) (cid:123) t C · H det d H det d t = Sys det → (cid:122) (cid:125)(cid:124) (cid:123) f D t H · Sys det − → C det (cid:122) (cid:125)(cid:124) (cid:123) t C · H det d H d t = Sys → (cid:122) (cid:125)(cid:124) (cid:123) f H t H · Sys − → R (cid:122) (cid:125)(cid:124) (cid:123) t R , h · H d H d t = Sys → (cid:122) (cid:125)(cid:124) (cid:123) f C t H · Sys − → C (cid:122) (cid:125)(cid:124) (cid:123) t C · H d H d t = Sys → (cid:122) (cid:125)(cid:124) (cid:123) f D t H · Sys − → C (cid:122) (cid:125)(cid:124) (cid:123) t C · H d C det d t = H det → (cid:122) (cid:125)(cid:124) (cid:123) t C · H det − → R (cid:122) (cid:125)(cid:124) (cid:123) t R , c · C det d C det d t = H det → (cid:122) (cid:125)(cid:124) (cid:123) t C · H det − → D (cid:122) (cid:125)(cid:124) (cid:123) t D · C det d C d t = H → (cid:122) (cid:125)(cid:124) (cid:123) t C · H − → R (cid:122) (cid:125)(cid:124) (cid:123) t R , c · C d C d t = H → (cid:122) (cid:125)(cid:124) (cid:123) t C · H − → D (cid:122) (cid:125)(cid:124) (cid:123) t D · C R d t = As det → (cid:122) (cid:125)(cid:124) (cid:123) t R , a · As det + As → (cid:122) (cid:125)(cid:124) (cid:123) t R , a · As + Sym det → (cid:122) (cid:125)(cid:124) (cid:123) t R , m · Sym det + Sym → (cid:122) (cid:125)(cid:124) (cid:123) t R , m · Sym + H det → (cid:122) (cid:125)(cid:124) (cid:123) t R , h · H det + H → (cid:122) (cid:125)(cid:124) (cid:123) t R , h · H + C det → (cid:122) (cid:125)(cid:124) (cid:123) t R , c · C det + C → (cid:122) (cid:125)(cid:124) (cid:123) t R , c · C d D d t = C det → (cid:122) (cid:125)(cid:124) (cid:123) t D · C det + C → (cid:122) (cid:125)(cid:124) (cid:123) t D · C Table 4:
The model parameters, with the unknown parameters to be estimated denoted in boldface.
The unknown param-eters K i , Sym , and d Sys are taken to be time-varying. The unknown parameters f sympt and f severe are taken to be intrinsicproperties of the disease and therefore constant numbers. The detection probability of asymptomatic cases is taken to be knownand zero. Units of time are days. Parameter Description ValueN
Total population 9,000,000 reduced
The property that a detected case is likely to transmit less, via successfulquarantine) 0.2 K i ( t ) Transmission rate
See
Appendix Bd As ( t ) Detection probability of asymptomatic cases 0.0 f sympt Fraction of positive cases that produce symptoms 0.6 [39] t in f ection Time from exposure to infection 4.0 [41] t R , a Time to recovery for asymptomatics 8.0
Assumed to besame as t R , m d Sym ( t ) Detection probability of mild symptomatics
See
Appendix B d Sys ( t ) Detection probability of severe symptomatics
See
Appendix B f severe Fraction of symptomatics that are severe 0.07 [40] t sympt Time to symptoms, for symptomatics 4.0 [41, 42] t R , m Time from symptoms to recovery, for mild symptomatics 8.0 [43] a f H Fraction of severe cases that are hospitalized and then recover: f H = − f C − f D f C Fraction of severe cases that require critical care and then recover 0.3 [44] f D Fraction of severe cases that die 0.04 [45] t H Time from symptoms to hospital, for severe symptomatics 5.0 [46] t R , h Time from entering hospital to recovery, for severe symptomatics thatdo not require critical care 10.0 [44, 45] t C Time from entering hospital to critical care, for severe symptomatics 5.0 [46] t R , c Time from entering critical care to recovery for severe symptomatics 10.0 [47] t D Time from entering critical care to death, for severe symptomatics 5.0 [48] a As described in [43], viral load can be high and detectable for up to 20 days. We choose a shorter duration of infectiousness to capturethe time during which transmissibility is highest.
Appendix B: Unknown time-varying parameters to be estimated
The unknown parameters assumed to be time-varying are the transmission rate K i , and the detection probabilities d Sym and d Sys for mild and severe symptomatic cases, respectively.The transmission rate in a given population for a given infectious disease is measured in effective contacts perunit time. This may be expressed as the total contact rate (the total number of contacts, effective or not, per unittime), multiplied by the risk of infection, given contact between an infectious and a susceptible individual. The totalcontact rate can be impacted by social behavior.In this first employment of SDA upon a pandemic model of such high dimensionality, we chose to represent K i asa relatively constant value that undergoes one rapid transition corresponding to a single social distancing mandate.As noted in Experiments , social distancing rules were imposed in New York City roughly 25 days following the first12eported case. We thus chose K i to transition between two relatively constant levels, roughly 25 days following time t . Specifically, we wrote K i (t) as: Ki ( t ) = − f · e ( T − t ) / s + + ξ .The parameter T was set to 25, beginning four days prior to the first report of a detection in NYC [34] to the impositionof a stay-home order in NYC on March 22 [38]. The parameter s governs the steepness of the transformation, andwas set to 10. Parameters f and ξ were then adjusted to 1.2 and 1.5, to achieve a transition from about 1.4 to 0.3.For detection probabilities d Sym and d Sys , a linear and quadratic form, respectively, were chosen to precludesymmetries, and both were optimistically taken to increase with time: dSym ( t ) = · tdSys ( t ) = · t Finally, each time series was normalized to the range: [0:1], via division by their respective maximum values.
Appendix C: Technical details of the inference experiments
The simulated data were generated by integrating the reaction equations (
Appendix A ) via a fourth-order adaptiveRunge-Kutta method encoded in the Python package odeINT. A step size of one (day) was used to record theoutput. Except for the one instance noted in
Results regarding Experiment i , we did not examine the sensitivity ofestimations to the temporal sparsity of measurements. The initial conditions on the populations were: S = N − N is the total population), As =
1, and zero for all others.For the noise experiments, the noise added to the simulated
Sym det , Sys det , and R data were generated by Python’s numpy.random.normal package, which defines a normal distribution of noise. For the “low-noise” experiments, weset the standard deviation to be the respective mean of each distribution, divided by 100. For the experiments usinghigher noise, we multiplied that original level by a factor of ten. For each noisy data set, the absolute value of theminimum was then added to each data point, so that the population did not drop below zero.The optimization was performed via the open-source Interior-point Optimizer (Ipopt) [49]. Ipopt uses aSimpsonâ ˘A ´Zs rule method of finite differences to discretize the state space, a Newtonâ ˘A ´Zs method to search, anda barrier method to impose user-defined bounds that are placed upon the searches. We note that Ipopt’s searchalgorithm treats state variables as independent quantities, which is not the case for a model involving a closedpopulation. This feature did not affect the results of this paper. Those interested in expanding the use of this tool,however, might keep in mind this feature. One might negate undesired effects by, for example, imposing equalityconstraints into the cost function that enforce the conservation of N .Within the annealing procedure described in Methods , the parameter α was set to 2.0, and β ran from 0 to 38 inincrements of 1. The inverse covariance matrix for measurement error ( R m ) was set to 1.0, and the initial value of theinverse covariance matrix for model error ( R f ,0 ) was set to 10 − .For each of the four simulated experiments, twenty paths were searched, beginning at randomly-generated initialconditions for parameters and state variables. All simulations were run on a 720-core, 1440-GB, 64-bit CPU cluster. References [1] IHME COVID-19 health service utilization forecasting Team and Christopher JL Murray. Forecasting COVID-19impact on hospital bed-days, ICU-days, ventilator-days and deaths by US state in the next 4 months. medRxiv ,page 2020.03.27.20043752, March 2020. Publisher: Cold Spring Harbor Laboratory Press.[2] Misty Heggeness. The need for data innovation in the time of covid-19. . Accessed: 2020-05-17. 133] Daniel Weinberger, Ted Cohen, Forrest Crawford, Farzad Mostashari, Don Olson, Virginia E Pitzer, Nicholas GReich, Marcus Russi, Lone Simonsen, Annie Watkins, et al. Estimating the early death toll of covid-19 in theunited states.
Medrxiv , 2020.[4] Ruiyun Li, Sen Pei, Bin Chen, Yimeng Song, Tao Zhang, Wan Yang, and Jeffrey Shaman. Substantial undocu-mented infection facilitates the rapid dissemination of novel coronavirus (sars-cov-2).
Science , 368(6490):489–493,2020.[5] Albert Tarantola.
Inverse problem theory and methods for model parameter estimation . SIAM, 2005.[6] Ryuji Kimura. Numerical weather prediction.
Journal of Wind Engineering and Industrial Aerodynamics , 90(12-15):1403–1414, 2002.[7] Eugenia Kalnay.
Atmospheric modeling, data assimilation and predictability . Cambridge university press, 2003.[8] Geir Evensen.
Data assimilation: the ensemble Kalman filter . Springer Science & Business Media, 2009.[9] John T Betts.
Practical methods for optimal control and estimation using nonlinear programming , volume 19. Siam,2010.[10] William G Whartenby, John C Quinn, and Henry DI Abarbanel. The number of required observations in dataassimilation for a shallow-water flow.
Monthly Weather Review , 141(7):2502–2518, 2013.[11] Zhe An, Daniel Rey, Jingxin Ye, and Henry DI Abarbanel. Estimating the state of a geophysical system withsparse observations: time delay methods to achieve accurate initial states for prediction.
Nonlinear Processes inGeophysics (Online) , 24(1), 2017.[12] Steven J Schiff. Kalman meets neuron: the emerging intersection of control theory with neuroscience. In , pages 3318–3321. IEEE,2009.[13] Bryan A Toth, Mark Kostuk, C Daniel Meliza, Daniel Margoliash, and Henry DI Abarbanel. Dynamicalestimation of neuron and network properties i: variational methods.
Biological cybernetics , 105(3-4):217–237,2011.[14] Mark Kostuk, Bryan A Toth, C Daniel Meliza, Daniel Margoliash, and Henry DI Abarbanel. Dynamicalestimation of neuron and network properties ii: path integral monte carlo methods.
Biological cybernetics ,106(3):155–167, 2012.[15] Franz Hamilton, Tyrus Berry, Nathalia Peixoto, and Timothy Sauer. Real-time tracking of neuronal networkstructure using data assimilation.
Physical Review E , 88(5):052715, 2013.[16] C Daniel Meliza, Mark Kostuk, Hao Huang, Alain Nogaret, Daniel Margoliash, and Henry DI Abarbanel.Estimating parameters and predicting membrane voltages with conductance-based neuron models.
Biologicalcybernetics , 108(4):495–516, 2014.[17] Alain Nogaret, C Daniel Meliza, Daniel Margoliash, and Henry DI Abarbanel. Automatic construction ofpredictive neuron models through large scale assimilation of electrophysiological data.
Scientific reports ,6(1):1–14, 2016.[18] Eve Armstrong. Statistical data assimilation for estimating electrophysiology simultaneously with connectivitywithin a biological neuronal network.
Physical Review E , 101(1):012415, 2020.[19] Luís MA Bettencourt, Ruy M Ribeiro, Gerardo Chowell, Timothy Lant, and Carlos Castillo-Chavez. Towardsreal time epidemiology: data assimilation, modeling and anomaly detection of health surveillance data streams.In
NSF Workshop on Intelligence and Security Informatics , pages 79–90. Springer, 2007. 1420] CJ Rhodes and T Déirdre Hollingsworth. Variational data assimilation with epidemic models.
Journal oftheoretical biology , 258(4):591–602, 2009.[21] Loren Cobb, Ashok Krishnamurthy, Jan Mandel, and Jonathan D Beezley. Bayesian tracking of emergingepidemics using ensemble optimal statistical interpolation.
Spatial and spatio-temporal epidemiology , 10:39–48,2014.[22] Luis MA Bettencourt and Ruy M Ribeiro. Real time bayesian estimation of the epidemic potential of emerginginfectious diseases.
PLoS One , 3(5), 2008.[23] Jörn Lothar Sesterhenn. Adjoint-based data assimilation of an epidemiology model for the covid-19 pandemicin 2020. arXiv preprint arXiv:2003.13071 , 2020.[24] Philip Nadler, Shuo Wang, Rossella Arcucci, Xian Yang, and Yike Guo. An epidemiological modelling approachfor covid19 via data assimilation. arXiv preprint arXiv:2004.12130 , 2020.[25] Henry Abarbanel.
Predicting the future: completing models of observed complex systems . Springer, 2013.[26] Henry DI Abarbanel, P Bryant, Philip E Gill, Mark Kostuk, Justin Rofeh, Zakary Singer, Bryan Toth, ElizabethWong, and M Ding. Dynamical parameter and state estimation in neuron models.
The dynamic brain: anexploration of neuronal variability and its functional significance , 2011.[27] Jingxin Ye, Paul J Rozdeba, Uriel I Morone, Arij Daou, and Henry DI Abarbanel. Estimating the biophysicalproperties of neurons with intracellular calcium dynamics.
Physical Review E , 89(6):062714, 2014.[28] Daniel Rey, Michael Eldridge, Mark Kostuk, Henry DI Abarbanel, Jan Schumann-Bischoff, and Ulrich Parlitz.Accurate state and parameter estimation in nonlinear systems with sparse observations.
Physics Letters A ,378(11-12):869–873, 2014.[29] J Ye, N Kadakia, PJ Rozdeba, HDI Abarbanel, and JC Quinn. Improved variational methods in statistical dataassimilation.
Nonlinear Processes in Geophysics , 22(2):205–213, 2015.[30] Nirag Kadakia, Eve Armstrong, Daniel Breen, Uriel Morone, Arij Daou, Daniel Margoliash, and Henry DIAbarbanel. Nonlinear statistical data assimilation for
HVC RA neurons in the avian song system. BiologicalCybernetics , 110(6):417–434, 2016.[31] Jun Wang, Daniel Breen, Abraham Akinin, Henry DI Abarbanel, and Gert Cauwenberghs. Data assimilation ofmembrane dynamics and channel kinetics with a neuromorphic integrated circuit. In
Biomedical Circuits andSystems Conference (BioCAS), 2016 IEEE , pages 584–587. IEEE, 2016.[32] Eve Armstrong, Amol V Patwardhan, Lucas Johns, Chad T Kishimoto, Henry DI Abarbanel, and George MFuller. An optimization-based approach to calculating neutrino flavor evolution.
Physical Review D , 96(8):083008,2017.[33] Jingxin Ye, Daniel Rey, Nirag Kadakia, Michael Eldridge, Uriel I Morone, Paul Rozdeba, Henry DI Abarbanel,and John C Quinn. Systematic variational method for statistical nonlinear state and parameter estimation.
Physical Review E , 92(5):052901, 2015.[34] First reported confirmation of coronavirus in New York City. . Accessed: 2020-05-19.[35] RN Thompson, JE Stockwin, RD van Gaalen, JA Polonsky, ZN Kamvar, PA Demarsh, E Dahlqwist, S Li, EveMiguel, T Jombart, et al. Improved inference of time-varying reproduction numbers during infectious diseaseoutbreaks.
Epidemics , 29:100356, 2019. 1536] Anne Cori, Neil M Ferguson, Christophe Fraser, and Simon Cauchemez. A new framework and software toestimate time-varying reproduction numbers during epidemics.
American journal of epidemiology , 178(9):1505–1512, 2013.[37] Jacco Wallinga and Peter Teunis. Different epidemic curves for severe acute respiratory syndrome reveal similarimpacts of control measures.
American Journal of epidemiology , 160(6):509–516, 2004.[38] PAUSE order in New York City takes effect 2020 March 22. . Accessed: 2020-05-19.[39] Daniel P. Oran and Eric J. Topol. Getting a handle on asymptomatic SARS-CoV-2 infection. . Ac-cessed: 2020-05-24.[40] Henrik Salje, Cécile Tran Kiem, Noémie Lefrancq, Noémie Courtejoie, Paolo Bosetti, Juliette Paireau, AlessioAndronico, Nathanaël Hoze, Jehanne Richet, Claire-Lise Dubost, et al. Estimating the burden of sars-cov-2 infrance.
Science , 2020.[41] Ruiyun Li, Sen Pei, Bin Chen, Yimeng Song, Tao Zhang, Wan Yang, and Jeffrey Shaman. Substantial undocu-mented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2).
Science , 368(6490):489–493, May 2020.[42] Qin Jing, Chong You, Qiushi Lin, Taojun Hu, Shicheng Yu, and Xiao-Hua Zhou. Estimation of incubation perioddistribution of COVID-19 using disease onset forward time: a novel cross-sectional and forward follow-upstudy. medRxiv , page 2020.03.06.20032417, March 2020. Publisher: Cold Spring Harbor Laboratory Press.[43] Roman WÃ˝ulfel, Victor M. Corman, Wolfgang Guggemos, Michael Seilmaier, Sabine Zange, Marcel A. MÃijller,Daniela Niemeyer, Terry C. Jones, Patrick Vollmar, Camilla Rothe, Michael Hoelscher, Tobias Bleicker, SebastianBrÃijnink, Julia Schneider, Rosina Ehmann, Katrin Zwirglmaier, Christian Drosten, and Clemens Wendtner.Virological assessment of hospitalized patients with COVID-2019.
Nature , 581(7809):465–469, May 2020.[44] Joseph A Lewnard, Vincent X Liu, Michael L Jackson, Mark A Schmidt, Britta L Jewell, Jean P Flores, Chris Jentz,Graham R Northrup, Ayesha Mahmud, Arthur L Reingold, et al. Incidence, clinical outcomes, and transmissiondynamics of hospitalized 2019 coronavirus disease among 9,596,321 individuals residing in california andwashington, united states: a prospective cohort study. medRxiv , 2020.[45] Dawei Wang, Bo Hu, Chang Hu, Fangfang Zhu, Xing Liu, Jing Zhang, Binbin Wang, Hui Xiang, ZhenshunCheng, Yong Xiong, Yan Zhao, Yirong Li, Xinghuan Wang, and Zhiyong Peng. Clinical Characteristics of138 Hospitalized Patients With 2019 Novel Coronavirusâ ˘A ¸SInfected Pneumonia in Wuhan, China.
JAMA ,323(11):1061–1069, March 2020. Publisher: American Medical Association.[46] Chaolin Huang, Yeming Wang, Xingwang Li, Lili Ren, Jianping Zhao, Yi Hu, Li Zhang, Guohui Fan, JiuyangXu, Xiaoying Gu, et al. Clinical features of patients infected with 2019 novel coronavirus in wuhan, china.
Thelancet , 395(10223):497–506, 2020.[47] Qifang Bi, Yongsheng Wu, Shujiang Mei, Chenfei Ye, Xuan Zou, Zhen Zhang, Xiaojian Liu, Lan Wei, Shaun ATruelove, Tong Zhang, et al. Epidemiology and transmission of covid-19 in shenzhen china: Analysis of 391cases and 1,286 of their close contacts.
MedRxiv , 2020.[48] Xiaobo Yang, Yuan Yu, Jiqian Xu, Huaqing Shu, Hong Liu, Yongran Wu, Lu Zhang, Zhui Yu, Minghao Fang,Ting Yu, et al. Clinical course and outcomes of critically ill patients with sars-cov-2 pneumonia in wuhan, china:a single-centered, retrospective, observational study.
The Lancet Respiratory Medicine , 2020.[49] Andreas Wächter. Short tutorial: getting started with ipopt in 90 minutes. In