Modeling The Impact of Vaccination On COVID-19 And Its Delta Variant

Vaccination is an important means to ﬁght against the spread of the SARS-CoV-2 virus and its variants. In this work, we propose a general susceptible-vaccinated-exposed-infected-hospitalized-removed (SVEIHR) model and derive its basic and effective reproduction numbers. We set Hong Kong as an example to prove the validity of our model. The model shows how the number of conﬁrmed COVID-19 cases in Hong Kong during the second and third waves of the COVID-19 pandemic would have been reduced had vaccination been available then. We then investigate the relationships between various model parameters and the cumulative number of hospitalized COVID-19 cases in Hong Kong for the ancestral and Delta strains of the virus. Next, we compare the evolution of the SVEIHR model to the traditional “herd immunity” threshold where the proportion of vaccinated individuals is static and no further vaccination takes place after model initialization. Numerical results for Hong Kong demonstrate that the static herd immunity threshold corresponds to a cumulative hospitalization ratio of about one percent (assuming the current hospitalization rate of infected individuals is maintained). We also demonstrate that when the vaccination rate is high, the initial proportion of vaccinated individuals can be lowered for while still maintaining the same proportion of cumulative hospitalized individuals.


Introduction
The COVID-19 pandemic has persisted for nearly two years since the end of 2019, resulting in huge human and socioeconomic costs. Since the discovery of the original SARS-CoV-2 virus in late 2019, many measures have been taken in an attempt to regain control over the pandemic, including vaccine development and deployment, changes to healthcare systems operations, and non-pharmaceutical interventions (NPIs) such as mask-wearing, social distancing, and, in extreme situations, lockdowns. However, in the meantime, the original virus has mutated to become even more infectious, especially the Delta variant which has become dominant in many parts of the world. Mathematical modeling 1, 2 is therefore required to capture the dynamics of the COVID-19 pandemic in the presence of vaccinations, NPIs, and the new Delta variant. Important and urgent questions which epidemiological modeling can help address include whether herd immunity can be achieved and the pandemic eventually be eradicated, or whether long-term co-existence with the virus should be pursued instead; whether and when to reopen international or national (e.g., state or provincial) borders; whether additional surge capacity is required at hospitals; and who to prioritize during vaccination (e.g., providing second doses or booster shots versus increasing efforts to give more people their first dose). In the following subsection, we summarize the existing literature on epidemiological modeling of COVID-19 in the presence of vaccination.

Background and related work
A common approach for epidemiological modeling is the use of compartmental models, where each compartment represents a possible state of individuals in the population. The dynamics of the number of individuals in each compartment can then be described as a system of differential equations. A key concept in epidemiological modeling is the reproduction number, defined as the average number of individuals an infected individual will infect. In particular, the basic reproduction number R 0 assumes a fully susceptible population, while the effective reproduction number changes in accordance with the disease dynamics. Generally, a basic reproduction number less than 1 is associated with herd immunity, in which case a disease will eventually die out.
If the basic reproduction number of a disease exceeds 1, vaccination can be used to reduce the effective reproduction number. The vaccination threshold, defined as the minimum proportion of individuals to be vaccinated in order to reach herd immunity, is generally given as where v e is the (mean) vaccine efficacy. However, a problem with this classical formula is that it fails to describe the dynamic interactions between disease transmission and ongoing vaccination during an epidemic, which is important for novel diseases such as COVID-19 where vaccines are initially not available. Matrajt et al. 3 proposed a multi-dose vaccination model for COVID-19 and found that if the single-dose efficacy of a vaccine against COVID-19 is high, then prioritizing the first dose, even at the possible expense of delaying the second dose of two-dose vaccines, may help to contain the pandemic more quickly. Meanwhile, References 4-6 consider optimal prioritization of vaccination groups (based primarily on age). However, 3 and 6 both consider relatively low basic reproduction numbers (R 0 ) for COVID-19, even when all NPIs are relaxed, compared to current estimates of R 0 ≈ 5 for the Delta variant (which is the focus of this paper). Furthermore, the works 4, 5 do not provide estimates of herd immunity thresholds for COVID-19.
Giordano et al. 7 show that NPIs have a higher effect than vaccination alone and advocate for the need to keep NPIs in place during the first phase of the vaccination campaign. However, the authors do not consider the relaxation of NPIs later on during the pandemic when the vaccination ratio is high. We shall consider this latter scenario in this paper. In contrast, Reference 8 considers the optimization problem of minimizing the total impact of social distancing measures over time given certain healthcare capacity constraints. Therefore, the gradual relaxation of NPIs as vaccination ratios increase occurs naturally.
A shared feature of References [4][5][6]8 is that they all employ complex compartmental models, with stratification based on age and possibly other factors such as vaccination status. In this paper, we instead use a much simpler model without the need for many stratified compartments, allowing for easier sensitivity analysis and interpretation of the results. Although many conclusions and recommendations (such as the prioritization of certain groups for vaccination) can only be obtained using stratified models, we demonstrate that some important insights can still be obtained using our simpler model.

