Abstract
The introduction of COVID-19 vaccination passes (VPs) by many countries coincides with the Delta variant fast becoming dominant across Europe. A thorough assessment of their impact on epidemic dynamics is still lacking. Here, we propose the VAP-SIRS model that considers possibly lower restrictions for the VP holders than for the rest of the population, imperfect vaccination effectiveness against infection, rates of (re-)vaccination and waning immunity, fraction of never-vaccinated, and the increased transmissibility of the Delta variant. Some predicted epidemic scenarios for realistic parameter values yield new COVID-19 infection waves within two years, and high daily case numbers in the endemic state, even without introducing VPs and granting more freedom to their holders. Still, suitable adaptive policies can avoid unfavorable outcomes. While VP holders could initially be allowed more freedom, the lack of full vaccine effectiveness and increased transmissibility will require accelerated (re-)vaccination, wide-spread immunity surveillance, and/or minimal long-term common restrictions.
Introduction
In the past, governments have required proof of vaccination for travel, with yellow fever being the best-known example, and the only disease for which a certificate is needed as a precondition of entry to a country in compliance to the International Health Regulations [1]. However, the idea that proof of vaccination will become a prerequisite for crossing borders or to enter facilities, visit businesses premises, participate in events, and generally enjoy more freedom, has only arisen in the context of combatting the COVID-19 epidemic. Despite technical challenges, scientific uncertainties, and ethical and legal dilemmas, the idea of VPs, i.e., documents issued on the basis of vaccination status, is now receiving unprecedented attention [2, 3, 4]. The Commission of the European Union (EU), in an effort to ensure a uniform pan-European approach, as similar initiatives for VPs were emerging at national level, put forth a proposal for a framework of issuing, verifying and accepting interoperable vaccination certificates to be implemented across the EU [4], along with a corresponding proposal for third-country nationals residing in the EU [5]. The proposal, in its amended form, for the ‘Digital COVID Certificates’ (DCCs), took effect on July 1, 2021. Many consider the EU DCCs, and other forms of VPs in general, as tools to restore people’s freedoms and increase well-being, whilst allowing economies to reopen. Finally, even without VPs, vaccinations alone may result in less stringent behavior. Those vaccinated may feel more secure and restrict themselves less from contacts they would refrain from when not being vaccinated.
The introduction of VPs and consequent changes in behavior coincided with the emergence of new variants of concern of the virus [6]. Notably, the Delta variant (B.1.617.2) was detected in many countries across Europe, causing a resurgence of COVID-19 in the United Kingdom at a startling pace [7, 8]. Delta has been estimated to be 50% more transmissible than the Alpha variant (B.1.1.7), already estimated to be 50% more transmissible than the parental strain [9, 10, 11].
Evidence indicates vaccine effectiveness can greatly vary [12, 13] and it may be compromised due to escape variants [14] and waning immunity [15, 16, 17, 18]. Preliminary data from several countries indicate reduced vaccine effectiveness against the infection with the Delta variant compared to the Alpha variant [19, 20, 21], even as low as 64% for the Comirnaty (Pfizer-BioNTech) vaccine according to unpublished data from Israel [22]. Emerging evidence suggests that the vaccines are still highly effect at preventing serious illness and hospitalization [11, 20, 21].
Still, avoiding another COVID-19 infection resurgence remains a valid and potentially attainable goal [23]. An estimated 10% of COVID-19 infections will have long-term sequelae (long COVID), posing an increasing threat to national health systems [24, 25]. Finally, large numbers of infected create a large pool of virus hosts, resulting in more replications of the virus and higher chances of emergence of mutations conferring evolutionary advantage, including increased transmissibility and antigenicity. To detect the emerging variants, wide-spread surveillance of genetic and antigenic changes in the virus population has to be conducted, together with experiments elucidating their phenotypic implications [26]. Such needed comprehensive surveillance and experiments may become stalled for a large population of infected. Given these circumstances, it is critically important to understand the impact of key risk factors such as: vaccine ineffectiveness, slow vaccination rate, waning immunity, fraction of individuals in the population who will never become vaccinated, and finally the levels of restrictions, on infection dynamics. Not being aware of the risks and their consequences, and a false sense of security, including when approaching higher vaccination coverage, may result in policymakers opting to select suboptimal levels of restrictions.
Various models have been developed to inform vaccination strategies [27, 28, 29, 30, 31, 32, 33, 34, 35]. One such effort indicates lower vaccine effectiveness coupled with an increase in social contact among those vaccinated (behavioral compensation) may undermine vaccination effects, even without considering immunity waning [36, 37]. Scenarios for the post-vaccination era were also considered by Sandmann and colleagues (2021), finding that under realistic scenarios periodic epidemics are likely [38]. So far, there has been no model to focus on the medium- and long-term impact of relaxing restrictions for VP holders, with due consideration to vaccine effectiveness, durability of response, and vaccine hesitancy, especially in the context of the increased transmissibility of the Delta variant. Given the implementation of the EU DCC, and emerging heterogeneous measures on utilising the VPs for different purposes at national level by establishing different levels of freedom for VP holders in terms of accessing premises, facilities, travelling within a country, etc., it is important to examine the broad parameters determining how to optimise the implementation of measures such as the EU DCC and other VPs.
To address these needs, we propose a mathematical VAP-SIRS model, which accounts for key parameters affecting current infection dynamics, such as for lower restrictions for VPs holders than for the rest of the population, imperfect vaccination effectiveness against infection, rates of (re-)vaccination and waning immunity, and the increased transmissibility of the Delta variant, which all impact the effective reproduction number of the virus. The model predicts the impact of restrictions for VP holders and the rest of the population on epidemic thresholds for various parameter settings, and delivers a systematic framework to assess key considerations for policymaking.
Results
The VAP-SIRS model of the impact of COVID-19 VPs
The VAP-SIRS model extends the classical SIRS model [39] (red arrows in Figure 1a) with additional states and parameters that describe the dynamics of vaccination rollout in a population (green arrows in Figure 1a). To this end, we consider the following subpopulations: (i) initially susceptible SN, who, if successfully vaccinated, populate the immune group V , with rate av, where v is the vaccination rate and a is the vaccination effectiveness, (ii) susceptible who were vaccinated but did not gain immunity (S1), (iii) vaccinated, whose immunity waned with rate ω and who became susceptible again (S2), (iv) susceptible, who are not and will never get vaccinated (SD). The SD compartment contains people who for health reasons cannot receive current types of vaccines, as well as individuals who do not get vaccinated because of hesitancy, beliefs or other individual reasons. The fraction of the population that will never be vaccinated is denoted by d. Additionally, revaccination of S2 populates V with rate avr. All recovered, unless in the recovered compartment RD, are also subject to vaccination. Before the recovered in the RV lose immunity, they might be revaccinated, and, thus, populate the V group with rate vr (similarly, RN are vaccinated with rate v). In this case, vaccination effectiveness is fixed to 1, which is substantiated on the basis of the fact that vaccination combined with a previous infection should confer a much stronger protection than only vaccination of a susceptible individual. Across the manuscript, we assume the revaccination and vaccination rates are equal, vr = v.
The presented model analysis is performed for carefully selected parameter setups. We consider two different vaccination rates a, 0.004 and 0.008 doses per person daily, chosen on the basis of the current rates observed in Europe [40, 41]. As vaccine effectivenesses for the Delta variant, we consider 0.6 and 0.79, which were reported as the effectivenesses of the most widely used vaccines: Vaxzevria (AstraZeneca) and Comirnaty (BioNTech/Pfizer) respectively for this variant [20]. For the Alpha variant, the effectivenesses of the same vaccines for that variant are considered instead, namely 0.79 (Vaxzevria) and 0.92 (Cominarty) [20]. We consider realistic fractions d of never-vaccinated equal to 0.12 (optimistic), and 0.3 (pessimistic), reported for the United Kingdom and France, respectively [https://ourworldindata.org/, as of June 15th, 2021]. Furthermore, two post-vaccination immunity waning rates ω are considered corresponding to optimistic (500 days; ω = 1/500) and pessimistic (200 days; ω = 1/200) average immunity duration periods, reflecting emerging data on large individual variation of immunity waning and other key factors influencing this process [15, 18, 42, 43, 44]. There remains uncertainty regarding the waning time for natural immunity, and whether it varies between the different SARS-CoV-2 variants, but early evidence indicates it lasts at least 180 days [45, 46, 47]. Hence, we consider an optimistic scenario of natural immunity lasting on average similarly long as the optimistic immunity gained via vaccination: 500 days (corresponding to natural immunity waning rate κ = 0.002). Based on the current studies, we fix the generation time to 6 days (γ= 1/6) [48, 49].
We assume that VP holders are all those who completed at least one complete vaccination cycle, i.e., one dose or two doses depending on the vaccine used (Fig. 1), which is also the basis on which the EU DCC is issued. The restriction level (ranging from 0 to 1) is introduced as a modulator of the SARS-CoV-2 reproduction number. Here, we consider that without any restrictions, the basic reproduction number for the B.1.617.2 variant (Delta) is equal to 6 (an optimistic estimate based on [50, 51]), while for the B.1.1.7 variant (Alpha) an optimistic estimate is equal to 4 [48, 49]. Two levels of restrictions are considered: restrictions fv for contacts among VP holders, as well as restrictions f for contacts of the VP holders with the rest of the population and for contacts within the rest of the population. The impact of VPs is studied assuming that fv < f: a VP holder has more freedom of contact with other VP holders, or is generally subject to fewer restrictions on the VP holders than the rest of the population. Importantly, in general f and fv should be interpreted as the net effect of all combined factors that reduce the reproduction number of the virus within the respective groups: all applied non-pharmaceutical interventions, including testing and isolation, together with the resulting changes in behavior. The situation where no VPs are implemented, hence the vaccinated have the same restrictions as the rest of the population, and there are no changes of behavior due to vaccination, is modeled by fixing fv = f. Finally, to analyze the impact of social behavior, we consider two types of mixing between subpopulations: proportional (typical for SIR models) and preferential, where the VP holders prefer contacts with other VP holders. See Methods for a detailed model description.
VAP-SIRS predicts a possible infection resurgence despite vaccinations
VAP-SIRS predicts unfavourable epidemic dynamics for a wide range of parameters, both for the Delta and the Alpha variants. As an example consider the Delta variant, and vaccine effectiveness a = 0.79 (the effectiveness of the Comirnaty vaccine against the Delta variant), (re-)vaccination rates v = vr = 0.004, low never-vaccinated fraction d = 0.12 (reported for the United Kingdom), low immunity waning rate ω = 0.002, low natural immunity rate loss κ = 0.002, and proportional mixing. This set of parameters corresponds to a seemingly safe setup, which we will call the reference setup. The impact of various parameter changes with respect to this reference will be considered below. For such a setup, consider medium-high restrictions level f = 0.77 for contacts of the VP holders with the rest of the population and within the rest of the population, along with a restrictions reduction for the VP holders compared to the rest of the population by around 30%, resulting in medium restrictions fv = 0.55 for VP holders. For these parameters, the model predicts a small wave of infections shortly after the vaccination program starts, followed by a large wave later (color curves Fig. 1b). This behavior is explained by the population structure (Muller plot, Fig. 1b) and can only happen due to the different levels of restrictions for the VP holders and the rest of the population. In this scenario, the first wave is driven by the unvaccinated susceptibles (SN) and suppressed by ongoing vaccination, as expected. Interestingly, the second, larger wave is driven by the SV group. The SV group is composed of the number of individuals for whom the vaccine was ineffective (S1) and those vaccinated who lose their immunity and are not yet revaccinated (S2).
Stability analysis identifies potential scenarios for the COVID-19 epidemic depending on the restrictions imposed on VP holders and the rest of the population
To assess the epidemic evolution in different scenarios, we analyse stability by linearising the model equations with I = R = 0 and introduce the instantaneous reproduction number R* (see Methods). R*(t) is the reproduction number that would be observed at time t, given the restrictions f = (f, fv) and the composition of the population, where the number of infected is very small. For R*(t) > 1, switching to f at time t results in an overcritical epidemic evolution, with an initially exponential growth of infections; for R*(t) < 1, switching to f at time t results in a subcritical epidemic evolution, where the number of active cases decreases to zero. The R* is more informative of epidemic thresholds than the standard effective reproduction number, as it does not depend on the actual number of infected and recovered.
Assuming the reference setup for the Delta variant, we consider five choices of restriction combinations (prototypical for five regions of the parameter space, see Figure 2), leading to different time profiles of R* (Fig. 1c). As control setups, we introduce two settings that represent policies in a given population without the implementation of the VPs, one with a common high restriction level f = fv = 0.92, which keeps the epidemic subcritical (scenario denoted -, green dot-dashed line in Fig. 1c), and one with a common medium-high restriction level f = fv = 0.77, which results in a time evolution fo R* from overcritical to subcritical (denoted +-, yellow dot-dashed line in Fig. 1c). In such settings, in case VPs are introduced, VP holders can gain low (<20%), medium (20-50%) or high (> 50%) restriction reduction with respect to the restrictions for the rest of the population. Granting high restriction reductions to VP holders, both with mid-high and with high restrictions enforced for the rest of the population, eventually leads to an overcritical epidemic (red and pink curve in Fig. 1c: the red curve shows a persistent overcritical epidemic, a scenario denoted -, while the pink curve shows an epidemic that is initially subcritical and then becomes overcritical, a scenario denoted -+). Medium restriction reductions for VP holders, along with high restrictions for the rest of the population, yield a subcritical epidemic evolution (another example of scenario -, blue curve in Fig. 1c). When mid-high restrictions are enforced for the rest of the population, a medium restriction reduction for VP holders leads to an epidemic that is initially overcritical, then becomes subcritical and after a few months switches to overcritical again, starting a new wave of infections (orange curve in Fig. 1c; denoted +-+, this is also the scenario shown in the simulation in Fig 1b). Finally, always with mid-high restrictions enforced for the rest of the population, a low restriction reduction for VP holders leads to an epidemic that is initially overcritical and then switches to subcritical (another example of scenario +-, cyan curve in Fig 1c).
In each scenario we computed the time evolution of the instantaneous doubling time D, capturing how fast the infections grow. For a given f, D(t) is the doubling time that would be observed for the growth of a small initial number of infections at time t, with enforced restrictions f. Very short doubling times, below 30 days, can be observed in three scenarios that are (eventually) overcritical: see the red, orange and pink curves in the Supplementary Fig. S1.
Flexible measures are required to avoid epidemic resurgence depending on parameter setups
The relevant f − fv parameter space, where fv ≤ f, can be divided into five regions, where the epidemic dynamics follows the distinct patterns exemplified in Fig. 1c. Fig. 2 shows the impact of changing specific single parameter values on the expected scenarios and on times to critical events, tracking time up to two years. The area occupied by each region changes depending on the parameter setups. For example, in the reference setup for the Delta variant (a = 0.79 - the Comirnaty effectiveness on the Delta variant, v = vr = 0.004, d = 0.12 - the fraction of never vaccinated in the United Kingdom, ω = 1/500, κ = 1/500, and proportional mixing) in Fig. 2a, the overcritical region (denoted +, with R* always above 1) occupies the lower left corner. This region is enlarged in the case of a lower vaccine effectiveness (a = 0.6 - the effectiveness of Vaxzevira on the Delta variant, Fig. 2b), and higher waning rate (Fig. 2f). In contrast, it shrinks with a higher vaccination rate (Fig. 2c), indicating that there is a concrete benefit from deploying efficient vaccination programs. The subcritical region (−, with R* always smaller than 1) lies in the opposite corner of the f − fv space, for larger restriction values, and, for a fixed fraction of never-vaccinated d, tends to decrease for setups where the overcritical region increases. As expected, the switch to a larger fraction of never-vaccinated (to d = 30%, corresponding to the reported fraction in France), increases the overcritical (+) region (Fig. 2e). But, at the same time, the larger fraction of never-vaccinated increases also the subcritical (−) region. This is due to the fact that the never-vaccinated are assumed to follow stricter restrictions, compared to VP holders, and therefore their larger fraction can constrain the emergence of the later waves, characteristic of the regions +−+ and −+. Still, a strategy relying on this effect might be difficult to implement due to the large + region and can lead to undesirable outcomes in practice.
Inside each of the three regions associated with the +-+, -+, +-scenarios in Fig. 1c, the specific parameter settings differ by the time to the critical threshold of interest for that region (the last observed switch between subcritical and overcritical epidemic, which for the +-+ region, for example, is the second critical threshold; see Methods for the computation of the times to critical thresholds). For the reference setup (Fig. 2a) and the +-+ region, the critical threshold is reached after a minimum ∼4 months. Decreasing the vaccine effectiveness from Comirnaty’s to Vaxzevira’s (Fig. 2b), as well as increasing the waning rate (Fig. 2f), leads to overcriticality sooner, after a minimum of ∼2 and ∼3 months respectively, for low fv values. Increasing vaccination rate (Fig. 2c) shrinks the +-+ region. The comparison between proportional and preferential mixing shows the impact of more intense interactions of the VP holders inside of their own group, and less intense contacts of the VP holders with the rest of the population. With preferential mixing (Fig. 2d), the +-+ region becomes larger and overcriticality is reached even sooner. This is due to the fact that preferential contacts among VP holders accelerate the emergence of the wave caused by infections of the VP holders. Seemingly counter-intuitively, increasing the number of never-vaccinated people (Fig. 2e) shrinks the +-+ region and delays the onset of overcriticality. This is due to the fact that the onset of overcriticality in the +-+ region depends not only on the intensity of contacts of the VP holders, but also on their fraction in the population; with a larger fraction of never-vaccinated, the fraction of VP holders in the population decreases.
The above analysis of the different regions predicts a possible switch to overcritical epidemic growth for a given parameter setup and, if there is a switch, it provides the time it happens, counting from the onset of the vaccination program. It does not, however, indicate how fast the overcritical growth will be. To inform about what growth rates can be eventually expected in the overcritical regime, we compute the asymptotic R* (the R* (t) for t → ∞, see Methods) for all parameter setups and all combinations of restrictions in the relevant f − fv space. For a given restriction combination f, the asymptotic R* indicates how quickly the infections grow shortly after the restrictions are set to f in the asymptotic state. For all considered parameter setups, except for the one with high (re-)vaccination rate, and for all except the +- and the - regions, large asymptotic R* can be expected, which corresponds to short doubling times (Fig. 2). This analysis highlights the importance of avoiding the overcritical (+) region, as there the asymptotic R* values can even exceed 2 when the restrictions are low.
Comparing Figure 2 to Supplementary Figure S2 shows how the Delta variant worsens all scenarios with respect to the Alpha variant: in all panels of Figure 2, the Delta variant leads to a considerable expansion of the overcritical region, shrinking of the safe subcritical region, and to consistently larger values of asymptotic R*. This is due not only to a higher transmissibility of the Delta variant, but also due to the fact that the considered vaccines have lower effectiveness for this variant, as compared to the Alpha variant.
We further investigate how the expected scenarios, times to critical events (tracking time up to two years), and asymptotic R* values are affected by changes of two parameters at once, compared to the reference setup, for the Delta (Fig. 3) and the Alpha variant (Supplementary Fig. S3). The double parameter changes give insight into the possible synergistic and compensatory effects between individual parameter changes. Compared to the effect of only decreasing the vaccine effectiveness from Comirnaty’s to Vaxzevira’s (Fig. 2b), the effect of jointly decreasing the vaccine effectiveness and increasing the vaccination rate (Fig. 3a) indicates that a higher vaccination rate can compensate to some extent for the loss of effectiveness. Similarly, an increased vaccination rate can counteract increased immunity waning rate (Fig. 3e). The combination of decreased effectiveness and increased immunity waning rate has the worst effect, as it largely increases the overcritical region (+), decreases the subcritical region (-) and accelerates the times to the overcriticality in all other regions (Fig. 3c). Finally, combinations of an increased never-vaccinated fraction with other parameter changes show an interesting mix of effects. When both the never-vaccinated fraction and the vaccination rate increase, the overcritical (+) region decreases and the subcritical region increases, while the times to overcriticality in the +-+ and the − + regions increase (Fig. 3d). Similarly, there is a synergistic effect of the combination of the increased never-vaccinated fraction and the increased immunity waning rate (Fig. 3f). For the Alpha variant, the effects of coupled parameter changes combine the same way as for the Delta variant, but once again it is apparent that, for all the parameter setups we considered, with the Alpha variant much less restrictions are required to avoid epidemic resurgence than with the Delta variant (Supplementary Fig. S3).
Taken together, these results indicate that, unless novel vaccines with higher effectiveness are invented and distributed, and unless much faster and wider vaccination programs are implemented, resulting in much more favorable parameter settings than the realistic ones analysed here (including those considered optimistic), highly unfavorable infection dynamics are likely to emerge for the Delta variant, and less, but still, for the Alpha variant. The -+ and +-+ regions in Figures 2 and 3 can seem attractive as restriction policies, because they entail larger freedom for the VP holders; both these regions, however, eventually result in epidemic resurgence and either should be avoided or the time spent in these regions should be very carefully regulated. For example, if sufficient restrictions are enforced for the rest of the population, the VP holders may initially be granted additional freedoms (larger if the Alpha variant is dominant in the population, and much lower if the Delta variant is dominant), which corresponds to the -+ region. In this way, an overcritical situation (region +) will be avoided. However, to prevent the epidemic from becoming overcritical after an initial decline in case numbers, restrictions on VP holders need to be timely increased and adapted, to avoid spending longer time in the -+ region than the time to overcriticality. Thus, moving out of the -+ region to the +-region with the right timing could be one of possible strategies. It may, however, be more practical to circumvent many changes of restriction policies over time and it may be fair for everyone to face the same restrictions. Safe common restrictions, however, corresponding to the parameters on the diagonal in the subcritical (-) region or Figures 2 and Fig. 3 and Supplementary Figures S2 and Fig. S3, are relatively high, especially those required by the Delta variant, and may therefore cause unrest in the population.
A minimum common restriction level can keep the epidemic subcritical in the long-term
We compute the minimum common restriction level fmin for the whole population that would guarantee to avoid an overcritical epidemic in the long-term (for time approaching infinity, Methods): where Vas as is the asymptotic fraction of the immunized in the population The resulting values differ depending on the setups of vaccine effectiveness a, revaccination rate vr, the fraction of never-vaccinated population d and immunity waning rate ω (Tab. 1). The minimum common restrictions for the reference setup are equal to fmin = 0.69. Out of parameter setups with single change compared to the reference, doubled (re-)vaccination speed leads to the lowest possible common restriction level. Even for this most optimistic setup (high a = 0.79, high vr = 0.008, low d = 0.12, low ω = 0.002; Tab. 1 third row) we obtain Vas = 0.6, and fmin = 0.62. The level of 0.62 restrictions is around twice as high as the level 0.29 that would be required for the Alpha variant (Supplementary Tab. 1), and is a considerable reduction of freedom compared to before the pandemic. It is noteworthy that in the long term, to avoid infections rising, minimum common restrictions have to be increased to 0.74 with the larger fraction of never vaccinated d. Thus, a scenario with a large fraction of the population without immunity gained via vaccination requires long-term high restriction levels, and as such seems politically unfeasible.
When changing two parameters simultaneously in order to assess synergies, we find that a decreased vaccine effectiveness or an increased share of never vaccinated or an increased waning rate can barely be offset by an increase in vaccination speed. Both a decreased vaccine effectiveness and an increase in the share of never vaccinated in combination with an increased waning rate considerably increase the minimum restriction level that is adequate to ensure resurgence can be avoided. The latter (increased d, increased ω as compared to the reference) is the most pessimistic of the considered scenarios, with fmin = 0.8.
This analysis highlights the importance of vaccine effectiveness, vaccination speed, but also of the fraction of the never-vaccinated. Such demanding requirements for stringent minimum common restrictions could be reduced if novel vaccines with higher effectiveness become available, if faster and wider vaccination programs are implemented, and finally, if the never-vaccinated fraction shrinks.
Endemic state analysis reveals the possibility of large daily infection numbers
For a given restriction combination f, the above analyzed asymptotic instantaneous reproduction number R* (Figures 2 and 3, Supplementary Figures S2 and S3) indicates how quickly the infections grow shortly after the restrictions are set to f in the asymptotic state; however, it does not provide insight into the daily infection numbers the system converges to. To this end, we compute the daily infection numbers both in the vaccinated and the unvaccinated subpopulations in the endemic state, as functions of the restrictions f for the Delta variant (Fig. 4) and compare it to the scenarios achieved with the Alpha variant (Supplementary Fig. S4). In contrast to the computation of the instantaneous reproduction number R* and its asymptotic values, which is based on the analysis of the linearized system of the ordinary equations in the VAP-SIRS model, the endemic state is based on the computation of the stationary point of the full set of the equations (Methods).
For all parameter setups, in all regions apart from the subcritical (-) region, the daily infections in the endemic state will exceed 10 per million, which is the tolerance threshold for efficient test, trace and isolation policy [52]. For the setups that correspond to low vaccination effectiveness or short waning time, the endemic state is most unfavorable, as the daily infections can exceed 1000 daily cases per million. A high (re-)vaccination rate is crucial to expand the safe region (Fig. 4c). A sharp transition can be seen between favorable and unfavorable parameter setups. In the endemic state, the daily infection numbers in the vaccinated subpopulation can exceed that of the unvaccinated subpopulation, which underlines the risks of waning immunity.
Considering the parameter setups that arise by changing two parameters at once with respect to the reference setup gives insights about their joint effects, shown in Fig. 5; the effect of simultaneous parameter variations is akin to that described earlier for the values of asymptotic R* and time to critical thresholds in Figure 3.
Again, comparison with the endemic infection numbers predicted for the Alpha variant (Supplementary Fig. S4 and Supplementary Fig. S5) shows that Delta has considerably narrowed opportunities to reduce restrictions for the VP holders, underlining the negative impact of the higher transmissibility of the Delta variant and lower effectiveness of the vaccines on this variant.
Discussion
Introducing VPs is widely seen as a means to opening up economies and societies, despite the ongoing epidemic. A recent complication in this respect is the rise of the Delta variant with its higher transmissibility and decreased vaccine efficacy. To inform this discussion, we extend a SIR model to reflect vaccination dynamics and possibly different restrictions for VP holders, with empirical parameters for both the Alpha and Delta variant.
VAP-SIRS deliberately keeps several aspects simple. The model is not compartmentalized for age groups and does not explore the impact on hospitalisation or intensive care unit utilization like some other models, albeit in the context of exploring different parameters than larger freedom for VP holders [27, 29, 31, 32, 33, 34]. Deaths could be taken into account through a straightforward modification of the model, which would however lead to more parameters. In this case, the fact that vaccination reduces the risk of death would have to be accounted for. The impact of deaths, need for intensive care, or long-COVID cases on public health and society, which is very important especially when the epidemic becomes overcritical, can be indirectly assessed on the basis of the number of cases. In general, features such as age groups do not add further insights into the questions which are the focus of our study, namely, the impact of VPs on infection dynamics and the parameters that increase the risks of overcritical dynamics, given the spread of the Delta variant. In this context, the advantage of our model is that it is enriched in features such as revaccinations and waning immunity, which are particularly relevant in the long term. Avoiding another wave is a prudent goal due to the threats it poses, in the form of long-term health effects, the deleterious impact on societies and the emergence of new variants.
Nevertheless, possible extensions to our model could include interindividual or age-dependent variations in immunity, which would render it relevant for people with severe chronic conditions or immunodeficiencies. The presented analysis has been performed assuming that without restrictions, the maximum reproduction number of the virus is Rmax = 6 or Rmax = 4 for Delta and Alpha variants, respectively. More transmissible variants could easily be modeled by fixing higher values of Rmax. Possible future variants, for which existing vaccines may potentially be less effective could be considered using our model by fixing smaller vaccine effectiveness parameter a than the values we considered. We also do not account for seasonality, which seems to have a dampening effect on epidemic dynamics during the summer months, when it is possible to temporarily reduce restrictions. Not all analyzed parameter values are exactly known, such as the post-vaccination or natural immunity waning time. We, however, fix optimistic values for such parameters, and show that unfavorable infection dynamics can still be obtained even under optimistic assumptions.
Despite limitations, our model accounts for key parameters influencing infection dynamics and gives valuable insights into policies pertaining to the introduction of VPs, contributing to render the valid goal of avoiding resurgence attainable. We find that a wide range of the VAP-SIRS model parameter choices, even optimistic ones, show the possibility of an epidemic resurgence for both variants. The risk of resurgence is higher in the case of implemented VP, i.e., with VP holders enjoying reduced restrictions, such as being exempt from wearing masks and testing before entering business, gastronomical, or tourist premises. The resurgence can be avoided in the short and in the long run only when the restrictions are kept high for the rest of the population, and the reduction for the VP holders is moderate or small, especially for the Delta variant. The main driver of this phenomenon is the potential lack of immunity of VP holders. With a VP, people enjoy lower restrictions, while some actually remain both susceptible and potentially contagious because the vaccine was ineffective or the immunity has waned.
For all analyses, a comparison between values for the Alpha and Delta variants shows that Delta has drastically worsened all scenarios. One illustrative finding is that the minimum level of common restrictions to avoid resurgence in the reference setup has doubled from a low 0.29 (Alpha) to medium 0.62 (Delta).
Changing key parameters such as vaccine effectiveness, (re-)vaccination rate, or waning immunity rate to realistic levels found in studies or certain countries shows the expected effect these changes would have on infection dynamics. We quantified these effects by evaluating the times to overcriticality, asymptotic instantaneous reproduction number R*, minimum necessary common restriction level that avoids resurgence in the long term, and numbers of cases per million in the endemic state, for the relevant range of possible restrictions for the VP holders and the rest of the population. As expected, the model shows that there is a larger selection of admissible restrictions’ combinations under high vaccine effectiveness, low share of never vaccinated, a higher (re-)vaccination rate, slowly waning immunity, and proportional social mixing. For the Delta variant, however, and even for optimistic of these parameter setups, the room for manoeuvre in terms of lowering the restrictions is very small. Moreover, not all of these parameters are amenable to policy action. In a nutshell, our results consistently suggest that with the Delta variant and with the way the vaccination program and introduction of VPs is currently implemented, unfavourable developments of the epidemic are likely, and to counteract these developments and to maximize possible freedoms for their citizens, decision makers should exploit all possibilities to enhance the development of effective vaccines, increase vaccination speed and the number of vaccinated.
It is noteworthy that VP holders are less likely to be tested, as they are assumed to be protected and they may exhibit milder symptoms. Therefore, their potential infection is more likely to remain undetected, resulting in an effect similar to that of lowering restrictions. To prevent undesirable outcomes, the testing and quarantine criteria should be applicable also to the VP holders. Testing should aim at detection of vaccinated people that have lost, or have never gained, immunity. Finally, temporary VPs could be considered, with their prolongation conditioned on high antibody level or recent (re-)vaccination.
The utilisation of tools such as the VAP-SIRS model, along with different tools available to policy-makers should be explored in the context of monitoring the implementation of VPs, including the EU DCC measures, to ensure optimisation of key parameters. In this manner, evidence-informed policy-making would be safeguarded as would the best possible outcomes in terms of effectively combatting the current pandemic.
Methods
Mathematical model
We introduce a modified susceptible-infectious-recovered-susceptible (SIRS) model [39] (Fig. 1a). The population is divided into two subpopulations: those who are not vaccinated (S, I, R) and those who got vaccinated at least once (SV, IV, RV, V). We assume that the group of non-vaccinated susceptible individuals S (and, similarly, infected I and recovered R) is divided into two subgroups: SN and SD. The SN compartment contains such susceptible who will eventually be vaccinated, while those in SD will not.
The SN population is vaccinated with rate v and effectiveness a. Consequently, the individuals from the SN group populate the vaccinated group V with rate av. The individuals in V are considered immune, and we assume that immunization prevents them both from getting infected and infecting others. The SV compartment is composed of S1 and S2 (and, similarly, vaccinated infected IV consists of I1 and I2). Due to vaccine ineffectiveness, people in S1 are perceived as immunized, but in fact are susceptible. S1 is populated from SN with rate (1 − a)v. The vaccinated from the V group move to the S2 group of susceptibles with immunity waning rate ω. The individuals from the S1 group move to S2 with the same rate ω. The S2 group is the group of vaccinated, but no longer immune, and thus, susceptible individuals. In contrast to S1 , we consider that the S2 group is subject to revaccination. Consequently, a fraction of size a of the population from S2 populates V with rate avr and a fraction of size (1 − a) populates S1 with rate (1 − a)vr. Across the manuscript, we assume vr = v, but the model is general and different values can be considered. Individuals from S1 move to S2 with rate ω to ensure that the ineffectively vaccinated are revaccinated with the same speed as the ones for which the vaccine was effective.
Some of the susceptibles in S1 (or, similarly, S2) may not get revaccinated fast enough and may become infected and populate I1 (or, I2). Then, as in the classical SIRS model, the I1 (or I2) population recovers and populates group RV with rate γ. We consider that the recovered in RV may also lose the immunity, and become susceptible again and move to S2 with rate κ. The remaining susceptible subgroups (the SN and SD) may undergo the same classical dynamics, i.e., become infected, recover, and either become susceptible again or, in case of the recovered in the RN subgroup, become vaccinated with rate v.
The following parameters are used to describe population dynamics in the model:
fv,f : restrictions level (for VP holders and others)
β0 : basic transmission rate
βv, β : transmission rate (for VP holders and others)
γ: recovery rate
κ: natural immunity waning rate
a: vaccination effectiveness
v: vaccination rate
vr: revaccination rate
ω: vaccine-induced immunity waning rate
d : fraction of population that will never get vaccinated
Finally, the following set of ordinary differential equations (ODEs) defines the dynamics where also the following relations hold with the constraint S, SV, I, TV, R, RV ≥ 0. Finally, to consider the subpopulation dynamics in terms of fractions of the entire subpopulation, we set and denote d to be the fraction of the never-vaccinated population
Modeling restrictions
We assume that the VP holders consist of the following subpopulations of vaccinated at least once: V, SV ,IV, RV. Recall that the net effect of all non-pharmaceutical interventions is modeled using parameters fv and f, called restrictions throughout the text. The parameter fv amounts to the level of restriction of contacts, and thus the ability to infect, within the group of VP holders. The parameter f satisfies f > fv and corresponds to restriction of contacts within the rest of the population, as well as between the VP holders and the rest of the population.
The restriction level fv for the VP holders is introduced in the model as a modulator of the transmission rate βv. Specifically, we assume that βv = β0(1 − fv), where β0 is the transmission rate of the SARS-CoV-2 virus without restrictions. We assume fv ranges from 0 to 1, where fv = 0 corresponds to no restrictions enforced on the VP holders, and fv = 1 corresponding to full restrictions. Given that for fv = 0 the reproduction number Rmax = β0/γ , and that the recovery rate γ = 1/6, we obtain the no-restriction transmission rate β0= Rmax/6. Thus, for the Delta variant, with Rmax= 6, β0= 1
Similarly, the transmission rate parameter β = β0(1 − f) describes the transmission rate within the rest of the population and between VP holders and the rest, given the restrictions f.
Proportional versus preferential types of social mixing
The above described model equations are based on the assumption that the social mixing between social groups in the population is proportional to the group sizes (the mass action principle). Instead, preferential mixing can be assumed, where the VP holders are more likely to contact other VP holders, since they have lower restrictions [53]. This preferential bias is proportional to the difference between the restrictions f and fv. To incorporate the preferential mixing effect in the ODE model (Equation 1) we rescale the interaction terms according to the following rules: where S + I + R is the non-immune population.
Model simulations
For simulations, we solve the model numerically by means of joint Adams’ and BDF methods, as implemented in the R package deSolve, lsoda method of the ode function [54]. The method monitors data in order to select between non-stiff (Adams’) and stiff (BDF) methods. It uses the non-stiff method initially [55].
To generate the data presented in Figure 1b, we use the reference setup of parameters for the Delta variant: β0 = 1, f = 0.77 (and thus β = 0.23), fv = 0.55 (and thus βv = 0.45), γ = 1/6, K = 1/500, a = 0.79, v = vr= 1/250, ω = 1/500, d = 0.12, with initial conditions I = 10−6, ID= d · I = 10−7; IN = (1 − d) · I = 0.9 · 10−6, R = 0, V = 0. Given I(t) resulting from the solution of the model’s ODE system, to present the final results as easier interpretable cases per million rather than fractions, we re-scale the results by 1M. Additionally, we compute a proxy for the daily incidence number of new cases from the following relation between I (t) and I * (t): Thus, the I* (t) is computed as We proceed similarly to obtain daily incidence numbers and for the sum of all infected, and again to make it interpretable in the figures we re-scale it by 1M.
Stability analysis
The vaccination dynamics can be solved explicitly in the absence of infections. Fixing I = IV = R = RV = 0, and assuming v = vr, we obtain For convenience, where it is not needed, we drop the time argument.
Taking an adiabatic approach we linearize the infection dynamics for small I, IV and R under the assumption of slowly varying S, SV and V. In that case, the infection dynamics decouples from the vaccination dynamics and the Jacobian submatrix Jsub for the equations for I and IV is given by: Given the Jacobian submatrix, we can approximate the dynamics in a small neighborhood of the I = IV = 0 state as
The instantaneous reproduction number R* and the instantaneous doubling time D
Since the largest and the second largest eigenvalues λmax and λ2 of Jsub are both real, the solution to Equation 3 providing the dynamics of infection numbers of the vaccinated and the rest of the population in time can be written in the following form where w1 and w2 are the respective eigenvectors, and c1 and c2 are constants depending on the initial conditions.
Since we have λ2 − λmax≤ 0, we can approximate the time evolution of infection numbers by The largest eigenvalue of Jsub is given by whereby it is convenient to express λmax as a function of and We then obtain Given the population fractions S (t) and SV (t) at a given time instant t, the linearized dynamics of infections given by Equation 3 has a corresponding two-type Galton-Watson branching process, which is a microscopic description of the dynamics. The two types of the process correspond to the I and IV groups. The type I individuals generate Pois (R1S) offsprings of type I and Pois (R1SV) offsprings of type IV. The type IV individuals generate Pois (R1S) offsprings of type I and Pois (R2SV) off-springs of type IV. The linearized dynamics (3) can then be understood as a mean field limit of the microdynamics described by such a branching process. Moreover, the spectral norm of the transition matrix of the branching process can be interpreted as the reproduction number of the branching process, since the expected number of infected in generation n grows like const · (R*)n [56]. We refer to R* as the instantaneous reproduction number. The term instantaneous comes from the fact that we are considering the linearized adiabatic dynamics in a small neighborhood of the I = IV = 0 (ref Eq. 3).
The above discrete branching process can be extended to a continuous time branching process by assuming a probability distribution on the generation time, denoted φ (γ). The growth of the continuous time branching process const · eαt is characterized by its Malthusian growth parameter, denoted a. The relation between the instantaneous reproduction number R*, the distribution φ(τ) and the Malthusian parameter α for such a branching process is given by where ℒφ(α)is the Laplace transform of the distribution φ [56]. Since the setting of ODE model (1) implies exponential distribution of the generation time, i.e, φ(γ) = Exp(γ), the following relation holds: α = γ(R* − 1).
By Equation 5, the Malthusian parameter a for our dynamics is given by the largest eigenvalue λmax. Hence we obtain the relation between the instantaneous reproduction R* and the λmax as λmax = γ(R* − 1). Note that since both S and SV are functions of time, so are λmax and R*.
It is noteworthy that in the above equations, all R1, R2, R1S and R2SV, and R* should be seen as reproduction numbers, but of a different nature [57]. R1 and R2 are reproduction numbers taking into account the restrictions f and fv, respectively. The R1S and R2SV are also group specific, but in addition incorporate the respective group sizes. Finally, R* combines all these factors together.
Having this and Equation 5, we define the instantaneous doubling time at time, denoted t D(t), as the solution D of . Such obtained doubling times are featured in Supplementary Figure S1.
The times of transitions between subcritical and overcritical epidemics
The analysis of the linearized dynamics around I = IV = 0 allows us to determine transitions between subcritical and overcritical epidemics. Such transitions occur at the time instants t at which λmax(t) = 0, or, equivalently, at R*(t) = 1. We thus find that for given values of S(t) and SV (t) the critical times t for transitions between subcritical and overcritical epidemics are the roots of the equation The obtained critical threshold times are plotted in the lower triangles of the panels in Figures 2 and 3 in the main text. In the case of proportional mixing the above equation is equivalent to:
Asymptotic structure of the population
The asymptotic structure of the population in terms of the sizes of the subpopulations V, SV and SD can be easily obtained by setting I = IV = R = RV = 0 and computing the stable stationary solution for Vas, Sas and of our ODE system 1: where can be seen as the actual immunization rate in the population, and is expressed as a function of vaccine effectiveness a and the ratio of the immunity waning rate ω and the revaccination rate vr. The obtained values correspond to the structure in the limit t → ∞ and represent the structure to which the population converges in the long term.
Having this, we obtain the asymptotic instantaneous reproduction number R* by inserting the asymptotic values Sas and into Equation 8. These values are plotted in the upper triangles in the panels of Figures 2 and 3 in the main text.
Finally, we solve for such minimum common restrictions f = fv = fmin, which will result in instantaneous reproduction number R* = 1 for the different vaccine effectiveness and vaccination rate setups. Hence fmin is found from as
Endemic state
The endemic state of the VAP-SIRS model is obtained by setting the derivatives of the ODE system 1 to 0. A straightforward computation reduces the endemic system of equations to the following: where for the proportional mixing we have: δ* = δ and δ+ = β, and for the preferential mixing we set and . The roots of the Equations 9 and 10 are plotted in Figures 4 and 5.
Data Availability
Not applicable
https://github.com/storaged/VAP-SIRS
Data and materials availability
The VAP-SIRS model was implemented using R version 4.0.2 along with the shiny package to build an interactive web application that allows to simulate the model. The code of the model is available online in the GitHub repository: https://github.com/storaged/VAP-SIRS, and the on-line tool is available http://bioputer.mimuw.edu.pl:85/VAP-SIRS/. The code to generate Figures 2-5 from the main text is available at https://github.com/eMaerthin/VAP_SIRS_Ananysis.
Authors contributions
AG, KG, TK, and ES conceived the VAP-SIRS model - with input and feed-back on the model and results from TC, GG, MP, EP and MR. TK performed the stability analysis. KG implemented model simulations and the Shiny application for visualizations. MB implemented the stability analysis. SC, TC, EP, and MR performed literature search. ES supervised the study. All authors wrote and provided critical feedback to the manuscript.
Competing interests
Other projects in the research lab of ES are co-funded by Merck Healthcare KGaA.
Acknowledgments
SC acknowledges support by University of Malta. EP acknowledges support by the University of Crete. TC has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101016233 (PERISCOPE). GG acknowledges support by the University of Trento within the COVID-19 Strategic Project MOSES (Models for Reasoning about the Spreading of Diseases). MP was supported by the Slovenian Research Agency (Grant Nos. P1-0403 and J1-2457). ES acknowledges funding by the Polish National Science Centre OPUS grant no 2019/33/B/NZ2/00956.
Footnotes
↵† Shared first authorship,
The entire Manuscript has been updated to the current situation when the Delta variant is becoming dominant across Europe.