Impacts of Vaccination and Severe Acute Respiratory Syndrome Coronavirus 2 Variants Alpha and Delta on Coronavirus Disease 2019 Transmission Dynamics in Four Metropolitan Areas of the United States

Abstract

To characterize Coronavirus Disease 2019 (COVID-19) transmission dynamics in each of the metropolitan statistical areas (MSAs) surrounding Dallas, Houston, New York City, and Phoenix in 2020 and 2021, we extended a previously reported compartmental model accounting for effects of multiple distinct periods of non-pharmaceutical interventions by adding consideration of vaccination and Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) variants Alpha (lineage B.1.1.7) and Delta (lineage B.1.617.2). For each MSA, we found region-specific parameterizations of the model using daily reports of new COVID-19 cases available from January 21, 2020 to October 31, 2021. In the process, we obtained estimates of the relative infectiousness of Alpha and Delta as well as their takeoff times in each MSA (the times at which sustained transmission began). The estimated infectiousness of Alpha ranged from 1.1x to 1.4x that of viral strains circulating in 2020 and early 2021. The estimated relative infectiousness of Delta was higher in all cases, ranging from 1.6x  to 2.1x. The estimated Alpha takeoff times ranged from February 1 to February 28, 2021. The estimated Delta takeoff times ranged from June 2 to June 26, 2021. Estimated takeoff times are consistent with genomic surveillance data.

Appendix

1.1 Imputation of Missing Daily Case Counts

By October 31, 2021, many regions in the US were not reporting new detected COVID-19 cases on a strictly daily basis. When one or more daily case counts were not available, we imputed daily case counts on the basis of a linear fit to the two nearest available cumulative case counts. This approach has the effect of evenly distributing case counts across the days for which daily reports are unavailable.

1.2 Equations of the Compartmental Model

The compartmental model, which is illustrated in Fig. 1 and Fig. 7 in Appendix, consists of the following 40 ordinary differential equations (ODEs):

$$\begin{aligned} \frac{d{S}_{M}}{dt} &=-\beta {Y}_{\theta }\left(t\right)\left(\frac{{S}_{M}}{{S}_{0}}\right) \left({\phi }_{M}\left(t,\rho \right)+{m}_{b}{\phi }_{P}\left(t,\rho \right)\right)\\ &\quad -{{U}_{\sigma }\left(t\right)\Lambda }_{\tau }\left(t\right)\left[{P}_{\tau }\left(t\right){S}_{M}-\left(1-{P}_{\tau }\left(t\right)\right){S}_{P}\right]-\mu \left(t\right){S}_{0}\left(\frac{{S}_{M}}{{\phi }_{V}(t)}\right) \end{aligned}$$

(1)

$$\begin{aligned} \frac{d{S}_{P}}{dt}&=-{m}_{b}\beta {Y}_{\theta }\left(t\right)\left(\frac{{S}_{P}}{{S}_{0}}\right)\left({\phi }_{M}\left(t,\rho \right)+{m}_{b}{\phi }_{P}\left(t,\rho \right)\right)\\ &\quad +{U}_{\sigma }\left(t\right){\Lambda }_{\tau }\left(t\right)\left[{P}_{\tau }\left(t\right){S}_{M}-\left(1-{P}_{\tau }\left(t\right)\right){S}_{P}\right]-\mu \left(t\right){S}_{0}\left(\frac{{S}_{P}}{{\phi }_{V}(t)}\right) \end{aligned}$$

(2)

$$\begin{aligned} \frac{d{E}_{1,M}}{dt}&=\beta {Y}_{\theta }(t)\left(\frac{{S}_{M}}{{S}_{0}}\right)({\phi }_{M}(t,\rho )+{m}_{b}{\phi }_{P}(t,\rho ))-{k}_{L}{E}_{1,M}\\ &\quad -{U}_{\sigma }(t){\Lambda }_{\tau }(t)\left[{P}_{\tau }(t){E}_{1,M}-\left(1-{P}_{\tau }(t)\right){E}_{1,P}\right] \end{aligned}$$

(3)

$$\begin{aligned} \frac{d{E}_{1,P}}{dt}&={m}_{b}\beta {Y}_{\theta }(t)\left(\frac{{S}_{P}}{{S}_{0}}\right)({\phi }_{M}(t,\rho )+{m}_{b}{\phi }_{P}(t,\rho ))-{k}_{L}{E}_{1,P}\\ &\quad +{U}_{\sigma }(t){\Lambda }_{\tau }(t)\left[{P}_{\tau }(t){E}_{1,M}-\left(1-{P}_{\tau }(t)\right){E}_{1,P}\right] \end{aligned}$$

(4)

$$\frac{d{E}_{i,M}}{dt}={k}_{L}{E}_{i-1,M}-{k}_{L}{E}_{i,M}-{k}_{Q}{E}_{i,M}-{U}_{\sigma }\left(t\right){\Lambda }_{\tau }\left(t\right)\left[{P}_{\tau }\left(t\right){E}_{i,M}-\left(1-{P}_{\tau }\left(t\right)\right){E}_{i,P}\right], \quad \hbox{for}\,\, i=2, 3, 4, 5$$

(5)

$$\frac{d{E}_{i,P}}{dt}={k}_{L}{E}_{i-1,P}-{k}_{L}{E}_{i,P}-{k}_{Q}{E}_{i,P}+{U}_{\sigma }\left(t\right){\Lambda }_{\tau }\left(t\right)\left[{P}_{\tau }\left(t\right){E}_{i,M}-\left(1-{P}_{\tau }\left(t\right)\right){E}_{i,P}\right], \quad \hbox{for}\,\,i=2, 3, 4, 5$$

(6)

$$\frac{d{E}_{2,Q}}{dt}={k}_{Q}({E}_{2,M}+{E}_{2,P})-{k}_{L}{E}_{2,Q}$$

(7)

$$\frac{d{E}_{i,Q}}{dt}={k}_{Q}\left({E}_{i,M}+{E}_{i,P}\right)+{k}_{L}{E}_{i-1,Q} -{k}_{L}{E}_{i,Q}, \quad \hbox{for}\,\,i=3, 4, 5$$

(8)

$$\frac{d{A}_{M}}{dt}={f}_{A}{k}_{L}{E}_{5,M}-{k}_{Q}{A}_{M}-{U}_{\sigma }\left(t\right){\Lambda }_{\tau }(t)[{P}_{\tau }(t){A}_{M}-\left(1-{P}_{\tau }(t)\right){A}_{P}]-{c}_{A}{A}_{M}$$

(9)

$$\frac{d{A}_{P}}{dt}={f}_{A}{k}_{L}{E}_{5,P}-{k}_{Q}{A}_{P}+{U}_{\sigma }\left(t\right){\Lambda }_{\tau }(t)[{P}_{\tau }(t){A}_{M}-\left(1-{P}_{\tau }(t)\right){A}_{P}]-{c}_{A}{A}_{P}$$

(10)

$$\frac{d{A}_{Q}}{dt}={f}_{A}{k}_{L}{E}_{5,Q}+{k}_{Q}\left({A}_{M}+{A}_{P}\right)-{c}_{A}{A}_{Q}$$

(11)

$$\frac{d{I}_{M}}{dt}=\left(1-{f}_{A}\right){k}_{L}{E}_{5,M}-\left({k}_{Q}+{j}_{Q}\right){I}_{M}-{U}_{\sigma }\left(t\right){\Lambda }_{\tau }\left(t\right)\left[{P}_{\tau }\left(t\right){I}_{M}-\left(1-{P}_{\tau }\left(t\right)\right){I}_{P}\right]-{c}_{I}{I}_{M}$$