Contributions of this paper
The main contributions of this paper are as follows: • We propose a mathematical model to characterize the epidemiological process of COVID-19 with vaccination, and derive formulas for its vaccine reproduction number (similar to R 0 but assuming a vaccinated population) and effective reproduction numbers. We show that, based on the current ratio of CoronaVac (Sinovac) and Comirnaty (Pfizer-BioNTech) vaccinations in Hong Kong (and the computation method can also applicable to other regions), herd immunity cannot be obtained via vaccination alone.
• Using our new model, we observe the impact of vaccination of various parameters (including the initial vaccination coverage at the start of a new outbreak, the rate of new vaccinations, the average vaccine efficacy, and NPI intensity) on the cumulative number of hospitalized cases, for both the original and Delta strains of the SARS-CoV-2 virus.
• We compute, for different vaccination rates and efficacies, the minimum initial vaccination coverage required such that the cumulative number of hospitalized cases is less than a given percentage of the population. Furthermore, we demonstrate how a high rate of vaccinations decreases the required initial vaccination coverage.
• We find that the traditional formula (1) for herd immunity, assuming a static vaccination coverage, corresponds to a cumulative hospitalization ratio of about one percent.
• We show how regions can achieve herd immunity at a lower vaccination coverage ratio than that predicted using (1), but at the cost of additional infections.
For consistency and to gauge the demand for healthcare services, we will maintain the assumption of hospitalization of all confirmed cases even when infeasible due to lack of actual healthcare capacity. Note that although we take the population of Hong Kong as an exapmle in this paper, our model is general and thus applicable to any approximately homogeneous region. Note that a previous version of our proposed model, without vaccination, has been applied to several different regions and found to be quite accurate at modelling the evolution of the COVID-19 pandemic there. In this regard, it is expected that the model used in this paper can be also applied to other cities and countries as well for capturing the dynamics of vaccination. In summary, rather than focusing on the conditions required to achieve herd immunity, we focus instead on the cumulative number of infections throughout the course of an epidemic. One justification for this is that the definition of herd immunity does Transmission rate of exposed individuals α Transmission rate of infected individuals ε V Rate of vaccinated individuals becoming exposed α V Rate of vaccinated individuals becoming infected β Rate of exposed individuals becoming symptomatic γ Hospitalization rate of infected individuals δ Recovery rate of exposed, infected, or hospitalized individuals k Vaccination rate of susceptible individuals v e Vaccine efficacy r V (0) Ratio of vaccinated individuals in the population at time 0, i.e. V (0)/N not discriminate between immunity derived via vaccination or via recovery from infection. Although it is theoretically possible to achieve herd immunity by simply letting the entire population become infected, due to the severity of some COVID-19 cases and the limitations on healthcare capacity, this is not a feasible approach.

Compartmental modeling for COVID-19
In References 9 and 10 an SEIHR model is proposed for the modeling of COVID-19 in Hong Kong and other regions in the absence of vaccinations and found to accurately model the number of confirmed cases in each region. Compared to the classical SEIR model 11 , which contains a susceptible (S), exposed (E), infected (I) and removed (R) compartment, the SEIHR model adds a "hospitalized" (H) compartment for confirmed cases in isolation and makes the E state contagious to model asymptomatic and pre-symptomatic transmissions. By modeling both unconfirmed cases (using the E compartment) and the isolation of confirmed cases (using the H compartment), the SEIHR model achieves better results than the classical SIR and SEIR models when fitting actual COVID-19 infection data. Figure 1. Diagram of the SVEIHR model. The transition rates shown between states are per individual.
In this paper, we modify the previous SEIHR model by adding a vaccinated (V) compartment. We call this new model the SVEIHR model. For simplicity, we assume a single-dose vaccine where vaccination provides immediate protection, with an efficacy of v e . We also assume that the number of daily vaccinations at time t is proportional to S(t), the number of unvaccinated susceptible individuals at time t. The full list of model parameters is given in Table 1 along with some additional notation. The transition rates of individuals between compartments are shown in Figure 1 and the system of differential equations describing 3/14 the SVEIHR model can therefore be written as: where and α, ε, γ, and δ are positive.

