key: cord-0565554-uo1g449a
authors: Frankenburg, Ian; Banerjee, Sudipto
title: A Compartment Model of Human Mobility and Early Covid-19 Dynamics in NYC
date: 2021-02-03
journal: nan
DOI: nan
sha: 3b2acd512ef2a7210c10bcbc4fff6b4adf4680b8
doc_id: 565554
cord_uid: uo1g449a

In this paper, we build a mechanistic system to understand the relation between a reduction in human mobility and Covid-19 spread dynamics within New York City. To this end, we propose a multivariate compartmental model that jointly models smartphone mobility data and case counts during the first 90 days of the epidemic. Parameter calibration is achieved through the formulation of a general Bayesian hierarchical model to provide uncertainty quantification of resulting estimates. The open-source probabilistic programming language Stan is used for the requisite computation. Through sensitivity analysis and out-of-sample forecasting, we find our simple and interpretable model provides evidence that reductions in human mobility altered case dynamics.

The global Covid-19 pandemic has underscored the importance of mathematical and statistical models in understanding disease dynamics, assessing policy efficacy, and examining counterfactual scenarios to formulate thorough cost-benefit analyses. Lockdown measures can have drastic impact on individual well-being as well as society and the economy at large (Bonaccorsi et al., 2020) , (Cutler and Summers, 2020) . Therefore a retrospective study of Covid-19 lockdown and mitigation measures can help policymakers and public health officials understand to what end such efforts were effective. The formulation of a mechanistic compartmental model is a pathway towards such goals. In this article, we review compartmental model methodology, construct our new Bayesian hierarchical model, and discuss numerical methods relevant for implementation and fitting to real-world data. The new compartmental model is a simple modification of the classical susceptible-infectious-removed (SIR) model and enables a mechanistic correspondence between smartphone mobility data and infection dynamics. This can provide evidence of how reduced mobility due to early lockdowns or mitigation measures within New York City influenced Covid-19 spread dynamics. Case count data is obtained from the official website of the City of New York, available at https://www1.nyc.gov/site/doh/covid/covid-19-data. page. Population transit mobility data is obtained from https://covid19.apple.com/mobility and consists of anonymized Apple iPhone transit usage reported as a percent relative to baseline. Starting from the day after Governor Andrew Cuomo declared a state of emergency in New York State on March 7th, 2020, both the raw case count and arXiv:2102.01821v2 [stat.AP] 22 Jun 2021 transit mobility time series are presented below. Our end goal is then to establish a relationship between the two series shown in Figure 1 . 

Mathematical modeling in epidemiology has a long history, famously dating back to the eighteenth century with the work of Bernoulli (Dietz and Heesterbeek, 2002) or the mid-nineteenth century through John Snow's modeling of the cholera outbreak in London (Shiode et al., 2015) . However, at the turn of the early twentieth century, mathematical epidemiology turned to the modern theory of dynamical systems analysis to understand outbreak evolution. In this section, we review the popular SIR model and subsequent mathematical analysis used to glean both qualitative and quantitative understanding of the dynamical system. This simple framework provides the necessary foundation for more complicated compartmental models with more population states. For examples of other compartment models designed to study early Covid-19 outbreak dynamics, see (Hao et al., 2020) or (Wang et al., 2020a) .

In modeling population-level data with a compartmental system, the population is typically subdivided into separate homogeneous groups. Here we focus on reviewing the simple SIR model developed in 1927 by A. G. McKendrick and W. O. Kermack to model a plague outbreak in Bombay (Kermack and McKendrick, 1927) . In such a model, the population is divided into susceptible (S), infectious (I), and removed (R) groups. Individuals then progress through the various states at certain rates over time. The mathematical description of this changing system is the coupled set of ordinary differential equations

(1)

The progression of the disease throughout the population depends upon the contact rate between susceptible and infectious individuals, the probability of transmission upon contact, and the prevalence of disease. To mathematically capture these factors and express the rate of new infections, let λ be defined as a per capita contact rate among individuals. In this way, λS(t) will give the average number of susceptible contacts over time. Now let p be the probability that a contact results in a new infection. Finally, the prevalence of the disease at time t is by definition I(t)/N . Combining these terms gives pλS(t)I(t)/N as the incidence rate. In defining the effective contact rate β as the product of the per capita contact rate λ and transmission probability p, the necessary form in equation (1) is recovered.