(12)

$$\frac{d{I}_{P}}{dt}=\left(1-{f}_{A}\right){k}_{L}{E}_{5,P}-\left({k}_{Q}+{j}_{Q}\right){I}_{P}+{U}_{\sigma }\left(t\right){\Lambda }_{\tau }(t)[{P}_{\tau }(t){I}_{M}-\left(1-{P}_{\tau }\left(t\right)\right){I}_{P}]-{c}_{I}{I}_{P}$$

(13)

$$\frac{d{I}_{Q}}{dt}=\left(1-{f}_{A}\right){k}_{L}{E}_{5,Q}+\left({k}_{Q}+{j}_{Q}\right)\left({I}_{M}+{I}_{P}\right)-{c}_{I}{I}_{Q}$$

(14)

$$\frac{d{H}_{U}}{dt}={f}_{H}{c}_{I}\left({I}_{M}+{I}_{P}+{I}_{Q}\right)-{c}_{H}{H}_{U}$$

(15)

$$\frac{dD}{dt}=\left(1-{f}_{R}\right){c}_{H}{H}_{U}+\left(1-{f}_{R}\right){c}_{H}{H}_{V}$$

(16)

$$\frac{d{R}_{U}}{dt}={c}_{A}\left({A}_{M}+{A}_{P}+{A}_{Q}\right)+\left(1-{f}_{H}\right){c}_{I}\left({I}_{M}+{I}_{P}+{I}_{Q}\right)+{f}_{R}{c}_{H}{H}_{U}-\mu \left(t\right){S}_{0}\left(\frac{{R}_{U}}{{\phi }_{V}(t)}\right)$$

(17)

$$\frac{d{V}_{1} }{dt}=\mu \left(t\right){S}_{0}\left(\frac{{S}_{M}+{S}_{P}}{{\phi }_{V}(t)}\right)-{k}_{V}{V}_{1}-\beta {Y}_{\theta }\left(t\right)\left(\frac{{V}_{1}}{{S}_{0}}\right)\left({\phi }_{M}\left(t,\rho \right)+{m}_{b}{\phi }_{P}\left(t,\rho \right)\right)$$

(18)

$$\begin{aligned} \frac{d{V}_{i}}{dt}&={k}_{V}{V}_{i-1}-{k}_{V}{V}_{i}-\beta {Y}_{\theta }\left(t\right)\left(\frac{{V}_{i}}{{S}_{0}}\right)\left({\phi }_{M}\left(t,\rho \right)+{m}_{b}{\phi }_{P}\left(t,\rho \right)\right),\\ & \quad \hbox{for}\,\,i=2, 3, 4, 5, 6 \end{aligned}$$

(19)

$$\frac{d{S}_{V,1}}{dt}=\left(1-{f}_{0}\right){{k}_{V}V}_{6}-\beta {Y}_{\theta }\left(t\right)\left(\frac{{S}_{V,1}}{{S}_{0}}\right)\left({\phi }_{M}\left(t,\rho \right)+{m}_{b}{\phi }_{P}\left(t,\rho \right)\right)$$

(20)

$$\frac{d{S}_{V,2}}{dt}=\left({f}_{0}-{f}_{1}\right){k}_{V}{V}_{6}-{U}_{{\theta }_{1}}\left(t\right)\beta {Y}_{\theta }\left(t\right)\left(\frac{{S}_{V,2}}{{S}_{0}}\right)\left({\phi }_{M}\left(t,\rho \right)+{m}_{b}{\phi }_{P}\left(t,\rho \right)\right)$$

(21)

$$\frac{d{S}_{V,3}}{dt}=\left({f}_{1}-{f}_{2}\right){k}_{V}{V}_{6}-{U}_{{\theta }_{2}}(t)\beta {Y}_{\theta }\left(t\right)\left(\frac{{S}_{V,3}}{{S}_{0}}\right)\left({\phi }_{M}\left(t,\rho \right)+{m}_{b}{\phi }_{P}\left(t,\rho \right)\right)$$

(22)

$$\frac{d{S}_{V,4}}{dt}={f}_{2}{{k}_{V}V}_{6}$$

(23)

$$\begin{aligned} \frac{d{E}_{V}}{dt}&=\beta {Y}_{\theta }\left(t\right)\left(\frac{1}{{S}_{0}}\right)\left({\phi }_{M}\left(t,\rho \right)+{m}_{b}{\phi }_{P}\left(t,\rho \right)\right)\\ &\quad \cdot \sum_{i=1}^{6}\left({V}_{i}+{S}_{V,1}+{U}_{{\theta }_{1}}(t){S}_{V,2}+{U}_{{\theta }_{2}}{\left(t\right)S}_{V,3}\right)-\left(\frac{{k}_{L}}{5}\right){E}_{V} \end{aligned}$$

(24)

$$\frac{d{A}_{V}}{dt}={f}_{A}\left(\frac{{k}_{L}}{5}\right){E}_{V}-{c}_{A}{A}_{V}$$

(25)

$$\frac{d{I}_{V}}{dt}=\left(1-{f}_{A}\right)\left(\frac{{k}_{L}}{5}\right){E}_{V}-{c}_{I}{I}_{V}$$

(26)

$$\frac{d{H}_{V}}{dt}={m}_{h}{f}_{H}{c}_{I}{I}_{V}-{c}_{H}{H}_{V}$$

(27)

$$\frac{d{R}_{V}}{dt}=\mu \left(t\right){S}_{0}\left(\frac{{R}_{U}}{{\phi }_{V}(t)}\right)+{c}_{A}{A}_{V}+\left(1-{m}_{h}{f}_{H}\right){c}_{I}{I}_{V}+{f}_{R}{c}_{H}{H}_{V}$$

(28)

Fig. 7

