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 [15]. 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

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, $\lambda 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\lambda 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 (2.1) is recovered.

The remaining parameter γ is interpreted by considering that $1/\gamma $ 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 (2.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 (2.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 ${\mathcal{R}_{0}}$, popularized recently through media coverage of the Covid-19 pandemic. ${\mathcal{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 ${\mathcal{R}_{0}}$. In the next section, we highlight the computation of ${\mathcal{R}_{0}}$ for a general class of compartmental models.

A steady state or equilibrium of a dynamical system is a point ${\mathbf{x}^{\ast }}$ where the system of differential equations evaluated at ${\mathbf{x}^{\ast }}$ is zero for all t. In this way, compartmental contents within the system are not changing over time. In the SIR model of (2.1), an important steady state is the so-called disease-free equilibrium of $\{({S^{\ast }},0,0):{S^{\ast }}\ge 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 (2.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^{\ast }},0,0)$ to linearize about the steady state results in eigenvalues of ${\lambda _{1}}=0$ and ${\lambda _{2}}=\beta {S^{\ast }}/N-\gamma $. The sign of these eigenvalues then determine the stability of the equilibrium point. The eigenvalue of ${\lambda _{1}}$ is ignored, as it corresponds to a line of equilibrium values ${S^{\ast }}$. In this way, the second eigenvalue of ${\lambda _{2}}$ is of main interest. If $N/{S^{\ast }}<\beta /\gamma $, then ${\lambda _{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 $\beta /\gamma $ 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^{\ast }}$ be the population size so that ${S^{\ast }}=N$, a simple relation emerges: if ${\mathcal{R}_{0}}>1$ the disease continues to spread but dies out otherwise. The effective reproductive number ${\mathcal{R}_{t}}$ then extends ${\mathcal{R}_{0}}$ by accounting for a changing susceptible population over time and is defined as ${\mathcal{R}_{t}}:={\mathcal{R}_{0}}S(t)/N$. A general method to compute ${\mathcal{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. [5] and [12] describe a general method to compute ${\mathcal{R}_{0}}$ called the Next Generation Matrix or Spectral Radius Method, which we briefly review and compute for the SIR model.

Let $d\mathbf{X}(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 $\mathbf{F}(\mathbf{X}(t))$ to be a function where each component specifies flow rate into one of the respective m infected compartments. Similarly, define a function $\mathbf{V}(\mathbf{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. [5] prove the Jacobian with respect to each infectious state will take the form where A and B are $m\times m$ matrices. The next generation matrix is then defined as $\mathbf{A}{\mathbf{B}^{-1}}$. The basic reproductive value of ${\mathcal{R}_{0}}$ is subsequently the spectral radius or largest eigenvalue of the next generation matrix. In the case of the SIR model, $\mathbf{F}(\mathbf{X}(t)):=-\beta SI/N$, while $\mathbf{V}(\mathbf{X}(t)):=\gamma I$. It follows that the necessary Jacobians evaluated at the disease-free equilibrium of $(N,0,0)$ results in $\mathbf{A}{\mathbf{B}^{-1}}=\beta /\gamma $ and agree with the previous section.

We present our methodology in a general hierarchical framework to facilitate Bayesian inference of compartmental system parameters in equation (3.1). 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 visualized 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 (3.1) be denoted by F, where and and we have suppressed the dependence of $\{{\mathcal{R}_{0}},\gamma ,a,b\}$ in $\mathbf{F}(\mathbf{X}(t))$. The system is reparameterized in terms of ${\mathcal{R}_{0}}$ rather than β to be more epidemiologically interpretable and results from a simple transformation of $\beta =\gamma {\mathcal{R}_{0}}$. To recover the system states $\mathbf{X}(t)$, the system of differential equations must be solved. Given fixed values of ${\mathcal{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 ${\phi _{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, ${\phi _{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\to 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 (3.1). 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.

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 in Figure 3.

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-Rubin $\hat{R}$ values, which compares the variation between chains to the variation within [9]. Ideally, the $\hat{R}$ value is close to one. These simulation results are displayed 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 real-world 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 [29].

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 assessing $\hat{R}$ values, shown in Table 2.

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 ${\mathcal{R}_{0}}$ estimate is cross-referenced with those of other popular online models. Using only death statistics as reported by Johns Hopkins University, [10] embeds a SEIR model in a machine learning framework for many regions across the United States and 70 countries. In this work, ${\mathcal{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, [14] provide an alternative methodology that fits a time series state-space model to death counts and explicitly accounts for reporting delays. ${\mathcal{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 [8] for evidence. Outside work estimates symptom onset time to also be around 5 days. For example, [17] 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 [18] 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 4, 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 5.

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.

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 (3.7) so that the numerical solution to $\frac{d\mathbf{X}(t)}{dt}$ is embedded within a likelihood function at the discrete observation time points. In this way, equation (3.7) is modified to include $p(\mathbf{X}(t)|\boldsymbol{\theta },{\mathcal{R}_{0}},\gamma ,a,b)$, where $\boldsymbol{\theta }$ 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 [21], [7], as well as to assess Covid-19 interventions in China [31] and Japan [16].

An alternative approach is to incorporate continuous-time stochasticity directly by expressing $d\mathbf{X}(t)$ as a stochastic differential equation (SDEs) [20]. In this case, one modification of equation (3.7) could be where $\mathbf{W}(t)={({W_{1}}(t),{W_{2}}(t),{W_{3}}(t),{W_{4}}(t))^{\top }}$ is a vector of standard Wiener processes or Brownian motions, $\boldsymbol{\mu }$ is some vector-valued function, and L is a compatible matrix. Depending on the structure of $\boldsymbol{\mu }$, there may or may not exist an explicit solution. A SDE that affords a closed form solution is the Ornstein-Uhlenbeck (OU) process. See [3] 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 [32] 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 (5.1) exists, inference in such a system is closely related to Bayesian filtering and smoothing, amenable to such methods as the Kalman filter and extensions [27]. 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.