The remaining parameter γ is interpreted by considering that 1/γ is the average sojourn time of an individuals within compartment I. It should also be noted the population is fixed throughout time since N = S + I + R. This can alternatively be seen since adding the terms in (1) gives zero. With meaning associated to the compartmental parameters, we next turn to an overview of the mathematical analysis involved in analyzing the basic SIR dynamical system.

The system of differential equations in (1) are nonlinear, arising from the term S(t)I(t). Linear stability analysis is the workhorse to understand the behavior of nonlinear dynamical systems and has a fundamental connection to the basic reproductive number R 0 , popularized recently through media coverage of the Covid-19 pandemic. R 0 is defined roughly to be the expected number of subsequent infections resulting from a single infected individual. In this section, we briefly review linear stability analysis and make the connection to R 0 . In the next section, we highlight the computation of R 0 for a general class of compartmental models.

A steady state or equilibrium of a dynamical system is a point x * where the system of differential equations evaluated at x * is zero for all t. In this way, compartmental contents within the system are not changing over time. In the SIR model of (1), an important steady state is the so-called disease-free equilibrium of {(S * , 0, 0) : S * ≥ 0}. A natural subsequent question is the behavior of the system around small perturbations of the equilibrium. In computing the Jacobian about such a point, we linearize and are afforded tractable analysis. A heuristic justification arises from considering a Taylor expansion of the system about the disease-free equilibrium and ignoring high-order terms since the perturbation is assumed small. After dropping the explicit dependence on time t from equation (1) to avoid clutter, the Jacobian of the system is computed as

Through elementary matrix operations, J can be transformed to a triangular matrix. The eigenvalues are then found from inspection. Evaluating the Jacobian J at (S * , 0, 0) to linearize about the steady state results in eigenvalues of λ 1 = 0 and λ 2 = βS * /N − γ. The sign of these eigenvalues then determine the stability of the equilibrium point. The eigenvalue of λ 1 is ignored, as it corresponds to a line of equilibrium values S * . In this way, the second eigenvalue of λ 2 is of main interest. If N/S * < β/γ, then λ 2 > 0 and the steady state is unstable; otherwise it is stable. In epidemic terms, an outbreak occurs if the disease-free equilibrium is unstable. This ratio β/γ acts as a bifurcation parameter in determining if an outbreak will occur and is thus afforded the fancy title of basic reproductive number. Letting the equilibrium point S * be the population size so that S * = N , a simple relation emerges: if R 0 > 1 the disease continues to spread but dies out otherwise. The effective reproductive number R t then extends R 0 by accounting for a changing susceptible population over time and is defined as R t := R 0 S(t)/N . A general method to compute R 0 for more elaborate compartmental models will be discussed in the next subsection.

For general compartment models that extend the simplistic SIR framework, computing the basic reproductive number can be difficult. (Diekmann et al., 2009) and (Heffernan et al., 2005) describe a general method to compute R 0 called the Next Generation Matrix or Spectral Radius Method, which we briefly review and compute for the SIR model.

Let dX(t)/dt represent a general coupled system of differential equations describing a compartmental model with n components and m infectious states. Define a vector-valued F(X(t)) to be a function where each component specifies flow rate into one of the respective m infected compartments. Similarly, define a function V(X(t)) where each component specifies the flow rate out of a respective infectious compartment.

Next, F and V are linearized about the disease-free equilibrium point by computing the Jacobian. (Diekmann et al., 2009 ) prove the Jacobian with respect to each infectious state will take the form

where A and B are m × m matrices. The next generation matrix is then defined as AB −1 . The basic reproductive value of R 0 is subsequently the spectral radius or largest eigenvalue of the next generation matrix. In the case of the SIR model, F(X(t)) := −βSI/N , while V(X(t)) := γI. It follows that the necessary Jacobians evaluated at the disease-free equilibrium of (N, 0, 0) results in AB −1 = β/γ and agree with the previous section.

In this section, we detail the construction of our compartmental model designed to formulate an understanding of how reduction in human mobility might have altered early infection dynamics within New York City. The model is parsimonious, in that it consists of only four compartments and shares many well-established mathematical properties of the SIR model discussed in the background section. The dynamical system proposed comprises a closed population divided into susceptible, lockdown (L), infectious, and removed states. Since the time horizon under investigation is short, we choose to ignore demographic factors such as birth, death, and migration. Below in Figure 2 is a visualization of the population progression through the different states. The system of differential equations describing the changing system are