Expanded illustration of the new compartmental model. In the extended model, vaccination was considered by allowing susceptible and recovered persons to transition into a vaccinated compartment, either ${V}_{1}$ and ${R}_{V}$. Susceptible (recovered) persons who have completed vaccination move into the ${V}_{1}$ (${R}_{V}$) compartment. The susceptible persons who move into ${V}_{1}$ are drawn from ${S}_{M}$ (populated by susceptible persons who are mixing and unprotected by social-distancing) and from ${S}_{P}$ (populated by susceptible persons who are protected by social-distancing). After susceptible persons enter ${V}_{1}$, they can move through a series of additional compartments (${V}_{2}$ through ${V}_{6}$), which are included to capture the time needed for immunity to develop after completion of vaccination. We estimate that the time needed to acquire immunity after vaccination is approximately three weeks based on longitudinal studies of anti-spike protein IgG levels (Bo et al. 2021). Persons who exit the ${V}_{6}$ compartment without becoming infected enter one of the following compartments: ${S}_{V,1}$, ${S}_{V,2}$, ${S}_{V,3}$, or ${S}_{V,4}$. Persons in ${S}_{V,1}$ are taken to remain susceptible to productive infection by all SARS-CoV-2 strains of interest (Alpha, Delta, and ancestral strains). Persons in ${S}_{V,2}$ are taken to be susceptible to SARS-CoV-2 Alpha and Delta variants. Persons in ${S}_{V,3}$ are taken to be susceptible to Delta. Persons in ${S}_{V,4}$ are taken to be protected against all strains of interest. Infection of persons in ${S}_{V,3}$ is only allowed if Delta is present, i.e., at times $t>{\theta }_{2}$. Infection of persons in ${S}_{V,2}$ is only allowed if Alpha or Delta is present, i.e., at times $t>{\theta }_{1}$. Vaccinated persons in compartments ${V}_{1}$ through ${V}_{6}$ and compartment ${S}_{V,1}$ are allowed to become infected at any time, at which point they transition to compartment ${E}_{V}$, consisting of vaccinated persons who were exposed before development of vaccine-induced immunity. Persons in compartment ${S}_{V,2}$ are allowed to become infected if $t\ge {\theta }_{1}$. Similarly, persons in compartment ${S}_{V,3}$ are allowed to become infected if $t\ge {\theta }_{2}$. Possible outcomes for persons in ${E}_{V}$ are taken to be the same as those for unvaccinated exposed persons; however, the incubation period is taken to be distinct. Persons in ${E}_{V}$ can experience asymptomatic disease (upon entering ${A}_{V})$ or they can become symptomatic (upon entering ${I}_{V}$). Persons in ${A}_{V}$ eventually recover, entering compartment ${R}_{V}$. Persons in ${I}_{V}$ can progress to severe disease (upon entering ${H}_{V}$) or recover (upon entering ${R}_{V}$). Persons in ${H}_{V}$ either recover (moving into ${R}_{V}$) or die (moving into $D$). Persons who have recovered from infection, in the ${R}_{U}$ compartment, move directly into the ${R}_{V}$ compartment upon vaccination. Persons in the ${R}_{U}$ and ${R}_{V}$ compartments are taken to have full immunity. The vaccination rate at which susceptible and recovered persons move into vaccinated compartments is updated daily for consistency with the empirical overall rate of vaccination, which we extract daily from the COVID Act Now Data API (Courtemanche et al. 2020) and the Democrat and Chronicle newspaper (Hsiang et al. 2020). The relative values of the vaccination rate are set such that each person eligible for vaccination has the same probability of being vaccinated. All unvaccinated persons are taken to be eligible for vaccination except symptomatic persons (in compartments ${I}_{M}$ and ${I}_{P}$), persons who are hospitalized or severely ill at home (in compartment $H$), quarantined persons (in the various compartments labeled with a Q subscript), and deceased persons (in compartment $D$). It should be noted that asymptomatic, non-quarantined persons (in compartments $A{}_{M}$ and ${A}_{P}$) and presymptomatic, non-quarantined persons (in the $E$ compartments) are taken to be eligible (and to influence the vaccination rate constants) but, as a simplification, these persons are not explicitly tracked as vaccinated or unvaccinated because each of these persons will eventually enter either the $D$ compartment or the ${R}_{U}$ compartment, at which point they will have immunity. In the model, the effects of SARS-CoV-2 variants are captured by a time-dependent dimensionless multiplier ${Y}_{\theta }(t)$ of the rate constant $\beta$. This rate constant, which appears in Appendix Eqs. 1–4, 18–22, and 24, determines the rate of disease transmission within the subpopulation unprotected by social-distancing behaviors when ${Y}_{\theta }\left(t\right)=1$. We take ${Y}_{\theta }\left(t\right)=1$ for times $t<{\theta }_{1}$, i.e., for the initial period of the COVID-19 pandemic in the US that we take to have started on January 21, 2020. We take ${Y}_{\theta }\left(t\right)$ to have the form of a step function with distinct values greater than 1 for periods starting at $t={\theta }_{1}$ and ${\theta }_{k+1}>{\theta }_{k}$ for $k=1,\dots ,m.$ Thus, the model allows for $m$ distinct periods of variant strain dominance delimited by a set of start times $\theta =\{{\theta }_{1},\dots ,{\theta }_{m}\}$. We considered $m=2$. We assume that variants differ only in transmissibility (Color Figure Online)

In these equations, the independent variable is time $t$, and the state variables (${S}_{M}$, ${S}_{P}$, ${E}_{1,M},\dots ,{E}_{5,M}$, ${E}_{1,P},\dots ,{E}_{5,P}$, ${E}_{2,Q},\dots ,{E}_{5,Q}$, ${A}_{M}$, ${A}_{P}$, ${A}_{Q}$, ${I}_{M}$, ${I}_{P}$, ${I}_{Q}, {H}_{U}$, $D$, ${R}_{U}$, ${V}_{1},\dots , {V}_{6}$, ${S}_{V,1},\dots , {S}_{V,4}$, ${E}_{V}$, ${A}_{V}$, ${I}_{V}$, ${H}_{V}$, and ${R}_{V}$) represent 40 (sub)populations (Table 5 in Appendix), which change over time. Thus, each ODE in Eqs. (1)–(28) defines the time-rate of change of a population, i.e., the time-rate of change of a state variable. Note that Eqs. (5), (6), (8) and (19) define 4, 4, 3, and 5 ODEs of the model, respectively. The model is formulated such that ${S}_{0}$, the total population, is a constant. Thus, the model does not account for birth, death for reasons other than COVID-19, immigration, or emigration.

The initial condition associated with Eqs. (1)–(28) is taken to be ${S}_{M}\left({t}_{0}\right)={S}_{0}$, ${I}_{M}\left({t}_{0}\right)={I}_{0}=1$, and all other populations (${S}_{P}$, ${E}_{1,M},\dots ,{E}_{5,M}$, ${E}_{1,P},\dots ,{E}_{5,P}$, ${E}_{2,Q},\dots ,{E}_{5,Q}$, ${A}_{M}$, ${A}_{P}$, ${A}_{Q}$, ${I}_{P}$, ${I}_{Q}, {H}_{U}$, $D$, ${R}_{U}$, ${V}_{1},\dots , {V}_{6}$, ${S}_{V,1},\dots , {S}_{V,4}$, ${E}_{V}$, ${A}_{V}$, ${I}_{V}$, ${H}_{V}$, and ${R}_{V}$) are equal to 0. Recall that the parameter ${S}_{0}$ denotes the total region-specific population size. Thus, we assume that the entire population is susceptible at the start of the local epidemic at time $t={t}_{0}>0$, where time $t=0$ corresponds to 00:00 h (midnight) on January 21, 2020. The parameter ${I}_{0}$ denotes the number of infectious symptomatic persons at the start of the regional epidemic.

In the model, the parameters $\beta$, ${k}_{L}$, ${k}_{Q}$, ${j}_{Q}$, ${c}_{A}$, ${c}_{I}$, ${c}_{H}$, and ${k}_{V}$ are positive-valued rate constants (all with units of ${\text{d}}^{-1}$), and the parameters ${m}_{b}$, ${m}_{h}$, ${f}_{A}$, ${f}_{H}$, ${f}_{R}$, ${f}_{0}\ge {f}_{1}$, ${f}_{1}\ge {f}_{2}$, and ${f}_{2}$ are (dimensionless) fractions. Brief definitions of parameters are given in Tables 1 and 2.

In the model, the quantities ${\phi }_{M}(t,\rho )$, ${\phi }_{P}(t,\rho )$, and ${\phi }_{V}\left(t\right)$ are functions of (time-dependent) state variables (as defined below), which represent the population of infectious persons who are mixing freely (i.e., not practicing social-distancing), the population of infectious persons who are practicing social-distancing (i.e., adopting disease-avoiding behaviors), and the population of persons eligible for vaccination, respectively. The quantities ${\phi }_{M}(t,\rho )$ and ${\phi }_{P}(t,\rho )$ are also functions of $\rho \equiv ({\rho }_{E},{\rho }_{A})$, where ${\rho }_{E}$ (${\rho }_{A}$) is a dimensionless ratio representing the infectiousness of persons in the incubation phase of infection (the infectiousness of asymptomatic persons in the immune clearance phase of infection) relative to the infectiousness of symptomatic persons with the same social-distancing behavior. The quantity ${\phi }_{V}(t)$ represents the population of persons eligible for vaccination.