Reproduction numbers and herd immunity thresholds
Recall that for the original SEIHR model proposed in Reference 10 , and herd immunity can be reached if 1 − 1/R 0 of the population achieves immunity, either via vaccination or disease recovery. We derive an expression for the reproduction number R V and effective reproduction number R e (t) of our SVEIHR model. There are two compartments, E and I, which are infectious (individuals in H are not considered infectious due to isolation). Let F (x i ), x i ∈ {E, I}, denote the net flow of individuals into state x i due to new transmissions and W (x i ) be the net flow of patients out of state x i due to all other means. For our system (2), we have: We thus obtain the next-generation matrices Let ρ(X) denote the spectral radius of X. Evaluating at the disease-free equilibrium X 0 = (0, N, 0, 0, 0, 0), we obtain the vaccine reproduction number Additionally, the effective reproduction number of system (2) is If NPIs are introduced with a control intensity of c, 0 ≤ c ≤ 1, then α, ε, α V , ε V , R V , and R e (t) are all scaled by a factor of 1 − c and written as:  13,14 for the Delta variant of SARS-CoV-2, we obtain v e 1 = 0.88, and v e 2 = 0.51. From these values, we obtain an average vaccine efficacy of v e = 0.732. (Due to a lack of efficacy data for CoronaVac against the Delta variant, we shall assume that it is unchanged from the ancestral strain of SARS-CoV-2.) Given an estimated R 0 of 5 for the Delta variant of COVID-19 15 , we obtain a vaccination threshold of ρ min = (1−1/R 0 )/v e = 1.092, which implies that herd immunity via vaccination alone is impossible with the current mix of vaccines. In fact, the minimum value of r 1 to achieve herd immunity, assuming a 100 percent vaccination rate, can be obtained by solving ρ min = 1 for r 1 : More generally, if 1 − 1/R 0 ≤ min(v e 1 , v e 2 ), then either vaccine can achieve herd immunity, whereas for 1 − 1/R 0 > max(v e 1 , v e 2 ), then herd immunity cannot be achieved via vaccination alone.
, then herd immunity depends on the proportion of individuals taking each vaccine.

Asymptotic Behavior
Let R(∞) = lim t→∞ R(t) denote the limiting number of removed individuals in system (2) and H c (∞) the cumulative number of hospitalized individuals as t → ∞. We now discuss the relationship between R(∞) and H c (∞). There are three paths from compartment E to compartment R, namely E → R, E → I → R, and E → I → H → R. Let the proportion of exposed individuals that follow each path be p 1 , p 2 , and p 3 , respectively. Since all paths from E lead to R and all paths to R come from E, the cumulative number of exposed individuals, E c (∞), is equal to R(∞). Since compartment H exists on path E → I → H → R only, we obtain H c (∞) = p 3 R(∞). Finally, from Figure 1, we see that

RESULTS
In this section, we present simulation results for our SVEIHR model under various scenarios. For the sake of example, we shall set N = 7394700, the approximate population of Hong Kong as of August 12, 2021.

Revisiting the second and third waves in Hong Kong
We now explore how the second and third waves in Hong Kong would have evolved differently if the Comirnaty and CoronaVac vaccines had been available, with v e 1 = 0.95 and v e 2 = 0.51. Based on Hong Kong vaccination data 12 , we set the ratio of Comirnaty and CoronaVac vaccines administered to r 1 = 0.6 and r 2 = 0.4, respectively. The cumulative number of hospitalized cases over time is shown in Figures 3 and 4 for different initial vaccination ratios r V (0) and vaccination rates k, with initial  Table 2. Table 2. Initial values for the second and third waves of COVID-19 in Hong Kong with vaccination added.

Value
Second wave Third wave The results demonstrate that vaccination has significant impact towards reducing the H c (∞) for both waves, even when the initial vaccination coverage r V (0) is zero (given a high vaccination rate starting at time zero) and despite the relatively lower efficacy of CoronaVac. The reduction in cumulative hospitalizations increases when the initial vaccination coverage or the total vaccination rate is increased.