with initial conditions of S(0) = N − i 0 , L(0) = 0, I(0) = i 0 , and R(0) = 0. An important qualitative feature of our model is that susceptible individuals are temporarily moved into the lockdown compartment to reflect social distancing, mitigation measures, and reduced mobility. Over time, individuals are reintroduced into the susceptible population out of the lockdown state. Through an application of the next generation matrix method described in section 2, R 0 can be seen to be equivalent to the standard SIR model, i.e. R 0 = β/γ.

We present our methodology in a general hierarchical framework to facilitate Bayesian inference of compartmental system parameters in equation (4). In this hierarchical formulation, we can be explicit about the role of the mechanistic system, requisite numerical integration, and underlying process parameters. This section will construct the statistical model piece-by-piece. The hierarchy of connected components in the model can be visualised bottom-up as follows,

We first establish notation to represent the mechanistic system of differential equations in the middle of the hierarchy. Let the equations in (4) be denoted by F, where d dt

and we have suppressed the dependence of {R 0 , γ, a, b} in F(X(t)). The system is reparameterized in terms of R 0 rather than β to be more epidemiologically interpretable and results from a simple transformation of β = γR 0 . To recover the system states X(t), the system of differential equations must be solved. Given fixed values of R 0 , γ, a, and b, the solution to the system of differential equations is a vector-valued function

This solution will be necessary to connect with the observed data.

The top of the hierarchy is described by formulating a measurement process for the two outcome variables. Let the first outcome of interest be labeled Y L (t) and represent the percent of the population removed from the susceptible compartment by adhering to mitigation protocol. The subscript L is used for a reminder that this data is used to gain information on the lockdown compartment. Likewise, the second outcome is labeled Y I (t) and denotes the observed case counts over time. To model observation error in Y L (t), we choose a Beta distribution dependent upon a parameter φ 1 to control dispersion, i.e.,

Notice L(t) is necessarily scaled by the population size N to respect the support of the beta distribution. As we seek to inform the L compartment through cell phone mobility data, the beta distribution is a natural choice because the data will be anonymized and reported as a percentage of nominal movement. This parameterization of the beta distribution has expectation and variance

To model observation noise in Y I (t), we use a negative binomial to account for overdispersion,

Stan provides an alternative parameterization of the negative binomial called neg_binomial_2 with first two moments of

In this way, φ 2 is viewed as a dispersion parameter. The full hierarchical description of the model can be completed by introducing prior distributions on the system parameters governing the differential equations. Writing the model in full,

The prior distributions on system parameters are weakly-informative. However, the prior distribution on a might at first appear suspect. Through prior predictive checks, we find that placing a uniform prior on a results in the SLIR model a priori favoring no epidemic breakout, as susceptibilities are removed from the population too quickly. Since the classical SIR model is a special case of our SLIR model as a → 0, we place mass closer to 0 through the Beta(1,5) distribution to ensure the model generates reasonable predictions before seeing the data. The hierarchical model of the previous section crucially depends upon the numerical solution to a coupled set of differential equations of (4). In the appendix section, we detail the internal workings of Stan's numerical optimization routines. Finally, efficient Bayesian analysis of parameters within the set of nonlinear differential equations relies upon the efficiencies gained through Hamiltonian Monte Carlo (HMC). A detailed review of this methodology is also included in the appendix section.

In this section, we present two simulation studies and conclude with the New York City analysis. First, the proposed SLIR compartmental model is used to simulate data from two lockdown scenarios that affect human mobility differently. We then fit our Bayesian model to assess whether the true parameter values are adequately recovered. After, we analyze the real-world mobility and case count data that initially inspired the model formulation.

To illustrate the nonlinear dynamics of which our model can capture, the first simulation reflects the idealized scenario of strict adherence to lockdown and mitigation measures, when population movement is quickly reduced in the early stages of the outbreak and remains reduced for the next 90 days. In this case, individuals move from the S compartment to the L compartment quickly and are slowly reintroduced into the susceptible population so that population movement decreases to about 60% relative to baseline within the first 20 days.

In the second scenario, we consider weak adherence to mitigation measures and illustrate the substantial change in dynamics by only altering the speed of flow back into the susceptible population from the L compartment. In this case, the peak percentage of the population in the lockdown compartment is 30% but quickly diminishes. In both simulations, the population size is fixed to N = 10, 000, i 0 = 1, and the true data-generating process is shown below as a dashed red line. The median of the posterior predictive distribution and 95% credible intervals are shown in blue. We fit our model to both scenarios using Stan's NUTS algorithm with 4 chains and 5,000 iterations each, the first half of which are discarded as warm-up. The convergence of the parameter chains are judged by inspecting the trace plots along with the Gelman-RubinR values, which compares the variation between chains to the variation within (Gelman and Rubin, 1992) . Ideally, theR value is close to one. These simulation results are displayed in in Table 1 . In both cases, the Bayesian hierarchical model is able to infer the structural parameters of the SLIR model. Notice the change in case counts of both scenarios resulting from different susceptible population sizes.