$${\phi }_{M}\left(t,\rho \right)={I}_{M}+{I}_{V}+{\rho }_{E}\left({E}_{2,M}+{E}_{3,M}+{E}_{4,M}+{E}_{5,M}+{E}_{V}\right)+{\rho }_{A}{(A}_{M}+{A}_{V})$$

(29)

$${\phi }_{P}\left(t,\rho \right)={I}_{P}+{\rho }_{E}\left({E}_{2,P}+{E}_{3,P}+{E}_{4,P}+{E}_{5,P}\right)+{\rho }_{A}{A}_{P}$$

(30)

$${\phi }_{V}\left(t\right)={{S}_{M}+S}_{P}+\sum_{i=1}^{5}\left({E}_{i,M}+{E}_{i,P}\right)+{A}_{M}+{A}_{P}+{R}_{U}$$

(31)

The state variables that appear in these equations represent time-varying populations. Recall that state variables are defined in Table 5 in Appendix.

Table 5 State variables of the compartmental model

In the model, the quantities ${U}_{\sigma }(t)$, ${U}_{{\theta }_{1}}(t)$, and ${U}_{{\theta }_{2}}(t)$ are unit step functions. The values of these functions change from 0 to 1 at the times indicated by the subscripts: $\sigma$, the onset time of the initial social-distancing period; ${\theta }_{1}$, the takeoff time of SARS-CoV-2 variant Alpha; and ${\theta }_{2}$, the takeoff time of SARS-CoV-2 variant Delta.