Sensitivity analysis for the spread of the Delta variant
In this subsection, we use contour plots to investigate the H c (∞) with respect to pairs of parameters in the SVEIHR model related to vaccination and NPIs. The other parameters of the model are set to ε = α = 1.3637, β = 0.2273, γ = 1, and δ = 0.1 for all examples in this subsection, thus obtaining a basic reproduction number of R 0 = 5. This is consistent with the estimate of R 0 = 5.08 obtained in Reference 15 for the Delta variant of SARS-CoV-2. Figure 5 shows H c (∞) with respect to the initial vaccination coverage r V (0) and vaccination rate k, for different control intensities c and a vaccine efficacy of v e = 0.88, as for Comirnaty 13 . Figure 6 shows the same, but with a vaccine efficacy of v e = 0.51, as for CoronaVac 14 . The results show that relatively low numbers of hospitalizations are possible under Comirnatyonly vaccination, even without NPIs for control, whereas a high control intensity is required for CoronaVac-only vaccination regardless of the number of vaccinations. This is consistent with Equation (4) where it is demonstrated that the minimum ratio of Comirnaty vaccinations to achieve herd immunity via vaccination alone is r min 1 = 0.7838, based on the static formula for herd immunity.
Furthermore, it is demonstrated that increasing either c or k decreases the minimum initial vaccination coverage r V (0) required for H c (∞) to remain under a certain limit. In other words, governments do not need to wait for the theoretical vaccination threshold ρ min = (1 − 1/R 0 )/v e to be reached to lift NPIs, as long as the vaccination rate is sufficiently fast.   Table 2. The total vaccination rate is k with a ratio of 0.6 to 0.4.   percent is still required to avoid a large outbreak, as shown by the green transition zone in Figure 5a.  Figure 7 shows the effective reproduction number R e (t) for various r V (0) and v e , given k = 0.002 and no NPIs. It is demonstrated that the lowest vaccine efficacy yields the lowest R e (t) values for sufficiently large t, despite yielding the highest R e (t) at the beginning of the outbreak. This can be explained by noting that a low vaccine efficacy results in a very large number of infections at the beginning of the outbreak, causing greater depletion of the susceptible and vaccinated compartments. In other words, dramatic drops in R e (t), especially early on in an outbreak, may not be a good thing.  Figure 8 shows H c (∞) with respect to r V (0) and v e , for different vaccination rates k and no NPIs. It is shown that a very high vaccine efficacy and initial vaccination coverage is required to control the cumulative number of hospitalizations. Note that for k = 0, the boundaries of the green transition arc, denoting a set of tipping points separating low and high levels of hospitalization, are around r V (0) = 0.8 (for perfect vaccine efficacy) and v e = 0.8 (for perfect initial vaccination coverage), corresponding to the expected value of 1 − 1/R 0 . On the other hand, increasing the vaccination rate k elongates the transition 9/14 arc in the r V (0) axis but not the v e axis. This is because when r V (0) = 1, the value of k has no effect, as all individuals are already vaccinated at the start of the outbreak. Figure 9. Cumulative hospitalizations H c (∞) with respect to the initial vaccination coverage r V (0) and average vaccine efficacy v e , for a simulated Delta outbreak and a vaccination rate of k = 0.005. Figure 9 shows H c (∞) with respect to v e and the control intensity c of NPIs, for different r V (0) and a vaccination rate of k = 0.005. It is shown that for lower vaccine efficacies, increasing c is the only way to control an outbreak of the Delta variant. Additionally, when r V (0) =0.7, there is no value of v e that can control the outbreak without at least some level of NPI control. As a Reference, 17 estimate that mask wearing alone can achieve a control intensity of c = 0.2 (i.e. 80% efficacy) "among compliant subjects".
Further increasing R 0 Figure 10 shows H c (∞) with respect to R 0 and the vaccine efficacy v e , for different r V (0) and a vaccination rate of k = 0.002. We increase the basic reproduction number R 0 by increasing α and ε in fixed proportion. The results demonstrate that the minimum vaccine efficacy required to avoid a large outbreak increases with R 0 (such as the new variant Omicron), until a certain threshold is reached upon which the outbreak cannot be controlled via vaccination alone. This threshold increases with r V (0) , implying that a more transmissible disease requires a larger initial vaccination coverage ratio or stricter NPIs to prevent a large outbreak.