The entire lead-up thus far was requisite background material for compartmental model inference and application to realworld data. As mentioned in the introduction, our main motivation for this article was to understand how a reduction in mobility affected early Covid-19 dynamics specifically within NYC. For convenience, we restate the data sources. Case counts are reported by the official website of the City of New York, available at https://www1.nyc.gov/site/doh/ covid/covid-19-data.page and mobility data is hosted at https://covid19.apple.com/mobility. Since the Apple mobility data reflects a percent decrease in movement, it must first be transformed by subtraction from unity to adhere with the SLIR compartmental model structure. In other words, to prepare the mobility data for use in the hierarchical model, it must first be subtracted from one so that it no longer represents a percent decrease in transit mobility but rather a percent increase in individuals adhering to mitigation measures. Finally, we mention again that although the first case of Covid-19 in New York City was recorded on February 29th, we align our movement and case data to begin on March 8th, 2020, the day after Governor Andrew Cuomo declared a state of emergency in New York State. Finally, we take the initial number of cases to be the cases recorded on March 8th, 2020, and the population is fixed at 8,336,817 as determined by the US Census (U.S. Census Bureau, 2019).

To achieve parameter calibration, we fit in Stan our full hierarchical model using four chains run for 10,000 iterations each. We discard the first 5,000 as warm-up. Sufficient posterior exploration is assessed by examining parameter chain plots below and assessingR values, shown in The fitted time series are presented below on the right, along with the raw data used to train the model. On the left, prior predictive distributions are included to illustrate the degree in which Bayesian learning occurs after observing the data. We also include the fit of an SIR model to illustrate its inability to capture the dynamics. To assess the structural fit of our hypothesized SLIR mechanism, we next interpret parameter values to ensure they are logical and consistent with outside literature. The R 0 estimate is cross-referenced with those of other popular online models. Using only death statistics as reported by Johns Hopkins University, (Gu, 2020 ) embeds a SEIR model in a machine learning framework for many regions across the United States and 70 countries. In this work, R 0 is estimated for NYC to be between 5.0 and 5.8, in close agreement with our model. As an additional point of reference, (Ives and Bozzuto, 2020) provide an alternative methodology that fits a time series state-space model to death counts and explicitly accounts for reporting delays. R 0 is estimated to be 6.3 with a 95% confidence interval of 4.5-9. Our model thus has the added benefit of mechanistic interpretability as well as a tighter credible interval.

The posterior of γ has a median value of 0.21 and a 95% credible interval of (0.191, 0.234). From this, we arrive at an estimate of approximately 5 days for the average infectious removal time. This estimate could be reasonable assuming asymptomatic or presymptomatic transmission is possible and that individuals isolate at the onset of symptoms; see (Gandhi et al., 2020) for evidence. Outside work estimates symptom onset time to also be around 5 days. For example, (Lauer et al., 2020) use 181 confirmed Covid-19 cases to estimate a median symptom onset time of 5.1 days with a 95% confidence interval (4.5, 5.8) days. A systematic review by (Madjid et al., 2020 ) estimates a mean symptom onset time of 5.2 days with a 95% confidence interval of (4.1, 7.0) days. These estimates provide credence in establishing a correspondence between population movement reduction and infection dynamics with the mechanistic SLIR model.

We conclude the New York City analysis by assessing the sensitivity of the model to changing mobility levels.

Additionally, we assess out-of-sample predictive capacity of the SLIR model. To perform the sensitivity analysis, we first generate a range of hypothetical mobility scenarios, from full mobility reduction through extremely stringent lockdown measures to a more mild decline. This is displayed in Figure 5 , with the true, real-world observed mobility levels highlighted as blue. It is important to note the explosive case growth with mobility reduction levels under 60% due to the nonlinear dynamics present in SIR-type models. This is evidence that a mobility reduction significantly altered infections within the city, assuming SIR-type dynamics.

