A Mathematical Study of COVID-19 Spread by Vaccination Status in Virginia

We introduce a novel n-stage vaccination model and corresponding system of differential equations that stratify a population according to their vaccination status. The model is an extension of the classical SIR-type models commonly used for time-course simulations of infectious disease spread and allows for the mitigation effects of vaccination to be uncoupled from other factors, such as changes in social behavior and the prevalence of virus variants. We fit the model to the Virginia Department of Health data on new COVID-19 cases, hospitalizations, and deaths broken down by vaccination status. The model suggests that, from 23 January through 11 September, fully vaccinated individuals were 89.8% less likely to become infected with COVID-19 and that the B.1.617.2 (Delta) variant is 2.08 times more transmissible than previously circulating strains of COVID-19. We project the model trajectories into the future to predict the impact of booster shots.


Introduction
In December 2019, a novel and highly-contagious coronavirus SARS-CoV-2 was identified in Wuhan, China. The associated disease, COVID-19, causes fever, muscles aches, fatigue, and shortness of breath and has been associated with a significant mortality rate, especially among older and immunocompromised individuals [1]. In the early stages of the pandemic, mitigation measures focused on non-pharmaceutical interventions including travel/border restrictions, social distancing/masking campaigns, work-from-home orders, school closures, and lockdowns. After the Emergency Use Authorization of effective vaccines in the United States at the end of 2020 (Pfizer-BioNTech 11 December 2020; Moderna, 18 December 2020; Johnson & Johnson, 27 February 2021), the response focused on vaccination drives and awareness campaigns. While these efforts have averted the worst-case forecast scenarios, COVID-19 continues to spread widely due in part to the more transmissible B.1.617.2 (Delta) variant, which became the dominant COVID-19 strain in the United States in July 2021 [2].
Mathematical modeling has played a vital role in public policy decision-making since the beginning of the pandemic. SIR-type compartmental models [3] have been used to estimate key epidemiological parameters of infectious disease spread, such as the transmission rate and basic reproduction number [4][5][6]; to evaluate the pathways of spread, such as spread through asymptomatic carriers [7,8]; and to forecast disease dynamics under a variety of non-pharmaceutical intervention strategies, such as lockdowns, contact tracing, and face masking campaigns [9][10][11][12][13][14]. More recent research has focused on aspects of pharmaceutical interventions, such as comparing COVID-19 vaccination strategies (e.g., which segment of society to vaccinate first) and post-vaccination re-opening scenarios [15][16][17][18]. While these and other studies incorporate vaccination into their modeling frameworks, time-course differential equation models have not been widely used to assess the realworld efficacy of vaccination itself. Rather, studies on vaccine efficacy have been limited to controlled clinical trials [19][20][21] and analysis of data aggregates [22][23][24].
In this paper, we introduce a novel n-stage vaccination model that partitions a population according to their vaccination status. This model allows us to concisely assess and forecast the spread of an infectious disease as the vaccination level of a population changes. We fit the model to COVID-19 incidence data from the Virginia Department of Health (VDH) [25] and vaccination uptake data from Our World in Data [26] to estimate (a) the real-world efficacy of vaccination in preventing COVID-19 infections, hospitalizations, and deaths, and (b) the transmissibility of the B.1.617.2 (Delta) variant relative to previously dominant strains. We also conduct a parameter sensitivity analysis and present forecasts of what the spread of COVID-19 could look like through 31 December 2021 given different vaccination coverage levels and with the incorporation of a booster shot.