$${U}_{\sigma }\left(t\right)=\left\{ \begin{array}{cc}0& t<\sigma \\ 1& t\ge \sigma \end{array}\right.$$

(32)

$${U}_{{\theta }_{1}}\left(t\right)=\left\{ \begin{array}{cc}0& t<{\theta }_{1}\\ 1& t\ge {\theta }_{1}\end{array}\right.$$

(33)

$${U}_{{\theta }_{2}}\left(t\right)=\left\{ \begin{array}{cc}0& t<{\theta }_{2}\\ 1& t\ge {\theta }_{2}\end{array}\right.$$

(34)

As indicated in Fig. 7 in Appendix, transitions from ${S}_{M}$ to ${S}_{P}$, for example, become possible at time $t=\sigma$, transitions from ${S}_{V,2}$ to ${E}_{V}$ become possible at time $t={\theta }_{1}$, and transitions from ${S}_{V,3}$ to ${E}_{V}$ become possible at time $t={\theta }_{2}$.

In the model, the quantities ${P}_{\tau }\left(t\right)$, and ${\Lambda }_{\tau }(t)$ are step functions that characterize changes in social-distancing. The value of ${P}_{\tau }(t)$ determines a setpoint steady-state fraction of susceptible persons who are practicing social-distancing. The value of ${\Lambda }_{\tau }\left(t\right)$ determines a time scale for approach to the setpoint steady state. Changes in the values of ${P}_{\tau }(t)$ and ${\Lambda }_{\tau }(t)$ occur coordinately. These changes occur at times $t=\sigma ,{\tau }_{1},\dots , {\tau }_{n}$, where $n$ is the number of distinct social-distancing periods beyond an initial social-distancing period. Initially, we took $n=7$ (i.e., 8 total social-distancing stages). The value of $n$ is decremented by 1 (at an inferred time) if $n\leftarrow n-1$ is determined to be admissible by a model-selection procedure, which is described below. It should be noted that ${{p}_{0}, p}_{1},\dots ,{p}_{n}$ are parameters of ${P}_{\tau }(t)$ and that ${{\lambda }_{0}, \lambda }_{1},\dots ,{\lambda }_{n}$ are parameters of ${\Lambda }_{\tau }(t)$. These parameters determine the non-zero values of the step functions over different periods. For example, ${p}_{1}$ is the value of ${P}_{\tau }(t)$ and ${\lambda }_{1}$ is the value of ${\Lambda }_{\tau }(t)$ for the period $t\in [{\tau }_{1},{\tau }_{2})$.

$${P}_{\tau }\left(t\right)=\left\{\begin{array}{c}\begin{array}{cc}{p}_{0}& \sigma \le t<{\tau }_{1}\\ {p}_{1}& {\tau }_{1}\le t<{\tau }_{2}\end{array}\\ \begin{array}{cc}\vdots & \vdots \\ {p}_{n}& {\tau }_{n}\le t<\infty \end{array}\end{array}\right.$$

(35)

$${\Lambda }_{\tau }\left(t\right)=\left\{\begin{array}{c}\begin{array}{cc}{\lambda }_{0}& \sigma \le t<{\tau }_{1}\\ {\lambda }_{1}& {\tau }_{1}\le t<{\tau }_{2}\end{array}\\ \begin{array}{cc}\vdots & \vdots \\ {\lambda }_{n}& {\tau }_{n}\le t<\infty \end{array}\end{array}\right.$$

(36)

The values of ${P}_{\tau }(t)$ and ${\Lambda }_{\tau }(t)$ are 0 for $t<\sigma$.

In the model, the quantity $\mu (t)$ is a piecewise linear interpolant to a function $\widetilde{\mu }\left(t\right)$ that characterizes the current rate of vaccination. The value of $\widetilde{\mu }\left(t\right)$ is determined by the empirical daily rate of vaccination, and thus, can vary from day to day. Daily vaccination data were extracted from the Covid Act Now database and the Democrat and Chronicle newspaper (footnote 8) using the Covid Act Now Data API (footnote 7) and web scraping. We will use ${\mu }_{i}$ to refer to the value of $\widetilde{\mu }\left(t\right)$ for $t\in \left[{t}_{i},{t}_{i+1}\right)$, where time ${t}_{i}$ corresponds to midnight on the $i$th day after January 21, 2020.

$$\widetilde{\mu }\left(t\right)={\mu }_{i}\text{ for }t\in [{t}_{i},{t}_{i+1})$$

(37)

Settings for ${\mu }_{i}$ were made such that ${\mu }_{i}{S}_{0}\times 1 {\text{d}}$ is the number of vaccinations completed in the nearest past 1-d period according to Covid Act Now data and the Democrat and Chronicle newspaper.

In the model, the quantity ${Y}_{\theta }(t)$ is a step function that quantifies how disease transmissibility increases upon emergence of SARS-CoV-2 variants Alpha and Delta. Initially, ${Y}_{\theta }\left(t\right)=1$. The value of ${Y}_{\theta }(t)$ is increased (by an inferred factor greater than 1 at an inferred time, ${\theta }_{1}$ or ${\theta }_{2}$) if the change is determined to be admissible by a model-selection procedure, which is described below. It should be noted that ${y}_{1}$ and ${y}_{2}$ are parameters of ${Y}_{\theta }(t)$. These parameters determine the values of the step function ${Y}_{\theta }(t)$ over different periods: ${y}_{1}$ is the value of ${Y}_{\theta }(t)$ for the period $t\in [{\theta }_{1},{\theta }_{2})$ and ${y}_{2}$ is the value of ${Y}_{\theta }(t)$ for the remaining period of concern (with Delta as the dominant circulating viral strain).

$${Y}_{\theta }\left(t\right)=\left\{\begin{array}{c}\begin{array}{cc}{y}_{0}& {\theta }_{0}\le t<{\theta }_{1}\\ {y}_{1}& {\theta }_{1}\le t<{\theta }_{2}\end{array}\\ \begin{array}{cc}\vdots & \vdots \\ {y}_{m}& {\theta }_{m}\le t<\infty \end{array}\end{array}\right.$$

(38)

where $m$ is the number of viral variants that have emerged up to the current time, ${\theta }_{0}\equiv {t}_{0}$, and ${y}_{0}\equiv 1$. Here, we consider $m= 2$.

1.3 Equations of the Auxiliary Measurement Model

As in the study of Lin et al. (2021), we assumed that only symptomatic persons are detected in testing. The accumulation of symptomatic persons is governed by

$$\frac{d{C}_{S}}{dt}=\left(1-{f}_{A}\right)\left[{k}_{L}\left({E}_{5,M}+{E}_{5,P}+{E}_{5,Q}\right)+{\left(\frac{{k}_{L}}{5}\right)E_V}\right]$$

(39)

where ${C}_{S}(t)$ is the cumulative number of symptomatic persons (cases) at time $t$. Here, unlike in the study of Lin et al. (2021), the expression for ${C}_{S}(t)$ accounts for exposed persons in quarantine. Initially, ${C}_{S}=0.$ We numerically integrated Appendix Eq. (39) together with the ODEs of the compartmental model. From the trajectory for ${C}_{S}$, we derive a prediction for the expected number of new COVID-19 cases reported on calendar date ${\mathcal{D}}_{{\text{i}}}$, $I\left({t}_{i},{t}_{i+1}\right)$, using the following equation:

$$I\left({t}_{i},{t}_{i+1}\right)={f}_{D}\left[{C}_{S}\left({t}_{i+1}\right)-{C}_{S}\left({t}_{i}\right)\right]$$

(40)

where ${f}_{D}$ is an adjustable region-specific parameter characterizing the time-averaged fraction of symptomatic cases detected and reported, ${t}_{i}$ corresponds to midnight on the $i$th day after January 21, 2020, and ${t}_{i+1}-{t}_{i}$ is the reporting interval (1 d). We compare $I({t}_{i},{t}_{i+1})$ to $\delta {C}_{i}$, the number of new cases reported on calendar date ${\mathcal{D}}_{{\text{i}}}$.

1.4 Definition of the Likelihood Function

Bayesian inference relies on the definition of a likelihood, which here serves the purpose of assessing the compatibility of available surveillance data with adjustable (free) parameter values. Let us use ${\left\{\delta {C}_{i}\right\}}_{i=0}^{d}$ to denote the daily case reporting data available between 0 and $d$ days after midnight on January 21, 2020 (the date of the first case report in the US) and let ${D=\left\{\delta {C}_{i}\right\}}_{i=0}^{d}$. Let us use ${\theta }_{F}(n,m)$ to denote the set of adjustable (free) parameter values. The number of adjustable parameters, $|{\theta }_{F}|$, depends on $n$, the number of social-distancing periods considered beyond an initial social-distancing period, and $m$, the number of SARS-CoV-2 variants under consideration. As in the study of Lin et al. (2021), we assume that $\delta {C}_{i}$, the number of new COVID-19 cases detected over a 1-d period and reported on calendar date ${\mathcal{D}}_{{\text{i}}}$ for a given region, is a random variable and its expected value follows a model-derived deterministic trajectory given by $I({t}_{i},{t}_{i+1})$ (Eq. 40). We assume that day-to-day fluctuations in the random variable are independent and characterized by a negative binomial distribution ${\text{NB}}\left(r,{q}_{i}\right)$, which has two parameters, $r>0$ and ${q}_{i}\in (\mathrm{0,1})$. Note that ${\mathbb{E}}\left[{\text{NB}}(r,{q}_{i})\right]=r(1-{q}_{i})/{q}_{i}$. We assume that this distribution has the same dispersion parameter $r$ across all case reports. With these assumptions, we arrive at the following likelihood function:

$$\mathcal{L}\left({\theta }_{F}(n,m);{\left\{\delta {C}_{i}\right\}}_{i=0}^{d}\right)=\prod_{i=0}^{d}{\mathcal{L}}_{i}({\theta }_{F}(n,m);\delta {C}_{i})$$

(41)

where

$${\mathcal{L}}_{i}\left({\theta }_{F}(n,m);\delta {C}_{i}\right)={\text{nbinom}}\left(\delta {C}_{i};r,{q}_{i}\right)=\left(\genfrac{}{}{}{}{\delta {C}_{i}+r-1}{\delta {C}_{i}-1}\right){q}_{i}^{r}{\left(1-{q}_{i}\right)}^{\delta {C}_{i}}$$

(42)

and

$${q}_{i}=\frac{r}{r+I({t}_{i},{t}_{i+1})}.$$

(43)

In these equations, $i$ is an integer indicating the date ${\mathcal{D}}_{i}$ and period $\left({t}_{i},{t}_{i+1}\right)$; ${\text{nbinom}}\left(\delta {C}_{i};r,{q}_{i}\right)$ is the probability mass function of the negative binomial distribution ${\text{NB}}(r,{q}_{i})$, and ${\theta }_{F}(n,m)=\{{t}_{0},\beta , \sigma , {\tau }_{1},\dots ,{\tau }_{n},{p}_{0},{p}_{1,}\dots ,{p}_{n},{\lambda }_{0},{\lambda }_{1},\dots ,{\lambda }_{n},{\theta }_{1},\dots ,{\theta }_{m}, {y}_{1},\dots ,{y}_{m},{f}_{D}, r\}$ for $n\ge 1$ and $m\ge 1$; ${\theta }_{F}(\mathrm{0,0})=\{{t}_{0},\beta ,{p}_{0},{\lambda }_{0},{f}_{D}, r\}$.

1.5 Parameters

Each model parameter is briefly described in Tables 1, 2, and 3. These parameters have either fixed values or adjustable values (i.e., values inferred from surveillance data). The fixed values may be universal (i.e., applicable to all MSAs of interest) or MSA-specific. All inferred parameter values are MSA-specific. In addition, the measurement model (Appendix Eqs.  39 and 40) has one adjustable MSA-specific parameter, ${f}_{D}$, and the likelihood function (Appendix Eqs.  41–43) has one adjustable MSA-specific parameter, $r$. Values of the other likelihood parameters, ${q}_{0},\dots ,{q}_{d}$, are constrained and are determined using Appendix Equation 43.

1.6 Original Model of Lin et al.

The model shares $19+3n$ parameters with the model of Lin et al. (2021), including parameters that define the initial condition (${t}_{0}$, ${I}_{0}$, and ${S}_{0}$). (Recall that $n$ is the number of social-distancing periods being considered beyond the initial social-distancing period.) The shared parameters are ${t}_{0}$, ${I}_{0}$, ${S}_{0}$, $\beta$, $\sigma$, ${\tau }_{1},\dots , {\tau }_{n}$, ${p}_{0},\dots ,{p}_{n}$, ${\lambda }_{0},\dots ,{\lambda }_{n}$, ${\rho }_{A}$, ${\rho }_{E}$, ${m}_{b}$, ${f}_{A}$, ${f}_{H}$, ${f}_{R}$, ${k}_{L}$, ${k}_{Q}$, ${j}_{Q}$, ${c}_{A}$, ${c}_{H}$, and ${c}_{I}$. As in the study of Lin et al. (2021), we inferred MSA-specific values for the following parameters: ${t}_{0}$, $\beta$, $\sigma$, ${p}_{0},\dots , {p}_{n}$, and ${\lambda }_{0},\dots , {\lambda }_{n}$. We also inferred MSA-specific values for ${\tau }_{1},\dots ,{\tau }_{n}$ provided that $n\ge 1$. As in the study of Lin et al. (2021), the remaining 14 parameters shared between the old and new models (${I}_{0}$, ${S}_{0}$, ${\rho }_{A}$, ${\rho }_{E}$, ${m}_{b}$, ${f}_{A}$, ${f}_{H}$, ${f}_{R}$, ${k}_{L}$, ${k}_{Q}$, ${j}_{Q}$, ${c}_{A}$, ${c}_{H}$, and ${c}_{I}$) were taken to have fixed values, and we adopted the settings of Lin et al. (2021) for these parameters (Table 3). These settings are universal except for the setting for ${S}_{0}$, the total population, which is MSA-specific.

1.7 Extension of the Model of Lin et al.

Our extension of the model of Lin et al. (2021) introduces $5+2\left(m+1\right)+(d+1)$ parameters, where $m$ ($=0$, 1 or 2) is the number of SARS-CoV-2 variants being considered and $d$ is the number of days since January 21, 2020: ${\theta }_{0}, \dots ,{\theta }_{m}$, ${y}_{0},\dots ,{y}_{m}$, ${m}_{h}$, ${f}_{0}$, ${f}_{1}$, ${f}_{2}$, ${k}_{V}$, and ${\mu }_{0},\dots ,{\mu }_{d}$. The $\theta$ and $y$ parameters are variant takeover times and transmissibility factors, respectively, except that the value of ${\theta }_{0}$ is defined as ${t}_{0}$ and the value of ${y}_{0}$ is defined as 1. The Alpha transmissibility factor ${y}_{1}$, the Alpha takeoff time ${\theta }_{1}$, the Delta transmissibility factor ${y}_{2}$, and the Delta takeoff time ${\theta }_{2}$ were inferred for each MSA with $m=2$ (cf. Table 1). The transmissibility factors were each constrained to be greater than or equal to 1. The settings for ${\mu }_{0},\dots ,{\mu }_{d}$ are empirical and MSA-specific. Each ${\mu }_{i}$ is set such that ${\mu }_{i}{S}_{0}\times 1 {\text{d}}$ is the number of vaccinations completed over the past 1-d period nearest to the $i$th day after January 21, 2020. As noted earlier, the number of completed vaccinations was obtained for each MSA from Covid Act Now and the Democrat and Chronicle newspaper (footnote 8) using the Covid Act Now Data API (footnote 7) and web scraping. In the spreadsheet accessed daily from Covid Act Now, the ‘metrics.vaccinationsCompletedRatio’ column gives the percentage of the total population that has received the recommended number of doses: one dose for Ad26.CoV2.S (Janssen, Johnson & Johnson) or two doses for mRNA-1273 (Moderna) and BNT162b2 (Pfizer-BioNTech). As a simplification, we considered all completed vaccinations to be equivalent. The parameters ${m}_{h}$, ${f}_{0}$, ${f}_{1}$, ${f}_{2}$, and ${k}_{V}$ were assigned fixed universal estimates (Table 3). Each of these estimates is explained below.

1.8 Estimation of ${{{k}}}_{{{V}}}$

The rate constant ${k}_{V}$ characterizes the rate of transition out of compartment ${V}_{i}$ for $i=1,\dots ,{n}_{V}$. Recall that, in the model, susceptible persons enter ${V}_{1}$ upon vaccination (Fig. 1, Fig. 7 in Appendix). The values of ${n}_{V}$ ($=6$) and ${k}_{V}$ ($=0.3 \, {\text{d}}^{-1}$) were selected so that the time a person spends in ${V}_{1},\dots ,{V}_{{n}_{V}}$, which we will denote as ${t}_{V}$, is distributed approximately the same as ${\widetilde{t}}_{V}$, the waiting time between vaccination of a previously uninfected person and detection of vaccine-induced SARS-CoV-2-specific IgG antibodies (Korodi et al. 2022) (Fig. 8 in Appendix). According to the model, the time that a person spends in ${V}_{1},\dots ,{V}_{{n}_{V}}$ is distributed according to the probability density function $f\left({t}_{V};{n}_{V},{k}_{V}\right)={k}_{V}^{{n}_{V}}{t}_{V}^{{n}_{V}-1}{e}^{-{k}_{V}{t}_{V}}/\left({n}_{V}-1\right)!$, i.e., ${t}_{V}$ is Erlang distributed with shape parameter ${n}_{V}=6$ and rate parameter ${k}_{V}=0.3$ ${\text{d}}^{-1}$. As can be seen in Fig. 8 in Appendix, the cumulative distribution function of this Erlang distribution reasonably captures the empirical cumulative distribution of waiting times observed in the longitudinal study of Korodi et al. (2022). Thus, in the model, passage through ${V}_{1},\dots ,{V}_{6}$ with rate constant ${k}_{V}=0.3 \, {\text{d}}^{-1}$ accounts for the variable and significantly non-zero amount of time required for development of a protective antibody response after vaccination.

Fig. 8

Comparison of an Erlang cumulative distribution function with shape parameter ${n}_{V}=6$ and rate parameter ${k}_{V}=0.3 {\text{d}}^{-1}$ and the empirical cumulative distribution of waiting times (${\widetilde{t}}_{V}$ values) observed in the longitudinal study of Korodi et al. (Bo et al. 2021). The waiting time ${\widetilde{t}}_{V}$ is the time between vaccination of a previously uninfected person and detection of vaccine-induced SARS-CoV-2-specific IgG antibodies (color figure online)

1.9 Estimation of ${{{f}}}_{0}, {{{f}}}_{1},$ and ${{{f}}}_{2}$

The parameters ${f}_{0}>{f}_{1}$, ${f}_{1}>{f}_{2}$, and ${f}_{2}$ are fractions that characterize the average effectiveness of vaccines used in the US and that determine the sizes of (mutually exclusive) subpopulations of vaccinated persons having different susceptibilities to productive infection (i.e., an infection that can be transmitted to others): ${S}_{V,1}$, ${S}_{V,2}$, ${S}_{V,3}$, and ${S}_{V,4}$ (Fig. 1, Fig. 7 in Appendix). We assume that persons in the ${S}_{V,1}$ subpopulation are susceptible to productive infection by any of the viral strains considered, and in contrast, we assume that persons in the ${S}_{V,4}$ subpopulation are susceptible to productive infection by none of the viral strains considered. Persons in the ${S}_{V,2}$ subpopulation are taken to be susceptible to productive infection by the Alpha and Delta variants but not viral strains in circulation before the emergence of Alpha. Persons in the ${S}_{V,3}$ subpopulation are taken to be susceptible to productive infection by the Delta variant but not the Alpha variant or viral strains in circulation before the emergence of Alpha. The quantity $1-{f}_{0}$ defines the fraction of vaccinated persons who enter the ${S}_{V,1}$ subpopulation after exiting ${V}_{6}$, the quantity ${f}_{0}-{f}_{1}$ defines the fraction of vaccinated persons who enter the ${S}_{V,2}$ subpopulation after exiting ${V}_{6}$, the quantity ${f}_{1}-{f}_{2}$ defines the fraction of vaccinated persons who enter the ${S}_{V,3}$ subpopulation after exiting ${V}_{6}$, and ${f}_{2}$ defines the fraction of vaccinated persons who enter the ${S}_{V,4}$ subpopulation after exiting ${V}_{6}$. We take ${f}_{0}$ to characterize vaccine effectiveness before the emergence of Alpha. According to Thompson et al. (Thompson et al. 2021), vaccine effectiveness was initially 90%. Thus, we set ${f}_{0}=0.9$. We take ${f}_{1}$ to characterize vaccine effectiveness after the emergence of Alpha but before the emergence of Delta. According to Puranik et al. (2024), in May 2021, vaccine effectiveness was 81%. Thus, we set ${f}_{1}=0.81$. We take ${f}_{2}$ to characterize vaccine effectiveness after the emergence of Delta. According to Tang et al. (2021), the effectiveness of two doses of the Pfizer-BioNTech vaccine (BNT162b2) against Delta is 53.5% and the effectiveness of two doses of the Moderna vaccine (mRNA-1273) against Delta is 84.8%. Taking the average of these figures, we set ${f}_{2}=0.69$.

1.10 Estimation of ${{{m}}}_{{{h}}}$

The parameter ${m}_{h}$ characterizes the reduced risk of severe disease for a vaccinated person in the case of a breakthrough infection. We set ${m}_{h}=0.04$, i.e., we assumed that there is a 25-fold reduction in the risk of severe disease for infected persons who have been vaccinated, which is consistent with the observations of Lopez Bernal et al. (2021).

1.11 Notable New Modeling Assumptions

It should be noted that we treat the incubation period for newly infected (exposed) vaccinated persons differently than for newly infected (exposed) unvaccinated persons (Fig. 1, Fig. 7 in Appendix). For unvaccinated persons, as in the study of Lin et al. (2021), we divide exposed persons in the incubation period of infection into five subpopulations: ${E}_{1,M},\dots ,{E}_{5,M}$ for exposed persons who are mixing (i.e., persons who are not practicing social-distancing), ${E}_{1,P},\dots , {E}_{5,P}$ for exposed persons who are practicing social-distancing, and ${E}_{1,Q},\dots , {E}_{5,Q}$ for exposed quarantined persons. Persons move through the five stages of the incubation period sequentially. In contrast, as a simplification, for vaccinated persons, we consider only a single exposed population: ${E}_{V}$. We take persons to exit ${E}_{V}$ with rate constant ${k}_{L}/5$ (Fig. 7 in Appendix). With this choice, the duration of the incubation period of infection is the same, on average, for both vaccinated and unvaccinated persons. The average duration is $5/{k}_{L}$ (about 5 d) in both cases. The difference is that the duration of the incubation period is Erlang distributed for unvaccinated persons, as discussed by Lin et al. (2021), but exponentially distributed for vaccinated persons.

As indicated in Eq. (29), we take vaccinated persons with productive infections to be equally as infectious as unvaccinated persons.

As noted earlier, we take all vaccinated persons to be mixing (i.e., to not be practicing social-distancing). Thus, populations of infected vaccinated persons (${E}_{V}$, ${I}_{V}$, and ${A}_{V}$) contribute to ${\phi }_{M}(t)$ (Appendix Eq. 29) but not ${\phi }_{P}\left(t\right)$ (Appendix Eq. 30).

As indicated in Appendix Eq. 31, we consider pre-symptomatic exposed and asymptomatic unvaccinated persons to be eligible for vaccination and, thus, these persons contribute to the consumption of vaccine doses (i.e., these persons account for a portion of the number of completed vaccinations on a given day $i$, ${\mu }_{i}{S}_{0}\times 1 {\text{d}}$). However, we do not move these persons to vaccinated compartments. The reason is that exposed and asymptomatic persons are expected to develop immunity faster through recovery from infection (i.e., movement to ${R}_{U}$) than from vaccination.

As indicated in Appendix Eq. 31, we do not consider symptomatic, quarantined, severely ill and hospitalized/isolated-at-home, and deceased persons to be eligible for vaccination.

1.12 Inference Approach

Recall that ${\theta }_{F}$ denotes the set of all adjustable parameters. As in the study of Lin et al. (2021), for each MSA, we inferred MSA-specific adjustable parameter values ${\theta }_{F}$ using all MSA-specific surveillance data available up to a specified day of inference ${\mathcal{D}}_{{\text{d}}}$ (i.e., the $d$th day after January 21, 2020). We took a Bayesian inference approach, meaning that, for a given dataset, we generated parameter posterior samples (a collection of ${\theta }_{F}$’s) through Markov chain Monte Carlo (MCMC) sampling. The parameter posterior samples provide a probabilistic characterization of the adjustable parameter values consistent with the dataset used in inference. By drawing samples from the parametric posterior distribution, we generated a posterior predictive distribution for $I({t}_{i},{t}_{i+1})$ for each $i$ of interest. We considered all days from January 21, 2020 to October 31, 2021. In other words, for each $i$ of interest, a prediction for $I({t}_{i},{t}_{i+1})$ was generated for each of many ${\theta }_{F}$’s drawn randomly (with uniform probability) from the parametric posterior distribution. The resulting distribution of $I({t}_{i},{t}_{i+1})$ values is the posterior predictive distribution for $I({t}_{i},{t}_{i+1})$. Recall that $I({t}_{i},{t}_{i+1})$ is given by Appendix Eq. 40 and corresponds to ${\mathbb{E}}\left[\delta {C}_{i}\right]$, the expected number of new COVID-19 cases detected over a 1-d surveillance interval and reported for the $i$th day after January 21, 2020. Observation noise was injected into the posterior predictive distributions by replacing each sampled value for $I({t}_{i},{t}_{i+1})$ with ${X}_{i}\sim {\text{NB}}(r,{q}_{i})$, where $r$ is a member of the sampled set of parameter values ${\theta }_{F}$ used to generate the prediction of $I({t}_{i},{t}_{i+1})$ and ${q}_{i}$ is given by Eq. 43.

According to Bayes’ theorem, given surveillance data $D={\left\{\delta {C}_{i}\right\}}_{i=0}^{d}$, the parametric posterior is given as

$${\mathbb{P}}\left\{{\theta }_{F}|D\right\}=\frac{{\mathbb{P}}\left\{D|{\theta }_{F}\right\}{\mathbb{P}}\left\{{\theta }_{F}\right\}}{Z}$$

(44)

where ${\mathbb{P}}\left\{{\theta }_{F}\right\}$ is the prior (which is formulated to capture knowledge of ${\theta }_{F}$ external to $D$ or to express lack of such knowledge), ${\mathbb{P}}\left\{D|{\theta }_{F}\right\}$ is the likelihood defined by Appendix Eqs. 41–43, and $Z$ is a normalizing constant. We assumed a proper uniform prior, i.e., for each adjustable parameter, we assumed that all values between specified lower and upper bounds are equally likely before consideration of $D$. We used the same bounds as in the study of Lin et al. (2021). Then, we used an adaptive MCMC (aMCMC) algorithm (Andrieu and Thoms 2008) to generate samples from ${\mathbb{P}}\left\{D|{\theta }_{F}\right\}{\mathbb{P}}\left\{{\theta }_{F}\right\}$, which is proportional to ${\mathbb{P}}\left\{{\theta }_{F}|D\right\}$. Thus, the relative probabilities of parameter sets ${\theta }_{F}$ according to ${\mathbb{P}}\left\{{\theta }_{F}|D\right\}$ are correctly represented by the samples.

Specifically, the adaptive MCMC algorithm (Andrieu and Thoms 2008) generates samples from the multivariate parametric posterior for adjustable model parameters (${t}_{0}, \beta$, and parameters for variant emergence and social-distancing), the measurement model parameter ${f}_{D}$, and the likelihood parameter $r$ (Tables 1 and 2). This algorithm is available within the PyBioNetFit software package (Neumann et al. 2022). Use of the algorithm was performed as described by Lin et al. (2021). The report of Neumann et al. (2022) includes helpful general usage advice, which was followed in this study. Inference jobs were executed on a computer cluster. For each inference job, a total of 25 chains were generated, and the chain with the best mixing and convergence properties was selected for subsequent analyses.

Each inference was conditioned on the compartmental model of Fig. 1 (Appendix Eqs. 1–38), settings for the structural parameters $m$ (the number of SARS-CoV-2 variants under consideration) and $n$ (the number of social-distancing periods under consideration beyond an initial social-distancing period), the measurement model (Appendix Eqs. 39 and 40), and fixed parameter estimates (Tables 1 and 2), including empirical daily per capita vaccination rates (i.e., the settings for ${\mu }_{i}$ in Appendix Eq. 37). We assumed a proper uniform prior for each adjustable parameter (Lin et al. 2021) and a negative binomial likelihood function (Appendix Eqs. 41–43). Use of proper uniform priors means that MAP estimates are maximum likelihood estimates (MLEs). In each inference, the data entering the likelihood function, $D={\left\{\delta {C}_{i}\right\}}_{i=0}^{d}$ (Appendix Eq. 41), were MSA-specific daily reports of newly detected COVID-19 cases available up to the date of inference ${\mathcal{D}}_{{\text{d}}}$ (i.e., the $d$-th day after January 21, 2020). Thus, all inferences are region-specific and time-dependent.

1.13 Use of Model Selection to Determine Intervals of Step Functions

Variant takeover times, $\theta =({\theta }_{1},{\theta }_{2})$, and start times of social-distancing periods, $\tau =({\sigma , \tau }_{1},\dots , {\tau }_{n})$, were inferred from data; however, changes of the associated time-dependent step functions, ${Y}_{\theta }\left(t\right)$, ${P}_{\tau }(t)$, and ${\Lambda }_{\tau }(t)$, were introduced only when an increase in model complexity was deemed to be justified. Each decision to introduce variant takeover or start of a new social-distancing period (beyond the initial period) was made using a model-selection procedure, which is described below. It should be noted that ${y}_{1}$ and ${y}_{2}$ are parameters of ${Y}_{\theta }(t)$, ${{p}_{0},p}_{1},\dots ,{p}_{n}$ are parameters of ${P}_{\tau }(t)$, and ${{\lambda }_{0},\lambda }_{1},\dots ,{\lambda }_{n}$ are parameters of ${\Lambda }_{\tau }(t)$. These parameters determine the values of the step functions over different periods. For example, for $n\ge 1$, ${p}_{1}$ is the value of ${P}_{\tau }(t)$ and ${\lambda }_{1}$ is the value of ${\Lambda }_{\tau }(t)$ for the period $t\in [{\tau }_{1},{\tau }_{2})$. Similarly, ${y}_{1}$ is the value of ${Y}_{\theta }(t)$ for the period $t\in [{\theta }_{1},{\theta }_{2})$, and ${y}_{2}$ is the value of ${Y}_{\theta }(t)$ for the period $t\in [{\theta }_{2},{t}_{\text{final}})$, where ${t}_{\text{final}}$ corresponds to the date of inference (October 31, 2021).

We started with a setting of $n=7$ for each MSA of interest (i.e., 8 total social-distancing stages). To determine if $n$ could be reduced, we conducted parsimony checks. In a parsimony check, 100 MLE curves, each constituting a fit to the data, were generated via optimization jobs. In these jobs, the total number of social-distancing stages in the model ($n+1$) was set at 1 less than for the current proposed best fit. Each of the 100 fits was visually inspected to determine whether or not the following criteria held: (1) the quality of fit is acceptable (i.e., comparable to what is obtained with the current proposed best fit); (2) Alpha and Delta surges (identified by sequencing data) are explained at least partly by increased transmissibility; (3) social-distancing setpoint parameter values are feasible (i.e., each setpoint parameter takes on a value between 0.2 and 0.8); and (4) social-distancing changes proximal to an Alpha or Delta surge (if any) occur only after an increase in transmissibility. If one or more of these conditions was not satisfied, we accepted the proposed best fit as the most parsimonious fit to the data. If all conditions held, we updated the proposed best fit, and the parsimony check was repeated for a model with one less social-distancing stage. In the Appendix, we show examples of violations of parsimony-check criteria for the MSA surrounding Houston (Fig. 9 in Appendix).

Fig. 9

Example of a parsimony check with criteria that are violated for the MSA surrounding Houston. A Alpha and Delta surges (identified by sequencing data) are not explained, at least in part, by increased transmissibility, while all other criteria are satisfied. Each broken vertical black line from left to right indicates the date of onset of a social-distancing stage (i.e., with n = 3 for this illustration). Each broken vertical red line from left to right indicates the takeoff date of a variant (namely, Alpha or Delta). B Not all social-distancing setpoint parameter values are feasible, while all other criteria are satisfied. C Social-distancing changes proximal to an Alpha or Delta surge (if any) precede an increase in transmissibility, while all other criteria are satisfied. In Panel A, we see that the orange segment between the two broken vertical red lines, which denotes the relative transmissibility of Alpha, is no larger than the orange segment to its left, which denotes the relative transmissibility of ancestral strains of SARS-CoV-2. In Panel B, the rightmost blue segment, which denotes the social-distancing setpoint parameter value for the final social-distancing stage, has a value of roughly 0.16, which is deemed infeasible by our criterion. Finally, in Panel C, the rightmost broken vertical black line is proximal to the leftmost broken vertical red line, which corresponds to the takeoff time of the Alpha surge, and the social-distancing setpoint parameter value changes prior to the Alpha surge. Note that all three panels show MLE curves, each of which constitutes an acceptable fit to the data (Color Figure Online)

1.14 Simulations

After specification of parameter values (Tables 1, 2, and 3), we used the SciPy (https://www.scipy.org) interface to LSODA (Hindmarsh and Stepleman 1983) to numerically integrate the system of coupled ODEs consisting of the 40 ODEs of the compartmental model and the 1 ODE of the measurement model (Appendix Eqs. 1–39). The initial condition was defined by settings for ${t}_{0}$, ${I}_{0}$, and ${S}_{0}$ (Tables 1, 2, 3). Integration combined with use of Appendix Eq. 40 yielded a prediction of the expected number of new cases detected for each 1-d surveillance period of interest in the past or future: $I({t}_{i},{t}_{i+1})$, where ${t}_{i}$ corresponds to midnight on the $i$th day after January 21, 2020. To account for randomness in case detection and reporting, we replaced $I({t}_{i},{t}_{i+1})$ with ${X}_{i}\sim {\text{NB}}(r,{q}_{i})$, where ${q}_{i}$ is given by Appendix Eq. 43.