To conclude this section, we analyse the out-of-sample forecasting ability of the SLIR model when trained on only a subset of the ninety day period. We consider for illustration a two, three, and four week training intervals. These are indicated by vertical dashed lines in Figure 6 . The predictive median is illustrated with a solid line. With only two weeks observed, the predictive median is reasonably representative of the future trajectory but with very large uncertainty. The upper 95% predictive curve for the two week window reaches roughly 400,000 infections, but is clear that the predictive interval quickly contracts as more data is observed (cf. Figure 4) . Both the three week and four week training periods have contracted credible interval forecasts that contain the observed data. Finally, we mention that we were unable to fit the standard SIR model to the New York City data. Using the quasi-Newton optimization functions available in Stan, we found the SIR model fit was highly unstable across a range of starting values indicating extreme multi-modality in the likelihood surface.

5 Discussion and Future Work

In this work, we have formulated a basic extension of the classical SIR model to jointly fit cell phone transit mobility data and case counts. By jointly modeling two outcome variables, we establish a mechanistic correspondence between reduced mobility and infection dynamics. The applied analysis and findings, however, are limited to NYC during the first 90 days. The disease progression throughout the region was well-approximated by deterministic dynamics and modeling with a simple compartment system. In other geographic regions, the dynamics might not be as suitably well-behaved or understood. It is difficult to capture the myriad of factors contributing to disease spread throughout a population, and often a stochastic model may be more appropriate. Additionally, the time horizon considered in our applied analysis is relatively short. As a disease progresses throughout the population and becomes more widespread, the strategy of informing the susceptible population numbers through cell phone mobility data becomes more difficult. It is also worth consideration about whether in general smartphone mobility data can be considered a representative sample of individuals' adherence to mitigation protocols. It is therefore difficult to build a quantitative understanding of the dynamics between lockdown measures and disease progression over long periods of time. A direction for future research is to establish modeling strategies for such scenarios. This will introduce additional challenges of accounting for population immunity, vaccinations, and natural seasonality effects.

We briefly discuss for both discrete and continuous time how the deterministic SLIR model can be relaxed through a state-space (or Hidden Markov Model) framework, where stochasticity is introduced into the underlying driving epidemic process. We first discuss a discrete-time extension by modifying equation (9) so that the numerical solution to dX(t) dt is embedded within a likelihood function at the discrete observation time points. In this way, equation (9) is modified to include p(X(t)|θ, R 0 , γ, a, b), where θ is introduce to accommodate new parameters. In other words, stochasticity is introduced to the driving dynamics through the choice of suitable probability distribution that is a function of the numerical solution of the system of differential equations. This approach has been successfully applied to forecast seasonal influenza (Osthus et al., 2017) , (Dukic et al., 2012) , as well as to assess Covid-19 interventions in China (Wang et al., 2020b) and Japan (Kobayashi et al., 2020) .

An alternative approach is to incorporate continuous-time stochasticity directly by expressing dX(t) as a stochastic differential equation (SDEs) (Øksendal, 2003) . In this case, one modification of equation (9) 

where W(t) = (W 1 (t), W 2 (t), W 3 (t), W 4 (t)) is a vector of standard Wiener processes or Brownian motions, µ is some vector-valued function, and L is a compatible matrix. Depending on the structure of µ, there may or may not exist an explicit solution. A SDE that affords a closed form solution is the Ornstein-Uhlenbeck (OU) process. See (Chatzilena et al., 2019) for Stan code in the case where the infectious compartment of the SIR model is endowed an OU process to allow for random fluctuations. See also (Wang et al., 2018) for an analysis of the susceptible-infectious-susceptible model where time-varying parameters are introduced through an OU process. In both examples, the analytic nature of the solution to the SDE enables efficient implementation in Stan or alternative software. In cases where no closed-form solution to (10) exists, inference in such a system is closely related to Bayesian filtering and smoothing, amenable to such methods as the Kalman filter and extensions (Särkkä and Solin, 2019) . The New York City data was captured well by deterministic dynamics, but often elsewhere disease progression is not suitably captured within such a framework. Bayesian inference for SDEs modeling infectious disease dynamics is thus an attractive area of future research.

The hierarchical model of the previous section crucially depends upon the numerical solution to a coupled set of differential equations. We briefly review the popular numerical integration schemes and connect them with our hierarchical model of the compartmental system in (4). In practice, solutions to the SIR model are numerically approximated and entail specific computational challenges. Runge-Kutta (RK) is a classical numerical method (Struthers and Potter, 2019; Press et al., 2007) in which we employ to numerically solve this set of nonlinear differential equations. RK generalizes the well-known Euler method for iteratively solving systems of differential equations. In the following, we formulate the RK method in vector notation for the proposed SLIR compartmental model.