n-Stage Vaccination Model
We introduce the following n-stage vaccination model, where within each vaccination class i = 0, . . . , n, there are susceptible individuals (S i ), infectious individuals (I i ), hospitalized individuals (H i ), removed individuals (R i ), and deceased individuals (D i ): Note that i = 0 corresponds to unvaccinated individuals and that the removed classes R i may include those who have recovered from illness but also those who have been removed by other means, such as self-isolation or quarantine. Individuals are transferred from the (i − 1) st vaccination class to the ith vaccination class through the ith dose vaccination uptake rate curve V i (t), which is the derivative of the cumulative ith dose uptake curve V i (t). To account for there only being a finite number of vaccination classes, we set V 0 (t) = 0, and V n+1 (t) = 0. That is, we do not include an incoming edge to the 0 th vaccination class or an outgoing edge from the n th vaccination class. We allow for a time-dependent baseline transmission rate for unvaccinated individuals, β(t), so that the transmission rate may vary, for instance, by changes in public health policies or changes in prevalence of COVID-19 variants.
To model COVID-19 spread in the United States, we use the two-stage and three-stage vaccination models: (i) i = 0 for individuals who have not received any vaccine shots (unvaccinated); (ii) i = 1 for individuals who have received one shot of either the Pfizer or Moderna vaccine (partially vaccinated); (iii) i = 2 for individuals who have received two shots of either the Pfizer or Moderna vaccine, or one shot of the Johnson & Johnson vaccine (fully vaccinated); and (iv) i = 3 for individuals who are fully vaccinated and received a booster shot.
The model (1) makes several assumptions: (a) there is no latency period (i.e., susceptible individuals become infectious immediately upon becoming infected); (b) both the vaccine and infection provide long-lasting immunity (i.e., there is no waning immunity or reinfection); (c) only hospitalized individuals may die; and (d) previously infected individuals will not vaccinate. We note that COVID-19 has a mean latency period of 3-7 days [27,28], that studies suggest there is decreased immunity 6 months after becoming fully vaccinated [29][30][31], that an estimated 10-15% of COVID-19 deaths in the United States occur outside of the hospital [32], and that the CDC recommends vaccination for individuals who have previously been infected with COVID-19. We chose not to incorporate these effects for model simplicity; however, an extension of the n-stage vaccination model (1) which incorporates waning immunity and reinfection will be the focus of future work.

Governing System of Differential Equations
Corresponding to (1), we introduce the following governing system of ordinary differential equations for i = 0, . . . , n: To accommodate time-dependent changes in social behavior, public health policies, and the prevelance of SARS-CoV-2 variants, we allow the unvaccinated transmission rate function β(t) to take the following piecewise-constant form: where L is the number of days in the simulation period and m is the number of intervals. That is, we divide the interval of consideration into m parts of equal length and allow each interval to have its own baseline transmission rate for unvaccinated individuals, β j .
To model the vaccination uptake curve V i (t) and corresponding uptake rate curve dV i dt , we use the four-parameter Richards' curve [33]: Explanations and units for parameters for the system of differential Equation (2) with vaccination uptake rate (4) can be found in Table 1. Numerical simulation are carried out with the fourth-order Runge-Kutta method with a fixed step size of h = 0.1 (in days).