Dynamic versus static vaccination thresholds
We consider for various vaccine efficacies the minimum vaccination threshold required such that the cumulative number of hospitalizations is less than some value h. Note that this is different from the theoretical herd immunity threshold ρ min = (1 − 1/R 0 )/v e . In addition to ρ min , we compute two values: r min V (0) , the minimum initial vaccination coverage for a given vaccination rate k, and r HI , which we define as follows: for a given r min V (0) and k, let t * be the minimum time value such that, if new vaccinations were to cease at time t * , the new number H c (∞,t * ) of cumulative infections remains less than h. Then r HI = V (t * )/N.
The results are shown in Figure 11 for h = 0.01N, no NPIs (c = 0), and various values of the vaccination rate k. The same is shown in Figure 12 for h = 0.1N. As in Figure 10    The results show that when k is high, r min V (0) can be significantly less than r HI and ρ min , meaning that faster vaccination results in a lower initial vaccine coverage ratio required to maintain the size of an outbreak under a given limit. A consequence of this is that governments can consider lifting NPIs even before the theoretical vaccination threshold for herd immunity is reached. Additionally, for h = 0.01N, r HI and ρ min are approximately equal, whereas for h = 0.1N, r HI can be significantly less than ρ min . This is because when h/N is large, a significant proportion of individuals are allowed become infected and gain immunity from the disease in that way, reducing the vaccination requirement to reach herd immunity. However, a large number of hospitalizations may lead to a decline in the quality of medical care due to overcapacity, increasing the case fatality ratio.

DISCUSSION
In this paper, we proposed an SVEIHR model for COVID-19 to capture the dynamic processes of virus transmission and vaccination, and their interactions. Based on the numerical results in this paper, we present the following conclusions: • For the ancestral strain of SARS-CoV-2 (R 0 = 2.26) and a 6:4 ratio of Comirnaty and CoronaVac vaccination, as currently observed in Hong Kong, the theoretical vaccination threshold for herd immunity is ρ min = 0.762, implying that herd immunity is possible via vaccination alone.
• For the Delta variant (R 0 = 5), herd immunity is possible via vaccination alone using Comirnaty, but not using the current ratio of Comirnaty and CoronaVac in Hong Kong, even if the entire population is to be vaccinated -the proportion of Comirnaty vaccinations necessarily has to increase.
• Additional control measures can be used to prevent an outbreak from getting out of control until vaccination has reached a sufficiently high level. Moreover, if the basic reproduction number is high and the vaccine efficacy relatively low, then control measures are required to prevent a mass outbreak even when vaccine coverage is high (in some cases, even if the population is fully vaccinated).
• Faster vaccination results in a lower initial vaccine coverage ratio required to maintain the size of an outbreak under a given limit. Thus governments can consider lifting controls even before the theoretical vaccination threshold for herd immunity is reached.

12/14
• Increasing the number of individuals who are allowed to become infected reduces the requirement on the number of individuals that need to be vaccinated before herd immunity is reached. However, a large number of hospitalizations may lead to a decline in the quality of medical care due to overcapacity, increasing the case fatality ratio.

Limitations of this work
There are still many aspects of vaccination dynamics that need further investigation. For example, our model assumes that immunity via vaccination is immediate and does not consider the reduced vaccine efficacy between the first and second doses of two-dose vaccines or in the period immediately after taking the second dose. We also do not consider the waning of vaccine protection over time or the introduction of booster shots, i.e. additional doses of vaccine to enhance protection, the formulations of which may be also altered to provide extra protection against new variants. Studies 18 have shown that vaccine antibody levels start to wane around 2-3 months after injection of the Comirnaty or Vaxzervria (Oxford-AstraZeneca) vaccines; however, further study is required to show the relationship between vaccine antibody levels and protection against COVID-19 infection. Heterogeneities in the population are also ignored, whereas Reference 19 showed that the herd immunity threshold can be lower in a heterogeneous population than a homogeneous one. We also do not model the differences in efficacy of different vaccines, instead using the average efficacy only, nor do we consider scenarios with multiple disease variants in co-existence, mainly due to the current dominance of the Delta variant in most of the world. Migration between heterogeneous geographical regions, which may have different vaccination rates, vaccine efficacies, and control measures, is also not considered here. This will become more important in the near future as international travel resumes.
Finally, it may be useful to model not only total "hospitalizations", but also the number of patients requiring isolation and observation only, basic treatment, mechanical ventilation, and intensive care, due to the different healthcare requirements for each case.
In future work, we will extend our model to include the effects of heterogeneity with respect to age, social activity, and administered vaccines. Other effects for possible consideration include waning vaccine efficacy over time, severe or critical cases as a subset of hospitalized cases, and the co-evolution of multiple disease strains.