Consider X(t) in (5). From the fundamental theorem of Calculus, we know that

Given the initial value X(t 0 ) at time t 0 , numerically solving for X(t) amounts to constructing a sequence {X n : n = 0, 1, . . . , T }, where X n := X(t n ) and t n = t n−1 + ∆t are equispaced time points for n = 0, 1, . . . , T , by approximating the integral of F(X(t), t) in (11). Using this approximation, a sequence is generated starting from some t 0 as,

where the initial condition X 0 := X(t 0 ) is given. There are several customary choices for such approximations, but for most practical purposes one need not look beyond the following:

Euler's approximation (first equation in (13)) is the simplest of the approximations which yields Euler's Method. In this case, the sequence is constructed as

Starting with X 0 , each element of the sequence in (14) is computed since F(X n , t n ) is available. The spacing between the time points, ∆t, is specified by the user and controls the resolution of the numerical solution.

Euler's approximation is easy to execute and based upon a first order Taylor expansion. It is also the least accurate. The Trapezoidal and Modified Euler approximations in (13) are based upon numerical integration using trapezoidal areas and the midpoint approximation. Both these methods help improve Euler's method with the Modified Euler outperforming the trapezoidal rule in terms of accuracy. However, the Trapezoidal and the Modified Euler methods involve X n+1 and X n+1/2 := X(t n+1/2 ) with t n+1/2 := t n + (∆)t/2, respectively, which are unknown at iteration n. In fact, X n+1 is the very quantity we wish to compute at iteration n. Therefore, we substitute these unknown quantities with their first order (Euler) approximations in (14) that are available at iteration n. For each n we compute a n = F (X n , t n ), b n = F (X n + (∆t)a n , t n+1 ) and c n = F X n + (∆t/2)a n , t n+1/2 and transition from X n to X n+1 as X n+1 = X n + (∆t) (an+bn) 2 (Trapezoidal) ; X n + (∆t)c n (Modified Euler) for n = 0, 1, . . . , T − 1 ,

Both methods in (15) deliver noticeable improvements over Euler's method in (14). Numerical error from (15) are much smaller in magnitude and grow less quickly. This can be explained by observing that while (14) depends only upon the data available at t n (only one data point), the two methods in (15) use current data at t n along with estimates of the slope at a point that lies in the future. While these estimates are computed using only the currently available data, they still produce substantially improved estimates. Also see Treiber and Kanagaraj (2015) for comparisons among different numerical integration schemes.

Higher order Taylor expansions produce other iterative schemes. Thus, second order methods emerge from

A second order iterative scheme corresponding to (16) updates

where we have used the multivariable chain-rule of derivatives to evaluate the derivative of F(X(t), t) with respect to t, ∂F ∂X is the matrix with (i, j)-th element being the partial derivative of the i-th element of F(X(t), t) with respect to the j-th variable in X, and d dt X(t) t=tn = F(X n , t n ). Unfortunately, computing (17) requires the derivatives of F(X(t), t) and may, in general, be numerically cumbersome.

The Runge-Kutta methods are among the most conspicuous of numerical methods for solving systems of ordinary differential equations. The underlying idea is to achieve the same accuracy as Taylor series updates without requiring higher order derivatives of F(X(t)). We can motivate this approach from the earlier methods. In (15), the Trapezoidal and Modified Euler methods define the updates using a n , b n and c n that are completely specified. In particular, observe that the Trapezoidal method updates using a weighted average of a n and b n . Instead of prescribing a n and b n , the Runge-Kutta approach prefers to find weighted averages to ensure that the approximation matches that from a Taylor series expansion such as in (16) or (17). Therefore, a second-order Runge-Kutta method (RK2) writes X n+1 = X n + (∆t) {ω 1 a n + ω 2 F(X n + (∆t)βa n , t n + (∆t)α)}

and seeks to find ω 1 , ω 2 , α and β so that the approximation matches (17). Substituting the first-order expansion, F(X n + (∆t)βa n , t n + α(∆t)) = F(X n , t n ) + (∆t)β ∂F ∂X X=Xn a n + (∆t)α ∂F ∂t t=tn + O((∆t) 2 ) = a n + (∆t)β ∂F ∂X X=Xn a n + (∆t)α ∂F ∂t t=tn

