Risk assessment of COVID-19 epidemic resurgence in relation to SARS-CoV-2 variants and vaccination passes.

The introduction of COVID-19 vaccination passes (VPs) by many countries coincided 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, only arose 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, received unprecedented attention[2-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 was estimated to be 50% more transmissible than the Alpha variant (B.1.1.7), already estimated to be 50% more transmissible than the previously dominant strains[9-11].

Evidence indicates vaccine effectiveness can greatly vary[12,13] and it may be compromised due to escape variants[14] and waning immunity[15-18]. Preliminary data from several countries indicate reduced vaccine effectiveness against the infection with the Delta variant compared to the Alpha variant[19-21], even as low as 64% for the Comirnaty (Pfizer-BioNTech) vaccine according to data from Israel[22]. Emerging evidence suggests that the vaccines are effective in preventing serious illness and hospitalization[11,20,21].

Still, avoiding another COVID-19 infection resurgence remains a valid and potentially attainable goal[23]. Immunity against both infection and hospitalization wanes over time[15,18,24-28]. An estimated 10% of COVID-19 infections will have long-term sequelae (long COVID), posing an increasing threat to national health systems[29,30]. 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[31]. 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 against infection, 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 were developed to inform vaccination strategies[32-40]. 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[41,42]. Scenarios for the post-vaccination era were also considered by Sandmann and colleagues (2021), finding that under realistic scenarios periodic epidemics are likely[43]. 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 utilizing the VPs for different purposes at national level by establishing different levels of freedom for VP holders in terms of accessing premises, facilities, traveling within a country, etc., it is important to examine the broad parameters determining how to optimize the implementation of measures such as the EU DCC and other VPs.

To address these needs, we propose a mathematical model called VAP-SIRS, which accounts for key parameters that impact the effective reproduction number of the virus, and consequently, infection dynamics: vaccination effectiveness, rates of (re-)vaccination and waning immunity, and the differences between SARS-CoV-2 variants. We perform comprehensive analysis for different levels of restrictions for VP holders and the rest of the population, for various realistic setups of these key parameters, including the different effectivenesses of the Comirnaty and Vaxzevira vaccines on the Delta and Alpha variants, as well as fractions of never vaccinated in the United Kingdom and France. 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 policy making. VAP-SIRS predicts a possible infection resurgence despite vaccinations. The resurgence is due to the lowered levels of restrictions for the VP holders compared to the rest of the population, while for some fraction of those VP holders the vaccine was ineffective and for the others the immunity may wane before they become re-vaccinated. A thorough analysis of our model identifies the complete set of potential scenarios for the COVID-19 epidemic depending on the restrictions imposed on VP holders and the rest of the population. For these scenarios, we estimate daily infection as well as hospitalization numbers and identify flexible measures to avoid epidemic resurgence. In particular, we derive the minimum common restriction level for the VP holders and the rest of the population, which can keep the epidemic subcritical in the long-term. Finally, we estimate the social benefit of VPs and find its strong dependence on (re-)vaccination rates.

Methods

Mathematical model

We introduce VAP-SIRS (VAccination Passes in Susceptible-Infectious-Recovered-Susceptible model), as an extension to the classical SIRS model[44] (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 (S, I, R, V). We assume that the group of non-vaccinated susceptible individuals S (and, similarly, infected I and recovered R) is divided into two subgroups: S and S. The S compartment contains such susceptible who will eventually be vaccinated, while those in S will not.

Fig. 1

The VAP-SIRS model and its predicted scenarios.

a Graphical scheme of the VAP-SIRS model. b, c Predicted scenarios for the reference setup for the Delta variant, with vaccine effectiveness a = 0.79 (corresponding to the effectiveness of the Comirnaty vaccine against infection with the Delta variant), slow (re-)vaccination rate (υ = υ = 0.004; typical for many European countries), slow immunity waning ω = 0.002, low fraction of never-vaccinated (d = 0.12; corresponding to the fraction in the United Kingdom) and proportional mixing (see Methods). b Color curves: Timeline of daily incidence per 1 million inhabitants in different infected compartments for the combination of restrictions f = 0.77 and fv = 0.55. A variable with the asterisk (*) indicates that we consider a daily incidence over the corresponding variable. The dashed lines describe infected who are: non-vaccinated (I*, yellow), vaccinated who did not gained immunity (), and vaccinated who already lost immunity (). By (red, solid line) we mean the sum of all daily infected (). Color bands: Muller plot of the population structure (the width of the color band in the y axis) as a function of time (x axis) for the same parameter settings. Colors correspond to specific subpopulations: non-vaccinated susceptible (S, yellow), vaccinated susceptible who did not gained immunity (S1, light orange) vaccinated susceptible who already lost immunity (S2, dark orange), vaccinated immunized (V, green). Moreover, by IΣ (red) and RΣ (blue) we denote all infected and all recovered (independently of vaccination result), respectively. c Time evolution of the instantaneous reproduction number  (y axis) depending on the number of days counted from the start of the vaccination program (x axis), in five different scenarios describing the epidemic evolution: overcritical (+, red, f = 0.77 and fv = 0.38), subcritical (-, blue, f = 0.92 and fv = 0.71), initially and eventually overcritical (+ - +, orange, the same restrictions as in b: f = 0.77 and fv = 0.55), eventually overcritical (-+, pink, f = 0.92 and fv = 0.38), and eventually subcritical (+-, cyan, with f = 0.77 and fv = 0.71). As controls, two additional scenarios of the epidemic evolution are presented, corresponding to no implementation of VPs and no changes in behavior due to vaccination: subcritical (another example of - scenario, green) with f = f = 0.92 and eventually subcritical (another example of +- scenario, yellow) with f = f = 0.77, both plotted with dot-dashed line.

The S population is vaccinated with rate υ and effectiveness a. Consequently, the individuals from the S group populate the vaccinated group V with rate aυ. The individuals in V are considered immune, and we assume that immunization prevents them both from getting infected and infecting others. The S compartment is composed of S1 and S2 (and, similarly, vaccinated infected I 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 S with rate (1 − a)υ. 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 ω to ensure that the ineffectively vaccinated are revaccinated with the same speed as the ones for which the vaccine was effective. 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 aυ and a fraction of size (1 − a) populates S1 with rate (1 − a)υ. Across the manuscript, we assume υ = υ, but the model is general and different values can be considered.

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 R with rate γ. We consider that the recovered in R may also lose the immunity, and become susceptible again and move to S2 with rate κ. The remaining susceptible subgroups (the S and S) may undergo the same classical dynamics, i.e., become infected, recover, and either become susceptible again or, in case of the recovered in the R subgroup, become vaccinated with rate υ.

The following parameters are used to describe population dynamics in the model:where the transmission rates are expressed int terms of the basic transmission rate and the restriction level parameters:

Finally, the following set of ordinary differential equations (ODEs) defines the dynamics:where also the following relations holdwith the constraint S, S, I, I, R, R ≥ 0. Finally, to consider the subpopulation dynamics in terms of fractions of the entire subpopulation, we setand denote d to be the fraction of the never-vaccinated population

The endemic state of the VAP-SIRS model is computed in the Supplementary Note 1.

Modeling restrictions

We assume that the VP holders consist of the following subpopulations of vaccinated at least once: V, S, I, R. 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 , and that the recovery rate γ = 1/6, we obtain the no-restriction transmission rate . Thus, for the Delta variant, with , β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[45]. This preferential bias is proportional to the difference between the restrictions f and fv. Preferential mixing is a common, socio-psychological motivated mixing scheme alternative to proportional mixing. In this scheme the group interaction is still proportional, but biased by the relative degree of freedom given to the passport holders. Preferential mixing as a modulation of proportional mixing was previously studied in the context of infectious diseases by Glasser et al.[46]. To incorporate the preferential mixing effect in the ODE model (Eq. (1)) we rescale the interaction terms according to the following rules:where S + I + R is the non-immune population.

Numerical integration and parameter values

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[47]. The method monitors data in order to select between non-stiff (Adams’) and stiff (BDF) methods. It uses the non-stiff method initially[48].

To generate the data presented in Fig. 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, κ = 1/500, a = 0.79, υ = υ = 1/250, ω = 1/500, d = 0.12, with initial conditions I = 10−6, I = d ⋅ I = 10−7; I = (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 and :Thus, the I*(t) is computed asWe 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 = I = R = R = 0, and assuming υ = υ, we obtainFor convenience, where it is not needed, we drop the time argument.

Taking an adiabatic approach we linearize the infection dynamics for small I, I and R under the assumption of slowly varying S, S and V. In that case, the infection dynamics decouples from the vaccination dynamics and the Jacobian submatrix J for the equations for I and I is given by:Given the Jacobian submatrix, we can approximate the dynamics in a small neighborhood of the I = I = 0 state as

The instantaneous reproduction number and the instantaneous doubling time D

Since both the eigenvalues λ and λ2 ≤ λ of J are real, the solution to Eq. (3) providing the dynamics of infection numbers of the vaccinated and the rest of the population in time can be written in the following formwhere w1 and w2 are the respective eigenvectors, and c1 and c2 are constants depending on the initial conditions.

Since we have , we can approximate the time evolution of infection numbers by

The largest eigenvalue of J is given bywhereby it is convenient to express as a function of and . We then obtain

We now describe the relation of the analyzed system with the corresponding branching process, which motivates the notion of the instantaneous reproduction number and the derivation of the doubling time. It also allows a straightforward generalization to more complex systems of equations than the one considered here. Given the population fractions S(t) and S(t) at a given time instant t, the linearized dynamics of infections given by Eq. (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 I groups. The type I individuals generate offsprings of type I and offsprings of type I. The type I individuals generate offsprings of type I and offsprings of type I. 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 normof the transition matrixof 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 [49]. We refer to 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 = I = 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 is characterized by its Malthusian growth parameter, denoted α. The relation between the instantaneous reproduction number , the distribution and the Malthusian parameter α for such a branching process is given bywhere is the Laplace transform of the distribution φ[49]. Since the setting of ODE model (1) implies exponential distribution of the generation time, i.e, , the following relation holds:

By Eq. (4), the Malthusian parameter α for our dynamics is given by the largest eigenvalue . Hence we obtain the relation between the instantaneous reproduction and the as Note that since both S and S are functions of time, so are and .

It is noteworthy that in the above equations, all R1, R2, R1S and R2S, and should be seen as reproduction numbers, but of a different nature[50]. R1 and R2 are reproduction numbers taking into account the restrictions f and fv, respectively. The R1S and R2S are also group specific, but in addition incorporate the respective group sizes. Finally, combines all these factors together.

Having this and Eq. (4), we define the instantaneous doubling time at time, denoted tD(t), as the solution D of . Such obtained doubling times are featured in Supplementary Fig. S1.

Only a small change is needed in the derivation to extend to more complex systems than the considered SIRS model. For example, in a dynamics with exponential Exp(u)-distributed additional incubation time u and exponential Exp(c)-distributed duration of the infectious period c (a so called SEIRS model), and a given , we would have for the Malthusian growth parameter α the relation , from which one can easily compute the corresponding doubling time.

The times of transitions between subcritical and overcritical epidemics

The analysis of the linearized dynamics around I = I = 0 allows us to determine transitions between subcritical and overcritical epidemics. Such transitions occur at the time instants t at which , or, equivalently, at . We thus find that for given values of S(t) and S(t) the critical times t for transitions between subcritical and overcritical epidemics are the roots of the equationThe obtained critical threshold times are plotted in the lower triangles of the panels in Fig. 2 and Supplementary Fig. S3 in the main text. In the case of proportional mixing the above equation is equivalent to:

Fig. 2

Possible COVID-19 epidemic dynamics for different parameter setups for the Delta variant.

The relevant f − f parameter space, where f ≤ f, can be divided into five regions (delimited by black borders), each associated with a different behavior of the epidemics. On the diagonal (white dashed line), f = f, i.e., the restrictions for VP holders and for the rest of the population are the same - corresponding to the situation when VPs are not introduced at all. Lower triangles show the time until the last critical threshold: different color scales correspond to the time until the switch either from a subcritical to an overcritical epidemic (time until overcriticality, violet-green scale), or from an overcritical to a subcritical epidemic (time until subcriticality, yellow-pink scale). Upper triangles show the asymptotic , as a function of the values of f and f (blue-red scale, with blue associated with and red associated with ). a Reference setup, with a = 0.79 (corresponding to the effectiveness of the Comirnaty vaccine on the Delta variant), υ = υ = 0.004, ω = 0.002, d = 0.12 (fraction of never-vaccinated in the United Kingdom) and proportional mixing. The choices of (f, f) corresponding to the five scenarios exemplified in Fig. 1c are denoted by points of the same color. b Setup with a decreased vaccine effectiveness: a = 0.6 (corresponding to the effectiveness of the Vaxzevria vaccine on the Delta variant). c Setup with an increased vaccination rate: υ = υ = 0.008. d Setup with preferential (instead of proportional) mixing. e Setup with an increased fraction of people who will not get vaccinated: d = 0.3 (fraction of never-vaccinated in France). f Setup with an increased waning rate: ω = 1/200.

Asymptotic structure of the population and minimum common restrictions required to avoid epidemic resurgence

The asymptotic structure of the population in terms of the sizes of the subpopulations V, S and S can be easily obtained by setting I = I = R = R = 0 and computing the stable stationary solution for Vas, Sas and of our ODE system (1):wherecan 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 υ. 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 by inserting the asymptotic values Sas and into Eq. (6). These values are plotted in the upper triangles in the panels of Fig. 2 and Supplementary Fig. S3 in the main text.

Finally, we solve for such minimum common restrictions , which will result in instantaneous reproduction number for the different vaccine effectiveness and vaccination rate setups. Hence is found from as

Table 1

Asymptotic level of immunization Vas and minimum common restrictions for the Delta variant and different parameter setups.

Parameter setup	a	υr	d	ω	Vas	f_min
Ref. setup	0.79	0.004	0.12	0.002	0.46	0.69
Dec. a	0.6	0.004	0.12	0.002	0.35	0.74
Inc. υr	0.79	0.008	0.12	0.002	0.56	0.62
Inc. d	0.79	0.004	0.3	0.002	0.37	0.74
Inc. ω	0.79	0.004	0.12	0.005	0.31	0.76
Dec. a, inc. υr	0.6	0.008	0.12	0.002	0.42	0.71
Dec. a, inc. d	0.6	0.004	0.3	0.002	0.28	0.77
Dec. a, inc. ω	0.6	0.004	0.12	0.005	0.23	0.78
Inc. υr, inc. d	0.79	0.008	0.3	0.002	0.44	0.70
Inc. υr, inc. ω	0.79	0.008	0.12	0.005	0.43	0.71
Inc. d, inc. ω	0.79	0.004	0.3	0.005	0.18	0.80

The studied parameters are: vaccine effectiveness a, revaccination rate υ, fraction of never-vaccinated d, and waning immunity rate ω. The first row concerns the reference setup; rows below are setups with the same parameters as in the reference setup, but with either one parameter changed (in bold; rows 2–5; same as in Figs. 2 and 3, apart from preferential mixing, as it is not relevant for common restrictions) or two parameters changed (in bold; rows 6–11). Dec. - Decreased, Inc. - increased.