New Approach to the SIR Inversion Problem: From the 1905-1906 Plague Outbreak in the Isle of Bombay To the 2021-2022 Omicron Surge in New York City

Discrete-Time SIR Model

Following the lead of other investigators [4, 35, 38-40], we adopt a discrete-time approach. In our SIR epidemic model, the time axis is marked off in equally spaced intervals t = 0, 1, …, T, where the duration of each interval is sufficiently small as to adequately approximate the classical, continuous-time version [5, 49-52]. At any time t, individuals within this closed population can be in one of three mutually exclusive states: susceptible (S), infected (I), or removed (R). The latter state, which includes both recovered living individuals and decedents, is assumed to be absorbing. To minimize possible complications arising from the non-identifiability of multiple parameters [53-56], we assume a fixed, demographically closed population of size N.

Let S_t, I_t, and R_t denote the respective numbers of individuals in each of the three states at time t. The dynamic path of the epidemic is governed by the following deterministic, coupled difference equations:

Equations (1a) through (1d) represent the well-known forward Euler approximation to the underlying continuous-time SIR model of coupled differential equations [57]. Apart from the population size N, this dynamic system has two parameters: β and γ. In equations (1a) and (1b), β >0 is an infection transmission parameter, while in equations (1b) and (1c), the parameter γ >0 gauges the proportion of infected individuals who transition to the removed state at each time period. The multiplicative term βS_t−1I_t−1/N in equations (1a) and (1b) reflects the law of mass action [58], whereby susceptible individuals become infected in proportion to their frequency of contact with currently infected individuals. All individuals within the population are assumed to mix homogeneously, with no subgroup of individuals mixing preferentially with any other subgroup.

The final equation (1d) reflects the constant size N of the population and is consistent with equations (1a) through (1c). Strictly speaking, one should adjust the size N of the mixing population in equations (1a) and (1b) to take account of removals by death. Unless the overall death rate is substantial, this adjustment is usually ignored in model implementations.

To further simplify matters, we assume that initially all individuals are either susceptible or infected. Specifically, at t = 0, we have:

In equations (2a) and (2b), the additional parameter i₀ represents the proportion of the entire population of size N that is initially infected. We address the more general case where R₀ >0 below.

In the dynamic system (1), the mean duration of infection is time-invariant and equal to 1/γ. Given the initial conditions (2), the system will result in an epidemic wave so long as (1 − i₀) β/γ >1, a well-known result known as the epidemic threshold theorem [24, 51, 59, 60]. Assuming that 1 − i₀ ≈ 1, most authors write this epidemic threshold condition as ℛ₀ >1, where ℛ₀ = β/γ is the basic reproduction number [24, 52, 61, 62].

Parameter Estimation Problem

We do not have direct observations on the underlying state variables S_t, I_t, and R_t. If we had such data, our inversion problem would border on trivial [39, 63, 64]. Instead, we observe only the reported counts of new infections at various intervals. For now, we assume that new infection counts y_t are observed at each discrete time t. Below we address the more general case where such counts are observed less frequently. The counts y_t represent observations on the output variables X_t, t = 1, …, T, which from (1) correspond to:

It will be helpful to define the parameter α = 1 − γ as the proportion of infected individuals who remain infected from one time period to the next, so that the probability an infected individual remains infectious after τ time periods is α^τ. Given the definition of the output variable X_t in (3), our dynamic system (1) can be redefined as:

So long as the state variables adhere to the condition that that S_t + I_t + R_t = N, an explicit difference equation for R_t is unnecessary.

We can now characterize our parameter estimation problem. Given observations y_t on the output variables X_t, we want to estimate the unknown parameters β, α, i₀, and N. Our stumbling block is that we cannot express X_t as a closed-form function of these parameters. We know only that the output variables X_t adhere to the dynamic system defined by (2), (3) and (4), which in turn depends on the parameters β, α, i₀, and N.

Solution Strategy: Nonlinear Least Squares

To proceed, we need some additional notation. Let y = (y₁, …, y_T)′ and X = (X₁, …, X_T)′, respectively, denote column vectors of the observed incidence data and the corresponding output variables at each time t, where we use boldface symbols denote vectors or matrices. Let θ = (β, α, i₀)′ denote the column vector of the unknown parameters excluding the population size N. For now, we condition on N, but as we discuss below, this additional parameter can be separately identified once θ has been estimated. Let X(θ) represent the functional dependence of the output variables on the parameters, suppressing for now that the fact that X also depends on N. Finally, let the operator D denote the gradient of partial derivatives with respect to the parameters. For example, the column vector represents the gradient of partial derivatives of the output variable X_t at time t.

We want to optimize some objective function V(y, X(θ)) with respect to θ. While X(θ) has no closed-form expression, the time-specific gradient vectors DX_t of partial derivatives of X_t with respect to the elements of θ follow a system of difference equations and can thus be readily calculated iteratively. The same goes for the Hessian matrices D²X_t of second-order partial derivates with respect to the elements of θ. Once we have computed the gradient vectors DX_t and Hessian matrices D²X_t, we can use the Newton-Raphson algorithm [46] or, under more restrictive conditions, the Gauss-Newton [47] algorithm to optimize V.

To see how this solution strategy works, we specifically consider the nonlinear least squares optimization criterion V(y, X(θ)) = (y − X(θ))′(y − X(θ)), which can be written in summation notation as:

This criterion has been widely used in attempts to fit the SIR model to incidence data [26, 41-44]. It is well known that minimizing V is equivalent to computing the maximum likelihood estimate of θ under the assumption that the observations y_t are independently normally distributed N(X_t(θ), σ²) with respective means X_t(θ) and homoscedastic variance σ². In Appendix A, we adopt a more general maximum likelihood framework and specifically consider an alternative Poisson distribution.

Calculating the Gradient of the Least Squares Optimization Criterion V

The gradient of the least squares criterion V in equation (5) with respect to the parameter vector θ is:

From (6), we learn that the gradient DV depends on the gradients DX_t. Our next step is to derive specific expressions for the latter gradient terms. Taking the derivative of (3) with respect to θ, we get: where the notation DV represents the unit column vector . Taking the derivative of (1a), we get:

From (1b) and (3), we can write .

Taking the derivative of this expression gives: where Dα similarly represents the unit vector (0,1,0)′.

We next need to compute the initial values of the gradients DS₀ and DI₀. Since I₀ = i₀N and S₀ = (1 − i₀)N, we can write the gradients DS₀ = −NDi₀ and DI₀ = NDi₀, where the gradient Di₀ simplifies to (0,0,1)′. Given these initial gradient values, we can use (8) and (9) to iteratively compute the vectors DS_t and DI_t for all t = 1, …, T. Once these gradients have been computed, we can then apply (7) to iteratively compute the corresponding gradients DX_t. Those quantities in turn yield the gradient DV of our objective function through equation (6).

Calculating the Hessian Matrix

Let DX denote the T × 3 matrix whose t-th row is the vector DX_t′. Taking the derivative of the gradient in (6), we obtain:

From (10), we similarly learn that the Hessian D²V depends on the gradients DX_t, which we have already derived, as well the Hessian matrices D²X_t.

To develop the corresponding iterative formulas for the Hessian matrices D²X_t, we need some additional notation. For two column vectors A and B of dimension L × 1, we define the outer product A · B as the L × L matrix AB^′ formed by the cross-products of the elements of A and B. (Note that B · A is the transpose of A · B.) We further define A⨀B = (A · B) + (B · A) as the symmetric matrix formed by the sum of the outer product and its transpose. In the three-dimensional case, A = (a₁, a₂, a₃)′, and B = (b₁, b₂, b₃)′, we have:

Now taking the derivative of DX_t as displayed in (7), we get:

The Hessian matrices D²X_t depend in turn on the Hessian matrices D²S_t and D²I_t. The corresponding expressions for these Hessian matrices are:

To complete our calculations, we note that the initial values D²S₀ and D²I₀ are simply null matrices.

Taken together, equations (10), (12), (13) and (14), along with the initial null values of D²S₀ and D²I₀, permit us to iteratively compute the Hessian D²V, just as we outlined for the gradient DV above.

Multiple Local Optima; Highly Correlated Parameter Estimates

In the Newton-Raphson [46] approach to optimization, which we employed here, the current value of the parameter vector θ^(k) is repeatedly mapped into an updated value θ^(k+1) according to the well-known rule:

In equation (15), the gradient DV, as defined in (6), and the Hessian D²V, as defined in (10), are both computed at the current value of the parameter vector θ^(k). The step size 0 < q ≤ 1 is under control of the programmer. This approach works flawlessly when the objective function V is globally convex, for example, when each output variable X_t is a linear function of θ. But that’s hardly the situation in the SIR inversion context.

To begin with, there is good reason to suspect that V is not globally convex and may have multiple local optima. While the classical SIR model predicts a single-peaked wave of incident cases X_t so long as ℛ₀ >1, the observed data y_t often display multiple peaks over time. (A good example is the two-peak plot of the 2001 Dengue fever outbreak in Havana [60].) In that case, the search algorithm embodied in (15) may easily end up at a local rather than a global optimum of V as it gets trapped in a region of the parameter space that fits one of the multiple peaks.