into the right hand side of (18) and comparing with the expansion (17) we find that the two expansions are equivalent if

RK2 specifies ω 1 = ω 2 = 1/2 and α = β = 1, which, when substituted into (18), yields the Trapezoidal approximation.

More generally, the explicit RK methods of order s specify updating schemes

where k 1 = F(X n , t n + (∆t)α 1 ), and k i = F X n + (∆t) 

In particular, the RK4 method specifies the following values in (21):

The appropriate step size ∆t is hard to determine. An advantage of Stan's implementation is an adaptive step size is used by comparing the solution obtained using the four term approximation of above as well as a five term approximation. If these approximations agree, the algorithm proceeds, otherwise a new step size is calculated. More details can be found in the user manual (Stan Development Team, 2020), but the result is a fast, efficient procedure of high accuracy.

Our choice of implementation in Stan, as opposed to more traditional BUGS or JAGS (Plummer, 2003) , is pragmatic. First, the latter languages are declarative and built upon graphical models. In contrast, Stan is a fully imperative programming language. Additionally, built-in differential equation routines are included such as the Runge-Kutta numerical solver described in the previous section. This makes the software implementation of our model more natural and readable. More importantly, however, the parameters in a nonlinear compartmental model are often highly correlated, as demonstrated in below in Figure 7 . Hamiltonian Monte Carlo (HMC) is more equipped to sample from complex posterior distributions with high autocorrelations than standard Metropolis schemes. The most popular presentation of Hamiltonian Monte Carlo is by way of analogy with statistical mechanics. Let θ be an arbitrary d-dimensional parameter vector. To sample efficiently from the posterior of θ after conditioning on data, an idealized physical system is introduced to leverage the geometry of the underlying manifold on which θ lives. We will not embark upon a comprehensive development of HMC here, referring the reader to excellent introductory expository articles by (Neal, 2012) and (Betancourt, 2018) . Instead, we provide a heuristic account of the HMC algorithm and why it works.

We begin by recalling the more conspicuous Metropolis random walk and the concept of detailed balance. Let p(θ | Y) be the posterior distribution from which we wish to sample. As it may be difficult to directly sample from p(θ | Y), the Metropolis random-walk algorithm constructs a Markov chain with p(θ | Y) as its stationary distribution. Given an initial value θ (0) , at iteration t we draw a "proposed" value θ * from a symmetric distribution q(· | θ (t−1) ). A simple yet effective choice for many applications is q(· | θ (t−1) ) = N (· | θ (t−1) , V), where V is a fixed variance covariance matrix that helps tune the algorithm, although in general the proposal can be generated from any symmetric distribution. After generating θ * we simulate a coin with probability of heads min 1,

and set θ (t) = θ * if it is a head.

Otherwise we set θ (t) = θ (t−1) . That p(θ | Y) is indeed the desired stationary distribution can be seen as follows. Let us assume that the current state θ (t−1) = a is a draw from p(θ | Y) and consider the possibility of moving to θ (t) = b. This conditional probability is given by the transition probability that a value of b is proposed from q(· | a) and that this value is accepted. Hence,

The form of the transition probability in (24) implies time-reversibility (or detailed balance) in the following sense:

where we have used the symmetry q(a | b) = q(b | a) in the fourth equality in (25). It follows that the draw of θ (t) is also from p(θ | Y) because

The underlying idea behind HMC is that instead of generating the proposed value from a random distribution, we use a deterministic symplectic integrator to propose θ * . This symplectic integrator is designed based upon Hamiltonian dynamics. Suppose that we wish to sample from p(θ | Y), where θ ∈ R d . We introduce an auxiliary variable r ∈ R d so that we can efficiently sample from the joint density p(θ, r | Y). If (θ (t) , r (t) ) ∼ p(θ, r | Y), then P (θ (t) = b) = P (θ (t) = b, r (t) = u)du = p(b, u | Y)du = p(b | Y) .

Hence, sampling from the joint density p(θ, r | Y) results in samples from p(θ | Y).

The auxiliary variable, r, is also called the "momentum" in Hamilton dynamics. For our purposes, it suffices to specify that p(θ, r | Y) = p(θ | Y) × p(r). Therefore, p(r | θ, Y) = p(r) which means that r is independent of the data Y and the model parameters θ. More specifically, we assume that p(r) = N (r | 0, I d ) ∝ exp − 1 2 r r . Therefore, log p(θ, r | Y) = constant + log p(θ | Y) − 1 2 r r .