Basic and Effective Reproduction Number
The basic reproduction number (R 0 ) of an infectious disease is the expected number of secondary infections produced by a single infectious person entering a fully susceptible population, assuming minimal public health interventions [34]. The effective reproduction number (R t ) additionally takes into account mitigating factors such as public health interventions, natural immunity through infection, and vaccination. These values are considered key epidemiological parameters since R 0 > 1 (R t > 1) suggests that a disease has the capacity to spread in the population while R 0 < 1 (R t < 1) suggests that it will die out. We compute the basic and effective reproduction number of the n-stage vaccination model (2) using the next generation method [35,36]. Table 1. Parameters and rate functions for the n-stage vaccination model (1) and corresponding system of ordinary differential Equation (2) with the ith vaccination dose uptake curve V i (t) (4 none Proportion hospitalized to deceased (ith class) To compute the basic reproduction number, we use the disease free equilibrium . . , n − 1, and I * i = 0 for i = 0, . . . , n, where V i (0) is the number of people who have received the ith shot by the initial time. This gives the following: Note that this number assumes no vaccination beyond the initial time and does not account for natural immunity obtained by previous COVID-19 infections. To account for these mitigating effects, we also compute the effective reproduction number: In general, R t (6) is more challenging to compute since it requires time-course simulations of S i (t), i = 0, . . . , n, whereas R 0 (5) only requires the initial conditions. Note that the outbreak containment condition R t < 1 is explicitly dependent on the effectiveness of the vaccine shots (i.e., α i ), the recovery rates (i.e., γ i ), and the hospitalization rates (i.e., h i ). The condition is also implicitly dependent on the number of individuals who have received each dose of the vaccine (i.e., V i (t)) as these functions affect S i (t) through (2).

Data Fitting Techniques
We fit the n-stage vaccination model (2) with the vaccination uptake curve (4) to data in two steps. We first fit the vaccination uptake curve V i (t) (4) to daily vaccination uptake data by minimizing the sum of squared error. We then incorporate these parameters into the vaccination uptake curve dV i dt in the system of differential Equation (2) and minimize a customized error function. For both steps, we utilize a nonlinear least squares fitter written in Python. We now describe the customized error function for fitting to (2).
We let x = [x i,j,k ] denote the cumulative number of COVID-19 cases in vaccination class i of type j and in week k, and y = [y i,j,k ] denote the weekly increase in that class from the previous week (i.e., y i,j,k = x i,j,k − x i,j,k−1 ). We define θ to be the set of all parameters to be fit andx(t; θ) = [x i,j (t; θ)] to denote the model simulation of (2) for the cumulative cases in vaccination class i of type j and for parameter set θ. We utilize the following normalized error function, which incorporates both the weekly cumulative and rate numbers: where µ j and ν j , j ∈ {I, H, D}, are the normalization factors Note that the normalizations (8) in the error function (7) equalize the weight of the cumulative and weekly rate of increase in infections, hospitalizations, and deaths. For initial conditions, we take We set our initial susceptible classes according to the following: We assume that there are no initial recoveries or deaths of partially or fully vaccinated individuals, so that D 1 (0) = D 2 (0) = R 1 (0) = R 2 (0) = 0. We allow R 0 (0) to be fit into the data to account for the possibility that the best fitting model is one in which there is a different number of overall COVID cases than has been officially reported.

Parameter Sensitivity Analysis
A local sensitivity analysis is a derivative-based approach that is a key tool in exploring parameter identifiability, i.e., whether it is possible to infer a set of parameters from a given data set [37]. Let y(t, θ) be the output variable at time t with vector of parameters θ = (θ 1 , θ 2 , ..., θ k ) from a model with k parameters; then, the normalized sensitivity function of y with respect to parameter θ i is That is, the normalized sensitivity function is the derivative of the output with respect to a specific parameter that has been normalized based on the parameter. The normalized sensitivity functions allow us to determine which parameters when changed by the same percentage result in similar changes in the output variable. Parameters that yield similar changes in the output variable show that those parameters are correlated and will be hard to identify from the data. For the two-stage vaccination model (1), we expect many parameters to be correlated. A sensitivity analysis, however, helps us understand which parameters have the greatest and least influence on an output variable.

Fitting Two-Stage Vaccination Model to Virginia Data
We fit the vaccination uptake curve V i (t) (4) to the cumulative vaccination data from the Commonwealth of Virginia [26] using the methods described in Section 2.4. We used data points starting on 23 January and ending on 2 October and used linear interpolation to fill in missing data. The best fitting parameters and curves V 1 (t) and V 2 (t) are shown in Figure 1. On the left, we indicate the best fit parameters for the generalized logistic curve or Richards' curve (4) to Virginia Department of Health data [25] on vaccination uptake starting with t = 0 corresponding to 23 January 2021. The values M 1 /N and M 2 /N correspond to the saturation proportions for the first and second doses, respectively (i.e., the number of willing/able vaccine recipients). We use the population denominator N = 8,535,519 corresponding to the population of Virginia. On the right, we display the model fit (solid) to the observed data (circles).
We then incorporated these parameters into the vaccination uptake curve (4) and used the methods outlined in Section 2.4 to fit the two-stage vaccination model (2) to the time-course data set VDH-COVID-19-PublicUseDataset-Cases-by-Vaccination-Status (report date 8 October 2021 [25]). This data set is publicly available through the VDH Data Portal and breaks down new weekly COVID-19 infections, hospitalizations, and deaths by vaccination status (unvaccinated, partially vaccinated, and fully vaccinated). We used data from 23 January through 11 September 2021. We used the piecewise-constant transmission rate function β(t) (3) with L = 238 and m = 8 so that L/m = 29.75. This allows the transmission rate to change on a roughly monthly basis. The fit initial conditions and parameters can be found in Table 2. Plots are displayed in Figure 2. Table 2. Best fitting initial conditions (left) and parameters (right) for the two-stage vaccination model (1) and corresponding system of differential Equation (2) with piece-wise constant transmission function β(t) (3), and vaccination uptake curve V i (t) (4) with parameters from Figure 1.  We then incorporated these parameters into the vaccination uptake curve (4) and used the methods outlined in Section 2.4 to fit the two-stage vaccination model (2) to the time-course data set VDH-COVID-19-PublicUseDataset-Cases-by-Vaccination-Status (report date 8 October 2021 [25]). This data set is publicly available through the VDH Data Portal and breaks down new weekly COVID-19 infections, hospitalizations, and deaths by vaccination status (unvaccinated, partially vaccinated, and fully vaccinated). We used data from 23 January through 11 September 2021. We used the piecewise-constant transmission rate function β(t) (3) with L = 238 and m = 8 so that L/m = 29.75. This allows the transmission rate to change on a roughly monthly basis. The fit initial conditions and parameters can be found in Table 2. Plots are displayed in Figure 2. Table 2. Best fitting initial conditions (left) and parameters (right) for the two-stage vaccination model (1) and corresponding system of differential Equation (2) with piece-wise constant transmission function β(t) (3), and vaccination uptake curve V i (t) (4) with parameters from Figure 1.

Efficacy of Vaccination
The fit parameter values α 1 and α 2 in Table 2 suggest that, over the period of study, partially vaccinated individuals are 56.0% less likely to catch COVID-19 while fully vaccinated individuals are 89.8% less likely to catch COVID-19. We can furthermore estimate that, relative to unvaccinated individuals, partially vaccinated individuals are 52.7% less likely to be hospitalized with COVID-19 and 46.0% less like to die from COVID-19, and fully vaccinated people are 88.5% less like to be hospitalized with COVID-19 and 85.7% less likely to die from COVID-19. The reduction in the proportion of people who become hospitalized or deceased relative to the unvaccinated population can be calculated by the following formulas: We note that the proportion of breakthrough cases is expected to rise with a rise in the proportion of the population that is vaccinated; however, the model still supports the notion that, over the period of study, the pandemic was largely a pandemic of the unvaccinated. During the final week of our study from 5 to 11 September, for example, the model predicts that 78.0% of new cases were unvaccinated individuals despite only 33.6% of the population of Virginia (around 2,800,000 people) being unvaccinated.
While this analysis strongly supports the efficacy of vaccination at preventing COVID-19 infections, it does not support the hypothesis that vaccinated individuals benefit from additional protections against severe illness. Rather, the data suggest that the likelihood of an individual who has contracted COVID-19 becoming hospitalized or deceased does not depend significantly on vaccination status. This is contrary to several recent studies [23,24]. We caution that it would be premature to draw conclusions since the two-stage vaccination model and vaccination uptake curves V i (t) do not account for the heterogeneous impacts of COVID-19. We note in particular that a higher proportion of older individuals are vaccinated than younger individuals and that age is correlated with more breakthrough cases and more severe COVID-19 outcomes [1,22].

Transmissibility of Delta Variant
The two-stage vaccination model (1) predicts that the baseline transmission rate of COVID-19 for unvaccinated individuals, β(t), rose significantly from June to August 2021 (see Table 3). This coincides with the emergence in the United States of the Delta variant, which became the dominant strain on the week ending July 3 and made up 98.9% of cases by the week ending 11 September [2]. Although we cannot tell how much of this increase is due solely to the Delta variant, we can compute that the average transmission rate from 23 January through 20 June 2021 was β = 0.214 while the average transmission rate from 21 June through 18 September 2021 was β = 0.445. The time frame in which Delta became the dominant strain of COVID-19 therefore corresponds to a 2.08 times real-world increase in the baseline transmissibility of COVID-19 in Virginia.

Effect of Waning Immunity and Delta Variant on Vaccine Efficacy
Recent studies have suggested that the COVID-19 immunity conferred by vaccination wanes over time [29][30][31] and that the Delta variant may be more likely to produce breakthrough cases [38]. To test for these effects in the VDH data, we introduce the following piecewise-constant vaccination efficacy rates α 1 (t) and α 2 (t) to the two-stage vaccination model (2) with piecewise-constant transmission rate β(t) (3) and vaccination uptake curves V i (t) (4): where L is the number of days the simulation period. That is, we divide the interval of consideration into four parts of equal length, each with a constant rate reduction of partial or fully vaccinated individuals becoming infected with COVID-19. For the simulation period of January 23 through September 18, 2021, this corresponds to L = 238 and L/4 = 59.5 so that the transmission reduction rates is allowed to change roughly every second month. In order to establish population-level waning immunity in the model (2) with timedependent α 1 (t) and α 2 (t) (12), we expect to see 4 . We also expect that the fit would be significantly better than the one obtained with the model (2) and the parameter values in Table 2. In fact, we do not observe any significant trend in the best-fitting α i,j values (α 1,1 = 0.363, α 1,2 = 0.518, α 1,3 = 0.583 α 1,4 = 0.750, α 2,1 = 1.00, α 2,2 = 0.879, α 2,3 = 0.880, α 2,4 = 0.908). We also do not observe a significant reduction in the error function (7). Our analysis consequently does not support including waning immunity or changes in the efficacy of the vaccines against the Delta variant at the population level. We caution, however, against drawing the conclusion that waning immunity is not occurring since the model (2) is only able to estimate vaccine efficacy values at the population rather than individual level. Figure 3 shows the sensitivity results for the cumulative number of cases of COVID-19 in Virginia and the effective reproduction number (6) with respect to the vaccine parameters (M 1 , M 2 , α 1 , and α 2 ) and recovery parameters (γ 0 , γ 1 , and γ 2 ). We interpret the graph for R t in Figure 3 to mean that the perturbation of γ 0 exhibits its greatest influence over R t early in the simulation, with an initial expected decrease of .9% if γ 0 is increased by 1%. We find that cumulative cases and the effective reproduction number are highly sensitive to the efficacy of the second vaccination shot α 2 . Specifically, the model suggests that a 1% increase in α 2 results in a 1% decrease in the cumulative cases (roughly 6500 cases) around the beginning of August. Similarly, we see that a 1% increase in α 2 decreases R t by 1.7% around the beginning of August. The results also suggest that the recovery rate of the unvaccinated individuals, γ 0 , is highly influential. In particular, around the beginning of August, a 1% increase in γ 0 would result in a 2.3% decrease in cumulative cases. The nontrivial impact of vaccination parameters and recovery parameters on R t strongly suggests the impact of changing susceptible populations as vaccination programs are rolled out in the simulation.

Modeling with Different Vaccine Coverage Levels
The sensitivity analysis depicted in Figure 3 suggests that the cumulative cases and effective reproduction number (R t ) are highly dependent on the number of individuals who become vaccinated (M 1 and M 2 ). To further explore this, we increase/decrease the vaccination levels in the model (2) with parameters given in Table 2 (11) on the cumulative cases of COVID-19 in Virginia (left column) and R t (the effective reproduction number, right column) with respect to vaccination parameters (M 1 , M 2 , α 1 , and α 2 ) and recovery rates (γ 0 , γ 1 , and γ 2 ). Note the negative impact of the parameters on the cumulative cases, but a non-trivial impact on the effective reproduction number as population sizes change in each vaccination stage through time. Further notice the small influence that M 1 has on both outputs, until later in the simulation.

Modeling with Booster Shots
The sensitivity analysis depicted in Figure 3 also suggests that the cumulative cases and effective reproduction number (R t ) are highly dependent on the reduction from baseline transmission rate parameters for the two vaccination classes (α 1 and α 2 ). To further investigate the effect that a reduction in transmission due to vaccination may have on the spread of COVID-19, we extend from the two-stage to the three-stage vaccination model (2) where j = 3 corresponds to fully vaccinated individuals who have received the booster shot. We use the parameters from Table 2 for j ∈ {0, 1, 2} and set h 3 = h 2 , γ 3 = γ 2 , and δ 3 = δ 2 . That is, we assume the effects of contracting COVID-19 are similar for somebody with the booster shot as somebody who is fully vaccinated. We estimate a reduction in baseline transmission rate resulting from a booster shot of α 3 ≥ 0.95 [39,40]. We assume that the booster shot uptake curve is similar to the second vaccine dose, delayed 210 days, and that the proportional decrease in people obtaining the booster is similar to that of the first to second shot (i.e., a 3 = a 3 , r 3 = r 2 , ). For the purpose of forecasting COVID-19 spread past September 18, we use the transmission value of the final interval (β(t) = β 8 = 0.453).
In Figure 5, we show forecast simulations for the period 23 January through 31 December 2021 incorporating booster shots, a hypothetical booster uptake schedule, and the utilized parameter values. The booster shot uptake curve V 3 (t) predicts that 40.9% of the population will receive booster shots by 31 December 2021. The model predicts a 19.7-34.3% reduction in new COVID-19 cases from 18 September to 31 December 2021 depending on the efficacy of the booster shot in reducing transmission from baseline (α 3 = 0.95 − 1.00). It is notable that, over the same period, there is a 17.6-31.0% reduction in new cases in unvaccinated individuals and a 14.4-25.2% reduction in new cases for partially vaccinated individuals. This strongly supports the notion that the protection offered by the booster shot, or vaccination in general, is not restricted to the individual receiving the shot; rather, it is distributed across the population and even across individuals of different vaccination statuses.

Discussion
Our analysis suggests that, from 23 January through 11 September 2021, in Virginia, partially vaccinated individuals were 56.0%/52.7%/46.0% less likely to become infected/hospitalized/deceased due to COVID-19 than unvaccinated individuals, while fully vaccinated individuals were 89.8%/88.5%/85.7% less likely to become infected/ hospitalized/deceased due to COVID-19 than unvaccinated individuals (see Table 2 and Figure 2). This strongly supports the efficacy of vaccination in combating the spread and severe outcomes of COVID-19. The model furthermore suggests COVID-19 was 2.08 times more transmissible during the time frame during which the B.1.617.2 (Delta) variant became the dominant strain in Virginia (see Table 3). We did not detect any changes in the vaccine efficacy levels during the time frame during which Delta became dominant in Virginia. This suggests that it is the raw increase in transmissibility of the Delta variant that was the primary driver behind the increase in COVID-19 cases from June to August 2021, rather than a decrease in vaccine efficacy.
We forecast the spread under a variety of different vaccination uptakes levels (see Figure 4) and potential booster shot schedules (see Figure 5). This analysis suggests a five times increase in new cases from 18 September through 31 December 2021, if 10% less of the population had opted to become vaccinated. Our analysis furthermore suggests that a booster shot in the 95-100% efficacy range could reduce the number of new cases between 18 September and 31 December 2021 by 19.7-34.3%, and that this reduction is distributed across vaccination statuses (see Figure 5). These results reinforce the continued push for vaccination as an essential and effective means of containing the spread of COVID-19. It also emphasizes that the benefits of vaccination are distributed across the entire population and not merely restricted to those who are vaccinated.  Figure 5. Simulations of the three-stage vaccination model (2) with j = 3 corresponding to individuals who have received a booster shot. The parameters for the booster shot uptake curve V 3 (t) are shown, and the parameters for the vaccination uptake curves V j (t), j = {0, 1, 2}, are from Table 2. The reductions in the baseline transmission parameter for the booster shot are α 3 = 0.95 (low booster) and α 3 = 1.00 (high booster) [39,40], and the model runs from 23 January through 31 December 2021.

Vaccination Curve Parameters
We note that there are several limitations of our study. Although we did not find evidence that vaccination provides additional protection against severe COVID-19 outcomes such as hospitalization and death, we caution that we did not consider age-stratified differences in vaccination coverage levels, vaccine efficacy, and severe COVID-19 outcomes. In particular, we note that older populations are at higher risk for severe COVID-19 outcomes but also have higher vaccination uptake levels. We also did not incorporate waning vaccine immunity. Although we did not find evidence for decreased vaccine effectiveness over time in our study, we note that our model only considers population-level average immunity and cannot control for when an individual received their first or second vaccination. We also note that, while the vaccination uptake curves V i (t) (4) fit the data well on the chosen time intervals, they do not account for the age-dependent differences in availability of vaccines. The predicted total number of people who will become vaccinated are consequently likely to be below the true value.

Conclusions
We introduced a novel n-stage vaccination model (1) that stratifies a population according to their vaccination status. By controlling for a population's vaccination uptake level, the model is able to effectively estimate the real-world efficacy of vaccination against infections, hospitalizations, and deaths and is able to generate detailed dynamical predic-tions of infectious disease spread under a variety of public health and vaccination uptake scenarios. These features will inform epidemiological research and policy in the years to come as new SARS-CoV-2 variants, such as B.1.617.2 (Delta) and B.1.1.529 (Omicron), emerge and as COVID-19 becomes endemic in the global population. Future work will focus on adapting the n-stage vaccination model into a delay-differential equation model that will incorporate waning vaccine immunity, fitting the n-stage vaccination model to data from a wider variety of jurisdictions and assessing the real-world efficacy of the booster shot as more data become available. Vaccination_Fit.py fits daily vaccination uptake data to a Richard's curve (4) and Model_Fit.py fits the two-stage vaccination model (2) to weekly data on COVID-19 incidence, hospitalizations, and deaths.