What’s more, there will be regions of the parameter space where the Hessian D²V is not positive definite because the second term in equation (11) (that is, will not necessary be a positive definite matrix. In that case, the updated parameter rule (15) may not result in a decrease in the objective function V as the algorithm veers away from the optimum.

These difficulties may be partly addressed by the Gauss-Newton algorithm [47], which in the current context entails the following parameter updating rule:

The idea here is that matrix 2 DX′DX, which corresponds to the first term in our expression for the Hessian D²V in (10), is always positive definite, so that the rule (16) will consistently result in a decrease in V. However, the matrix DX′DX may be ill-conditioned, with eigenvalues spanning multiple orders of magnitude, so that the steps θ^(k+1) − θ^(k) jump around uncontrollably and fail to converge to zero. While various modifications such the Levenberg-Marquardt algorithm have been proposed to recondition the matrix DX′DX [65], in the final analysis our search procedure may end up being exquisitely sensitive to the choice of the initial parameter vector θ⁽⁰⁾.

Even when our objective function V has an interior global minimum, where it is at least locally convex in the neighborhood surrounding the optimum, we may still encounter an inverted ridge or ravine where the parameters are highly correlated. There may indeed be a unique value θ^* at the very bottom of the ravine that achieves the global minimum V^*. But along the floor of the ravine, there is a one-dimensional curve containing θ^* where V is almost equal to V^*. Strictly speaking, the parameters are structurally identified, but they are not practically identified [66].

To confront these problems of multiple equilibria and correlated parameters in the applications below, we took care to construct and inspect maps of the least-squares criterion V and the gradient DV as functions of the parameters. We also examined the path of successive updates θ^(k) in parameter space generated by our Newton-Raphson search algorithm.

Confidence Intervals; Reparametrizing

Our reliance on the nonlinear least squares criterion V in equation (5) permits us to estimate confidence intervals for the estimated parameters [48]. Conditional upon N, the variance-covariance matrix of the parameters θ = (β, α, i₀) is C = s²(DX′DX)⁻¹, where s² = V/T, and where DX and V are evaluated at the optimum. The standard errors are the square roots of the diagonal elements of C. The symmetric 95 percent confidence intervals, based upon the assumption of an asymptotic normal distribution, can be evaluated as ±1.96 standard errors about the estimates. We can then use the Delta method [67] to compute the standard errors and confidence intervals around nonlinear functions of the parameters, such as the basic reproduction number ℛ₀ = β/γ = β/(1 − α) and the mean duration of infection 1/γ = 1/(1 − α).

In some applications where the data y_t exhibit substantial variability and the resulting value of s² is large, the estimated symmetric confidence interval surrounding α may exceed the allowable range of 0 ≤ α < 1. In such cases, we can reparametrize, instead specifying θ = (β, ω, i₀) with α(ω) = 1/(1 + e^−ω), where the substitute parameter ω is unbounded. Our expression for DI_t in equation (9) will remain valid, but with the modification that . Our expression for D²I in equation (14) will become:

The additional term on the right-hand side of equation (17), which does not appear in (14), contains the matrix , where .

Recovering the Population Size Parameter N

While X(θ) has no closed-form expression in terms of θ, it turns out that each element X_t is a linear function of the population size parameter N. That is, X_t can be written in the form φ_t(θ)N, where φ_t(θ) is a function of the remaining parameters θ.

To prove this result, we only need to show that the state variables S_t and I_t are themselves proportional to N. Thus, if for all times t = 1, …, T, we can show that S_t = f_t(θ)N and I_t = g_t(θ)N for some functions f_t(θ) and g_t(θ), then from equation (3), we would have X_t = βS_t−1I_t−1/N = βf_t−1(θ)g_t−1(θ)N. Thus, φ_t(θ) = βf_t−1(θ)g_t−1(θ).

We can prove that S_t and I_t are proportional to N by mathematical induction. From equations (2a) and (2b), respectively, we known that S₀ = (1 − i₀)N and I₀ = i₀N. So, the proposition is true for t = 0. Now suppose that S_t = f_t(θ)N and I_t = g_t(θ)N, at time t. We claim that S_t+1 = f_t+1(θ)N and I_t+1 = g_t+1(θ)N necessarily hold for some functions f_t+1(θ) and g_t+1 (θ). From (1a), we can write . We have , and so . Similarly, , and so g_t+1 (θ) = g_t(θ)α + βf_t(θ). ▪

The fact that each output variable X_t is proportional to N permits us to employ an iterative procedure resembling the expectation-maximization (or EM) algorithm [68] to recover an estimate of that parameter. Let’s say that we have a parameter estimate N⁽ⁿ⁾ at iteration n of the algorithm. We employ the previously described procedure to estimate θ⁽ⁿ⁾ conditional upon N⁽ⁿ⁾. That corresponds to the maximization (or M) step. (From the coding standpoint, when apply the Newton-Raphson parameter search rule described in (15) to estimate θ⁽ⁿ⁾, we’re looping on k within a loop on n.) Given both N⁽ⁿ⁾ and θ⁽ⁿ⁾, we have output-variable estimates for all t. We now compute:

That is, we compute an adjustment factor κ⁽ⁿ⁺¹⁾ by regressing y on X⁽ⁿ⁾. For the expectation (or E) step, we update the population size as N⁽ⁿ⁺¹⁾ = κ⁽ⁿ⁺¹⁾N⁽ⁿ⁾. This iterative procedure will work for any optimization criterion V and not just for the least squares criterion described above.

Non-Identifiability of the Initially Resistant Population

We have thus far assumed in equation (2c) that no one in the population was initially resistant to infection, that is, R₀ = 0. Now let’s assume more generally that a non-negative fraction 1 >r₀ ≥ 0 of the population is already resistant at t = 0. We thus replace the restricted initial conditions (2) with the following more general initial conditions:

So long as the R₀ initially resistant individuals mix homogeneously with those in the susceptible and infected states, equation (3) governing the output variables X_t as well as the equations of motion (4) for the state variables S_t and I_t remain unchanged.

Apart from our original parameter vector θ = (β, α, i₀) and the population size N, it appears that we now have an additional parameter r₀ to be estimated. It turns out, however, that r₀ cannot be identified in our discrete SIR model from the observed data y alone [54].

To see why, let’s suppose that we estimate our more general model conditional upon some fixed value of r₀. We denote the resulting minimum least squares estimates as θ^*(r₀) and N^*(r₀) to show their dependence on the value of r₀ chosen. We denote the corresponding individual components of θ^*(r₀) as β^*(r₀), α^*(r₀), and . Our original parameter estimates conditional upon r₀ = 0 thus correspond to β^*(0), α^*(0), , and N^*(0). Then the following conditions hold for all 1 >r₀ ≥ 0:

What’s more, the paths of the output variable X_t and the state variables S_t and I_t at the optimum values of the parameters, but not the state variable R_t, will be independent of r₀. Accordingly, the optimum value of V will likewise be independent of r₀.

We again rely upon mathematical induction to prove that the state variables S_t and I_t at the optimum are independent of r₀. First consider t = 0. For any value of r₀, we have from (19a) that . Applying (20c) and (20d), we get , which does not depend on r₀. Similarly for any value of r₀, we have from (19b) that . Again applying (20c) and (20d), we get , which is likewise independent of r₀. Now assume that S_t and I_t are independent of r₀. We show that S_t+1 and I_t+1 must also be independent of r₀. From (1a), we have S_t+1 = S_t(1 − β^*(r₀)I_t/N^*(r₀)). Applying (19a) and (19d), this expression resolves to S_t+1 = S_t(1 − β^*(0)I_t/N^*(0)). From (1b) and (3), we have I_t+1 = I_t(α + β^*(r₀)S_t/N^*(r₀)). Again applying (20a) and (20d), this expression similarly resolves to I_t+1 = I_t(α + β^*(0)S_t/N^*(0)). ▪

Interpreting the Estimated Population Size N: Underreporting and Incomplete Mixing

While the additional parameter r₀ is not identifiable from the data y alone, we might still be able to identify it from other data. Suppose, for example, that we had sharp prior information on the basic reproduction number ℛ₀ = β/γ = β/(1 − α). Let’s denote this prior estimate . We first compute the basic reproductive number implied by our parameter estimates conditional upon r₀ = 0. So long as , we can then use our prior information to conclude that .

Why couldn’t we similarly take advantage of prior information on the population size to ascertain the initial proportion r₀ of resistant individuals? Suppose we knew from census data that the population contained individuals. We would first estimate N^*(0) conditional upon r₀ = 0, and then, so long as , we would have the estimate . Unfortunately, reliance solely on census data for prior estimates of population size N is complicated by two other phenomena: underreporting and incomplete mixing.

In many contexts, particularly in recent applications of the SIR and related compartmental models to COVID-19 incidence, it has been widely recognized that a significant number of incident cases may have gone unreported [69]. In the absence of reliable information on the temporal pattern of such underreporting, the most parsimonious approach to this phenomenon has been to assume a case identification ratio equal to a constant p < 1 [70-72]. In that case, our model would need to be modified to accommodate the reality that our reported incidence data y_t are in fact estimates of pX_t rather than X_t. We’ve already learned, however, that X_t = φ_t(θ)N, where φ_t(θ) is a function of the remaining parameters θ, and this is the case even in our more general model admitting r₀ >0. Accordingly, our reported incidence data y_t are really estimates of pX_t = φ_t(θ)(pN), and the estimated population size derived from our model is really an estimate of pN.

We thus have a knotty problem of confounding. We can estimate population size N^*(0) assuming that there are no initially resistant individuals and no underreporting. If we had census data , we could account for both phenomena, writing . Without more information, we cannot separately identify p and r₀.

We would ordinarily interpret the parameter N to gauge the size of the population at risk for contagion. This population would consist of all individuals who homogeneously mix with each other in accordance with the law of mass action embodied in equations (1) and (3). Many investigators, however, have properly recognized that the underlying assumption of homogeneous mixing may not apply to the entire population [19, 73-78].

Let’s say we’re analyzing an outbreak on a college campus with known student enrollment . When we run our model with r₀ = 0, we obtain an estimate N^*(0) that is, say, only about 10 percent of . This finding does not necessarily imply that 90 percent of the student body was resistant or that only one in ten cases was reported. Instead, it may mean that only small fraction of the student body was directly involved in the mixing process that generated the outbreak. A good example would be the COVID-19 outbreak at the campus of the University of Wisconsin-Madison in September 2020, where total student enrollment was , but where the large fraction of cases was concentrated in two on-campus student residence halls with a combined population of about 2,932 [79].

Any comparison between the estimated population size N^*(0) and the census statistic thus entails three confounded interpretations. Some cases have gone unreported. Some individuals may be initially resistant. And other individuals may be susceptible but remain outside the core population through which the infectious agent has spread.

Adapting the Model to Data on Removals

Our analysis has thus far assumed that each data point y_t is an observation on the corresponding output variable X_t capturing the number of individuals who transitioned from the susceptible to the infected state at time t. Alternatively, we can construct an analogous model where the data point y_t is instead an observation on a distinct output variable Z_t capturing the number of individuals who transitioned from the infected to the removed state at time t. This approach would be important if we had data on the number of deaths during an epidemic of a disease with a near 100 percent fatality rate, as KM observed in the case of 1905-1906 plague outbreak on the Isle of Bombay [5].

In this alternative model, the output variable Z_t is:

Continuing with the least squares framework, the objective function is now:

To minimize this objective function employing the same strategy above, we will need to compute the gradient vector DZ_t and the Hessian matrix D²Z_t. From (25), we immediately have:

The corresponding Hessian matrix is:

To compute these two expressions, we will need the gradient DV_t and Hessian D²V_t, but these respective quantities were already computed in equations (9) and (14) above. To complete our estimation procedure, we note that the output variables are similarly linear functions of the population size parameter N, so that the EM-type algorithm described above likewise applies.

Working with Aggregate Data

In some applications, we may have only aggregate data reported over a coarse time scale, rather than the finely divided time axis assumed so far. For concreteness, let’s assume that the underlying SIR epidemic model is valid when the time axis t = 0, 1, …, T is marked off in days, but we only have data y_m on the cumulative number of deaths during each 7-day week, indexed by m = 1, …, M. Days t are mapped into weeks m by the relation , where the floor operator ∥x∥ maps into the largest integer j such that j ≤ x. We further assume that T = 7M, so that week M, which corresponds to the last week, ends on day T.

Now define the M × T aggregation matrix W with element w_mt = 1 if and only if m = h(t). Otherwise, w_mt = 0. While our finely disaggregated model generates daily values Z_t of the output variable from the underlying parameters θ, our data y_m represent observations on the aggregated values , where w denotes row m of the matrix W and where Z denotes the T × 1 column vector with coordinates Z_t. (From the computational standpoint, we’re calculating the integral required to convert a continuous- to a discrete-time SIR model [40].) In vector notation, our least squares criterion becomes:

Let DZ denote the T × 3 matrix whose row t is the gradient DZ_t, and let W(DZ) denote the M × 3 matrix derived by multiplying the M × T matrix W by the T × 3 matrix DZ. The gradient of V is then given by:

Let w denote row m of the matrix W, so that represents the predicted value. The corresponding Hessian matrix becomes:

Accordingly, we can run our SIR model to generate the state variables Z_t and then calculate the derivatives DZ_t and D²Z_t in the disaggregated time scale and then apply the aggregation matrix W to weight the results.

Omicron Wave, December 2021 – March 2022, New York City

We studied the reported daily incidence of COVID-19 during the SARS-CoV-2 Omicron/BA.1 wave of December 2021 – March 2022 in New York City, NY, United States, a city of population 8.49 million. Our data consisted of daily counts of cases reported by the New York City department of health [80], where the date of report was intended to be the date when a positive test was performed or when the diagnosis of COVID-19 was otherwise made.

Our data showed systematic variation in case counts by day of the week, with many fewer cases diagnosed over the weekends. To account for these fluctuations, we converted the raw case counts c_t into centered 7-day moving averages, that is, , where t indexes the date of report. Figure 1 shows the raw counts of reported cases per day c_t (connected gray datapoints) as well as the daily case counts y_t adjusted for the day of the week (red datapoints). We relied upon the adjusted daily counts y_t to estimate our SIR model parameters.

• Download figure
• Open in new tab

Figure 1. Daily Reported Cases of COVID-19, 12/1/2021 – 3/15/2022, in New York City.

The connected gray datapoints show the raw case counts (c_t). The red datapoints, covering 12/4/2021 – 3/12/2022, show the centered 7-day moving averages (y_t).

Plague Outbreak, December 1905 – July 1906, Isle of Bombay

Figure 2 displays the average daily number of reported deaths from plague in the Isle of Bombay during each week over a 31-week period, beginning with the week starting Sunday, December 17, 1905, and continuing through the week starting Sunday, July 15, 1906. The average daily deaths were computed from the weekly totals plotted by KM in their classic 1927 paper [5]. The height of each bar is the average daily number of deaths during that week. The area of each bar equals the reported number of deaths during that week.

• Download figure
• Open in new tab

Figure 2. Average Number of Daily Deaths from Plague, Isle of Bombay, Week Starting 12/17/1905 Through the Week Starting 7/15/1906.

Each bar corresponds to one week. The height of each bar is the average daily number of deaths during that week. The area of each bar equals the number of deaths during that week. Data adapted from Kermack and McKendrick [5].

In our analysis of the plague death data, we interpreted the weekly death counts y_m, indexed by m = 1, …, 31, as observations on the aggregated values , where Z_t are the underlying output variables with daily index t = 1, …, 217, defined in equation (21), and where w_mt are elements of the aggregation matrix W described above. We estimated the parameters of our SIR model by minimizing the least squares criterion V defined in equation (25).

Data Processing and Computation

The raw COVID-19 case counts (c_t) were downloaded from the New York City data portal [80] in comma-separated-value (CSV) format and converted via Stata Statistical Software Release 17 [81] into internal Stata (DTA) format. The 7-day moving averages (y_t), displayed in Figure 1, were computed via Stata programming language [82]. The data on weekly deaths from plague in the Isle of Bombay were taken from the graph on page 714 of the original KM paper [5]. All code for parameter estimation described in the Statistical Methods section above was written in Mata [83], a matrix programming language embedded within Stata. Calculations were carried out on a MacBook Pro with a 2.3 GHz 8-Core Intel Core i9 processor.

Omicron Wave, December 2021 – March 2022, New York City

Figure 3 shows the predicted values of the output variable X_t as a connected curve superimposed on the observed datapoints y_t for the Omicron wave in New York City. The computation time to achieve convergence of our Newton-Raphson algorithm was 0.93 sec.

• Download figure
• Open in new tab

Figure 3. Daily Reported and Predicted Cases of COVID-19, 12/4/2021 – 3/12/2022, New York City.

The red points show the data y_t derived from Figure 1. The curve connects the predicted values of the output variable X_t.

Table 1 summarizes the resulting parameter estimates. At t = 0, corresponding to December 3, 2021, the infected proportion i₀ was an estimated 0.8 percent. The estimated population-size N was about 12 percent of the city’s total population. The estimated basic reproduction number ℛ₀ was on the order of 4. The estimated mean duration of infectivity 1/(1 − α) was in the range of 2 to 3 weeks. The 95% confidence intervals were estimated conditional upon N, as described above.

Figure 4A plots the least squares criterion V (blue curve, left axis) and the first partial derivative ∂V/∂β (red curve, right axis) as functions of the parameter β. The remaining parameters have been held constant at their estimated values. The criterion V reaches a minimum at the optimum β = 0.233, where ∂V/∂β = 0. The function V is convex in the interval from β = 0.169, where ∂V/∂β reaches a minimum, to β= 0.352, where ∂V/∂β reaches a maximum.

• Download figure
• Open in new tab

Figure 4.

Panel A. Least Squares Criterion V (Left Axis) and First Partial Derivative ∂V/∂β (Right Axis) as Functions of the Parameter β. Panel B. Least Squares Criterion V (Left Axis) and First Partial Derivative ∂V/∂α (Right Axis) as Functions of the Parameter α.

Figure 4B above displays the analogous plot of V and ∂V/∂α as functions of the parameter α, where the remaining parameters are similarly held constant at their optimum values. The criterion V reaches a minimum at the optimum α = 0.941, at which point ∂V/∂α = 0. The function V is convex in the interval from α = 0.868, where ∂V/∂α reaches a minimum, to α = 1, the boundary of admissible values of α, where ∂V/∂α remains positive.

Figure 5 plots the projection of V onto the (α, β) plane. Again, the remaining parameters (i₀, N) were held at the optimum values given in Table 1 above. The darkest area represents a ravine where V attains its lowest values. The yellow point in the center is the optimum (α, β) = (0.941, 0.233) where V is minimized.

• Download figure
• Open in new tab

Figure 5. Least Squares Criterion V Projected onto the (α, β) Plane.

The yellow point identifies the minimum where (α, β) = (0.941, 0.233). The parameters i₀ and N were held constant at their optimum values given in Table 1. The plot was produced in part from the Stata program surface [84].

Plague Outbreak, December 1905 – July 1906, Isle of Bombay

Figure 6 reports our analysis of data on deaths from plague during December 1905 – July 1906 in the Isle of Bombay, originally reported and studied by KM [5].

• Download figure
• Open in new tab

Figure 6. Reported and Predicted Deaths from Plague, Isle of Bombay, 12/17/1905 – 7/21/1906.

The time axis is marked off in weeks, with week m = 1 starting on 12/17/1905 and week m = 31 starting on 7/15/1906. As noted in Figure 2, the height of each bar is the average daily number of deaths reported each week, while the area of each bar equals the number of deaths during that week, adapted from Kermack and McKendrick [5]. The green curve represents the authors’ original fit to the data, while the red curve represents our estimates of the output variable Z_t.

As in Figure 2 above, we have displayed the raw data on deaths as a bar graph, where each vertical bar covers one week, indexed m = 1, …, 31, where the height of each bar is the reported average daily deaths that week, and where the area of each bar is the reported number of weekly deaths (y_m). The superimposed green curve represents the fit to the data, based upon the KM’s approximation , where sech represents the hyperbolic secant, and where we have divided their published formula by 7 to convert their estimates into daily units. The green curve thus corresponds to the authors’ estimates of the underlying output variable Z_t.

The red curve in Figure 6 connects our predictions of the daily values Z_t for t = 1, …, 217. The underlying parameter estimates were: β = 0.610, with 95% confidence interval (0.537, 0.683); α = 0.44 (0.376, 0.515); i₀ = 2.80 × 10⁻⁵ (1.26 × 10⁻⁵, 4.34 × 10⁻⁵); and N = 5.103 × 10⁴. The estimated basic reproduction number was ℛ₀ = 1.101 (1.094, 1.107), while the estimated mean duration in the infected state was 1/(1 − α) = 1.804 (1.579, 2.031).

Additional Results

In addition to our main results above, we performed various ancillary analyses, which are reported in appendices. Appendix B provides an approach to computing the unconditional confidence intervals around the parameter estimates for the New York City Omicron wave. In Appendix C, we display the two- and three-dimensional paths of successive parameter estimates derived from our application of the Newton-Raphson algorithm to the New York City data. In Appendix D, we report the results of a robustness test, where we re-estimate our SIR model for New York City after multiplying the original case counts y_t by mean-preserving, lognormally distributed noise.

In Appendix E, we account for the observation that some of the SARS-CoV-2 cases reported during the rising phase of the New York City epidemic wave of Figure 1 were in fact infections by the Delta variant, while some of the cases reported during the declining phase were in fact Omicron BA.2 infections. To that end, we reran our SIR model after multiplying the original case counts y_t by the estimated proportion of Omicron BA.1 cases at each date t. We also reran our SIR model on the original case counts y_t, but restricted the observation interval to those dates t during which the BA.1 variant was estimated to comprise at least 80% of all cases.

In Appendices F and G, we report estimates of our SIR model for two other jurisdictions during the COVID-19 epidemic. Appendix F shows our estimates for the SARS-CoV-2 Omicron wave in Los Angeles County, CA, during the same time period covered by our study of New York City above. Appendix G shows our estimates based upon an outbreak of SARS-CoV-2 at the University of Wisconsin-Madison during September 2020, originally reported in a study of the potential super-spreader influence of a nearby cluster of local bars [79].

Finally, Appendix H provides a diagnostic plot of the least squares criterion V as a function of the parameter β for the Isle of Bombay plague data, comparable to Figure 4A for New York City above.

Is SIR the Correct Structural Model for Omicron?

The mere fact that we have found a workable solution to the SIR inversion problem does not necessarily mean that the Susceptible-Infected-Removed model is the most appropriate description of the data under study. Our illustrative analysis of the SARS-CoV-2 Omicron/BA.1 wave in New York City during December 2021 – March 2022 cogently brings home the point.

Figure 3 shows a tight fit of the SIR model-predicted curve X_t to the case incidence datapoints y_t for New York City. Figures 4A, 4B, and 5 confirm that the estimated values of the parameters β and α are indeed situated at the minimum of a convex region of our least squares criterion function V. Appendix C confirms that successive updates of the parameter vector θ^(k) generated by the Newton-Raphson algorithm follow a path along the convex surface of V. Appendix D confirms that our estimates are robust to the inclusion of noise in the observed data.

Despite these indicators that the underlying model was not somehow ill-conditioned, the mean duration of infectivity, computed as 1/(1 − α), was estimated from the New York City data to be on the order of 17 days, with a 95% confidence interval of two to three weeks. That’s way out of line with what is known from direct clinical measurement. In one cohort of 55 symptomatic Omicron-infected patients, only 13.5 percent continued to shed virus ten days after infection [85]. Yet our New York City-based parameter estimate for α would give the proportion remaining infectious after ten days at α¹⁰ = 54.4%.

To be sure, some of our alternative analyses yielded lower estimates of the mean duration of infectivity. When we attempted to exclude Delta and BA.2 infections from the original New York City case counts y_t, the estimated mean duration dropped to about 10 days (Appendix E). When we reran our SIR model on the Omicron wave in Los Angeles County, the estimated mean duration was about 9 days (Appendix F). But these estimates are still don’t come close to those based on the available clinical data. When we reran our SIR model on an outbreak at the University of Wisconsin-Madison in September 2020, we obtained an estimated mean duration of about 7 days (Appendix G). But that outbreak was quite likely driven by the ancestral strain of SARS-CoV-2, and certainly not by Omicron.

Despite its remarkable fit to the data, SIR may thus fail as a structural model of the Omicron wave, even if it apparently succeeds as a reduced form model [86]. While its parameters β and α can indeed be backed out from the data and the inversion problem solved, these parameters do not necessarily warrant the structural interpretation that we assumed in our exposition of the SIR model in equations (1) through (5) above. That is not to say, however, that an adequate reduced form model is incapable of making accurate projections [87].

We should not, however, abandon the possibility that SIR may indeed be the correct structural model of the Omicron wave. To the contrary, the parameter estimates may reveal a critically important characteristic of the Omicron variant that clinical studies, which are necessarily biased toward symptomatic patients, have so far missed. In effect, the Omicron wave was so massive not because the variant had such a high per-contact infectivity (through the parameter β), but because a much larger proportion of asymptomatic infected individuals remained persistently infectious (through the parameter α). We also need to consider that the degree of persistent infectivity is not determined solely by biological characteristics of the virus, but also by the extent to which infected individuals isolated themselves from others. That behavioral component of the parameter α may have changed during the Omicron wave.

It might be argued that there is a third plausible interpretation, namely, that SIR may indeed be the correct structural model, but that the case incidence data analyzed here were inadequate to separately identify the key parameters β and α. That is, our four-parameter version of SIR is what some physicists have described as a “sloppy” model [53, 55]. Indeed, the dark ravine in the surface map of V in Figure 5 points to a high correlation between the two parameters in the topological neighborhood of the global optimum in the (β, α) plane, an observation that has been noted previously [71]. That finding suggests that the basic reproduction number ℛ₀ = β/γ = β/(1 − α) is the only identifiable or “stiff” parameter.

Still, the results in Appendix C, which follows the path of the Newton-Raphson search algorithm in parameter space, suggest that β and α are indeed separately identified from the available count data. In Figure C2, as the algorithm approaches convergence to the minimum of the criterion V, the implicit derivative dα^(k)/dβ^(k) does not approach a constant value, as we would expect in the case of non-identification.

What About the Results for the Plague?

Not only do our SIR predictions fit reasonably well to the observed data on Bombay plague deaths in Figure 6, but they virtually coincide with the curve of the approximate closed-form solution drawn by KM [5]. While others have fitted SIR-type compartmental models to the 1905–1906 Bombay plague data [88], our study appears to be the first replication of the authors’ century-old result.

The plague, caused by the bacterium Yersinia pestis, can spread through human populations by various transmission pathways. Bubonic plague, generally thought to be the cause of a series of Bombay outbreaks that began in 1896 [89-91], is transmitted to humans through rat fleas [92], but transmission between humans also appears to occur via infected human fleas or body lice [93]. Pneumonic plague, by contrast, is transmitted by direct human-to-human transmission via respiratory droplets and can occur as a complication of bubonic plague [94, 95]. Pneumonic plague is more lethal, but transmission requires closer person-to-person contact and is thus thought to have a lower basic reproduction number ℛ₀ [94].

So, does our SIR model similarly succeed as a reduced form model but fail as a structural model of the Bombay plague outbreak? Several structural compartmental models of plague transmission have been tested, typically involving separate states and parameters for humans and vectors [88, 89, 95-97]. In a rat-flea transmission model for bubonic plague, for example, the structural model contained a separate sub-epidemic module for rats with a distinct basic reproduction number for rat-to-rat transmission [97]. A human-flea-human structural model appeared to have the best fit to nine plague outbreaks during Europe’s Second Pandemic from 1348 to 1813 [95]. Our minimalist SIR model has none of these structural features.

These general observations, however, hardly exclude the possibility that the pneumonic form may have been responsible for the 1905-1906 outbreak studied by KM. In his classic 1926 Treatise on Pneumonic Plague [98], Wu Lien-Teh acknowledged the 1908 Indian Plague Commission’s conclusion that pneumonic constituted less than 3 percent of cases overall. (See [98] at p. 9.) But he did not hesitate to point out that some “well-characterized pneumonic outbreaks are on record.” (p. 90-92.) “In fact,” he added, “it was this great incidence of respiratory diseases in Bombay which attracted Childe’s attention and led to his finally establishing the pneumonic form of plague as a distinct entity.” (pp. 11-12). He further noted an observation bias in favor of bubonic plague due to the rapid mortality from pneumonic “…because, as a rule, such patients succumb so quickly that they seldom reach the hospital.” (p. 269).

What’s more, the parameter estimates derived from our solution to the inversion problem are strikingly consistent with what is known about pneumonic plague transmission. Lien-Teh’s data suggested a median duration of illness of about two days from the onset of symptoms to inevitable death from heart failure ([98] at p. 249, 260). Our point estimate of the parameter α = 0.446, implying an average infectious period of 1.8 days, is similarly consistent with more contemporary estimates of the parameters of pneumonic human-to-human transmission [95]. Lien-Teh further stressed a very low early infectivity due to the initial lack of cough (p. 296). Transmission through contact at close range, he stressed, appeared to have been required (p. 299). His observations, along with more contemporary estimates of pneumonic human-to-human transmission [94, 95], are consistent with our estimate of a basic reproduction number ℛ₀ = 1.1.

The plague may indeed have been initially imported into Bombay in its bubonic form through rat-flea-human transmission, and the disease may have been sustained between outbreaks either by repeated importations or by an urban reservoir of selectively immune rats [89]. Still, our data are consistent the intriguing hypothesis that responsive isolation measures by British authorities [90] may have altered the mode of transmission during the 1905–1906 outbreak from the bubonic rat-flea-human mode to the pneumonic human-human mode, a phenomenon observed in the 17^th century when the Derbyshire village of Eyan fell victim to the Black Death [99]. If so, our parsimonious SIR model is in fact the appropriate structural model of the outbreak.

Why are the Estimated Population Sizes N So Small?

Table 1 reports estimated population size parameters of N = 1.013 million for New York City, while Appendix Table F1 gives N = 1.046 million for Los Angeles County. As noted, these estimates of N were only 11.9 percent of the total census population of New York City and 14.5 percent of the corresponding census population of Los Angeles County. Where, then, did all the other millions of inhabitants go?

As discussed above, we cannot distinguish between three possible explanations of the shortfall. First, there is concrete evidence of significant underreporting of Omicron infections, due principally to the widespread availability of home rapid antigen testing [69], particularly in New York City [100, 101]. Second, despite the evidence of immune escape on the part of the Omicron variant [102], a substantial fraction of the population may still have retained long-term cellular immunity, particular through the administration of multiple vaccine doses [103]. Third, contrary to the underlying assumption of homogeneous mixing, there is strong evidence that a substantial proportion of the population avoided retail establishments, drinking and eating places, transportation venues, worksites and other high-risk locations [104].

Our analysis of the Bombay outbreak data gave an estimated size parameter N = 51 thousand individuals, which comes to about 6.5 percent of estimated Bombay population of 775 thousand in 1900 [105]. The same three factors could well be responsible for the shortfall. In particular, there is significant evidence for the acquisition of at least humoral immunity against Yersinia pestis among plague survivors [106].

Limitations and Extensions

Parameter search algorithms such as Newton-Raphson [46] and Gauss-Newton [47] have superior performance when they rely on closed-form expressions for the gradient vector and matrix of second derivates, rather than on numerical approximation [45]. Still, as Figures 4A and 4B teach us, there will nonetheless be a restricted region of the parameter space where the least squares criterion function is convex. Some of these regions may contain local rather than global optima. To avoid such problems of implementation, we plotted the least squares criterion V and its first partial derivatives as functions of the parameters. Still, finding the right initial values and keeping the parameter search within bounds remain unavoidable challenges.

As seen in Figure 1A and Appendix Figures F1 and G1, striking day-of-the-week effects on COVID-19 testing and case reporting required us to pretreat the raw case counts. Our centered, 7-day moving average appeared to be the most flexible nonparametric approach to this task. Alternative approaches that incorporate parametric models of day-of-the-week effects directly in our model of the output variable X_t have yet to be tried.

We have focused sharply on the original SIR model, rather than its numerous variations. Still, our basic approach can be extended to these more complex models. Our findings offer a caution, however, that such models as SEIR, which require an additional parameter governing the transition from an intermediate exposed state to the infected state, may very well turn out to be sloppy [71].

One exemption may be the well-studied SIRS model [59, 107], where individuals in the R (recovered) state can transition back to the S (susceptible) state as a result of waning immunity. The SIRS model still has three states, and the definitions of the output variables X_t and Z_t in equations (3) and (21), respectively, remain unchanged, but the parameter vector θ has an additional component governing the rate of transition from R back to S. Since the SIRS model is known to admit oscillations with a stable endemic equilibrium [59, 108], we would be in a position to run our parameter recovery algorithm on multiple waves of data y_t. That is a task for future research.