The above density can be looked upon as a physical system subject to Hamiltonian dynamics, where θ is a particle's position in R d and r is the particle's momentum.

In order to sample from (27), a simple HMC algorithm proceeds closely on the lines of the Metropolis random walk described earlier, but replaces the random generation of a proposed value for θ by a symplectic integrator constructed from Hamiltonian dynamics. With the current state (θ (t−1) , r (t−1) ), we begin iteration t by drawing the momentum variable r * ∼ N (0, I d ). Setting r (t) = r * we perform L steps of a symplectic integrator (also known as "leapfrog"), where each step comprises the following:

r (t+ /2) = r (t) + 2 ∇ θ L(θ); θ (t−1+ ) = θ (t−1) + r (t+ /2) ; r (t+ ) = r (t+ /2) + 2 ∇ θ L(θ) ,

where L(θ) = log p(θ | Y). Letθ andr be the output of (28) at the end of L steps. The values ofθ andr are considered the "proposed" values at iteration t and accepted as (θ (t) , r (t) ) = (θ, −r) with acceptance probability

. This last Metropolis step together with the negation of the momentum variable in the final update ensures time-reversibility as in (25) and, as seen in (26), maintains p(θ | Y) as the stationary distribution.

We provide some further intuition on the time-reversibility of the simple HMC algorithm. The key to this result is that the leapfrog iteration in (28) preserves volumes. To be slightly more precise, let D be a small region in the (θ, r) space and suppose the L leapfrog steps maps D to a regionD. Then D andD both have the same volume. We write the transition probability from D toD as 

where δV represents the volume of D andD, and D p(θ, r)dθdr = exp (−H(D)). If (θ (t−1) , r * is drawn from the joint density in (27) = P ((θ (t−1) , r * ) ∈D, (θ (t) , r (t) ) ∈ D) ,

where the last equality follows from the symmetry in the expression above it.

This procedure, while maintaining the stationary distribution through time-reversibility, introduces a host of complexities. Perhaps most importantly, tuning the many parameters needed in this process is inherently difficult. This motivated the development of an automatic procedure known as the No-U-Turn-Sampler (NUTS) (Hoffman and Gelman, 2011) . This algorithm achieves significant efficiency over the simple HMC algorithm described above by either explicitly avoiding a U-turn to previously explored region or terminating after a pre-defined number of exploration steps. In this way, the algorithm is guaranteed to only explore new areas of the space. This efficient exploration results in typically faster convergence and higher effective sample sizes per iteration as compared to classical MCMC.

A conceptual introduction to hamiltonian monte carlo

Economic and social consequences of human mobility restrictions under covid-19

Contemporary statistical inference for infectious disease models using stan

The covid-19 pandemic and the $16 trillion virus

The construction of next-generation matrices for compartmental epidemic models

Daniel bernoulli's epidemiological model revisited

Tracking epidemics with google flu trends data and a state-space seir model

Asymptomatic transmission, the achilles' heel of current strategies to control covid-19

Inference from Iterative Simulation Using Multiple Sequences

Covid-19 projections using machine learning

Reconstruction of the full transmission dynamics of covid-19 in wuhan

Perspectives on the basic reproductive ratio

The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo

State-by-state estimates of r0 at the start of covid-19 outbreaks in the usa. medRxiv

A contribution to the mathematical theory of epidemics

Predicting intervention effect for covid-19 in japan: state space modeling approach

The incubation period of coronavirus disease 2019 (covid-19) from publicly reported confirmed cases: Estimation and application

Potential effects of coronaviruses on the cardiovascular system: A review

Mcmc using hamiltonian dynamics. Handbook of Markov Chain Monte Carlo

Stochastic Differential Equations: An Introduction with Applications. Hochschultext / Universitext

Forecasting seasonal influenza with a state-space sir model

Jags: A program for analysis of bayesian graphical models using gibbs sampling

Numerical Recipes 3rd Edition: The Art of Scientific Computing

The mortality rates and the space-time patterns of john snow's cholera epidemic map

Stan Modeling Language Users Guide and Reference Manual

Differential Equations: For Scientists and Engineers

Applied Stochastic Differential Equations

Comparing numerical integration schemes for time-continuous car-following models

Evolving epidemiology and impact of non-pharmaceutical interventions on the outbreak of coronavirus disease 2019 in wuhan, china

An epidemiological forecast model and software assessing interventions on covid-19 epidemic in china. medRxiv

A stochastic differential equation sis epidemic model incorporating ornstein-uhlenbeck process