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.


Plain language summary
Assessing the impact of vaccines, other public health measures, and declining immunity on SARS-CoV-2 control is challenging. This is particularly true in the context of vaccination passes, whereby vaccinated individuals have more freedom of making contacts than unvaccinated ones. Here, we use a mathematical model to simulate various scenarios and investigate the likelihood of containing COVID-19 outbreaks in example European countries. We demonstrate that both Alpha and Delta SARS-CoV-2 variants inevitably lead to recurring outbreaks when measures are lifted for vaccination pass holders. High revaccination rates and a lowered fraction of the unvaccinated population increase the benefit of vaccination passes. These observations are important for policy making, highlighting the need for continued vigilance, even where the epidemic is under control, especially when new variants of concern emerge. COMMUNICATIONS MEDICINE | (2022) 2:23 | https://doi.org/10.1038/s43856-022-00084-w | www.nature.com/commsmed I n 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][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 wellbeing, 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][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 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][25][26][27][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][33][34][35][36][37][38][39][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 postvaccination 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 mediumand 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 longterm. 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 V , I V , R V , 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 N and S D . The S N compartment contains such susceptible who will eventually be vaccinated, while those in S D will not.
The S N population is vaccinated with rate υ and effectiveness a. Consequently, the individuals from the S N 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 V compartment is composed of S 1 and S 2 (and, similarly, vaccinated infected I V consists of I 1 and I 2 ). Due to vaccine ineffectiveness, people in S 1 are perceived as immunized, but in fact are susceptible. S 1 is populated from S N with rate (1 − a)υ. The vaccinated from the V group move to the S 2 group of susceptibles with immunity waning rate ω. The individuals from the S 1 group move to S 2 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 S 2 group is the group of vaccinated, but no longer immune, and thus, susceptible individuals. In contrast to S 1 , we consider that the S 2 group is subject to revaccination. Consequently, a fraction of size a of the population from S 2 populates V with rate aυ r and a fraction of size (1 − a) populates S 1 with rate (1 − a)υ r . Across the manuscript, we assume υ r = υ, but the model is general and different values can be considered.
Some of the susceptibles in S 1 (or, similarly, S 2 ) may not get revaccinated fast enough and may become infected and populate I 1 (or, I 2 ). Then, as in the classical SIRS model, the I 1 (or I 2 ) population recovers and populates group R V with rate γ. We consider that the recovered in R V may also lose the immunity, and become susceptible again and move to S 2 with rate κ. The remaining susceptible subgroups (the S N and S D ) 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 N 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: 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 (υ = υ r = 0.004; typical for many European countries), slow immunity waning ω = 0.002, low fraction of nevervaccinated (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 f v = 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 (I Ã 1 ), and vaccinated who already lost immunity (I Ã 2 ). By I Ã Σ (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 (S 1 , light orange) vaccinated susceptible who already lost immunity (S 2 , 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. Finally, the following set of ordinary differential equations (ODEs) defines the dynamics: where also the following relations hold with the constraint S, S V , I, I V , R, R V ≥ 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 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 V , I V , R V . Recall that the net effect of all non-pharmaceutical interventions is modeled using parameters f v and f, called restrictions throughout the text. The parameter f v 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 ≥ f v 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 f v for the VP holders is introduced in the model as a modulator of the transmission rate β v . Specifically, we assume that β v = β 0 (1 − f v ), where β 0 is the transmission rate of the SARS-CoV-2 virus without restrictions. We assume f v ranges from 0 to 1, where f v = 0 corresponds to no restrictions enforced on the VP holders, and f v = 1 corresponding to full restrictions.
Given that for f v = 0 the reproduction number R max ¼ β 0 =γ, and that the recovery rate γ = 1/6, we obtain the no-restriction transmission rate β 0 ¼ R max =6. Thus, for the Delta variant, with R max ¼ 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 45 . This preferential bias is proportional to the difference between the restrictions f and f v . 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: 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: 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 I Ã 1 , I Ã 2 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 V = R = R V = 0, and assuming υ = υ r , 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, I V and R under the assumption of slowly varying S, S V and V. In that case, the infection dynamics decouples from the vaccination dynamics and the Jacobian submatrix J sub for the equations for I and I V is given by: Given the Jacobian submatrix, we can approximate the dynamics in a small neighborhood of the I = I V = 0 state as The instantaneous reproduction number R Ã and the instantaneous doubling time D. Since both the eigenvalues λ max and λ 2 ≤ λ max of J sub 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 form where w 1 and w 2 are the respective eigenvectors, and c 1 and c 2 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 J sub is given by whereby it is convenient to express λ max as a function of 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 V (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 V groups. The type I individuals generate Pois R 1 S À Á offsprings of type I and Pois R 1 S V À Á offsprings of type I V . The type I V individuals generate Pois R 1 S À Á offsprings of type I and Pois R 2 S V À Á offsprings of type I V . 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 49 . 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 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 α. The relation between the instantaneous reproduction number R Ã , the distribution φ τ ð Þ and the Malthusian parameter α for such a branching process is given by Since the setting of ODE model (1) implies exponential distribution of the generation time, i.e, φðγÞ ¼ Exp γ À Á , the following relation holds: α ¼ γ R Ã À 1 ð Þ : By Eq. (4), the Malthusian parameter α 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 S V are functions of time, so are λ max and R Ã .
It is noteworthy that in the above equations, all R 1 , R 2 , R 1 S and R 2 S V , and R Ã should be seen as reproduction numbers, but of a different nature 50 . R 1 and R 2 are reproduction numbers taking into account the restrictions f and f v , respectively. The R 1 S and R 2 S V are also group specific, but in addition incorporate the respective group sizes. Finally, R Ã 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 e γ R Ã ðtÞÀ1 ð Þ Á D ¼ 2. 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 R Ã , we would have for the Malthusian growth parameter α the relation ðα þ cÞðα þ uÞ ¼ ucR Ã , 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 V = 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 S V (t) the critical times t for transitions between subcritical and overcritical epidemics are the roots of the equation λ max ðtÞ ¼ 0: The 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: 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 V and S D can be easily obtained by setting I = I V = R = R V = 0 and computing the stable stationary solution for V as , S as and S as V of our ODE system (1): 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 υ r . 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 S as and S as V 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 f ¼ f v ¼ f min , which will result in instantaneous reproduction number R Ã ¼ 1 for the different vaccine effectiveness and vaccination rate setups. Hence f min is found from Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Results
The VAP-SIRS model of the impact of COVID-19 VPs. The proposed VAP-SIRS model extends the classical SIRS model 44 (red arrows in Fig. 1a) with additional states and parameters that describe the dynamics of vaccination rollout in a population (green arrows in Fig. 1a). To this end, we consider the following subpopulations: (i) initially susceptible S N , who, if successfully vaccinated, populate the immune group V, with rate aυ, where υ is the vaccination rate and a is the vaccination effectiveness, (ii) susceptible who were vaccinated but did not gain immunity (S 1 ), (iii) vaccinated, whose immunity waned with rate ω and who became susceptible again (S 2 ), (iv) susceptible, who are not and will never get vaccinated (S D ). The S D 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 S 2 populates V with rate aυ r . All recovered, unless in the recovered compartment R D , are also subject to vaccination. Before the recovered in the R V lose immunity, they might be revaccinated, and, thus, populate the V group with rate υ r (similarly, R N are vaccinated with rate υ). 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, υ r = υ. Additional sub-compartments are used to model hospitalizations (see the Supplementary Note 1). 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 51,52 . 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 (BioN-Tech/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 nevervaccinated 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 duration of immunity against infection, reflecting emerging data on large individual variation of immunity waning and other key factors influencing this process 15,18,[24][25][26] . 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 15,53,54 . 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) 10, 55 .
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 refs. 56,57 ), while for the B.1.1.7 variant (Alpha) an optimistic estimate is equal to 4 10,55 . Two levels of restrictions are considered: restrictions f v 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 f v < 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 f v 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 f v = 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 unfavorable 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 υ = υ r = 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 f v = 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 (S N ) and suppressed by ongoing vaccination, as expected. Interestingly, the second, larger wave is driven by the S V group. The S V group is composed of the number of individuals for whom the vaccine was ineffective (S 1 ) and those vaccinated who lose their immunity and are not yet revaccinated (S 2 ).
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 linearizing 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, f v ) and the composition of the population (the population fractions of the groups: susceptible with VPs S V , susceptible without VPs S, and the immune group V), 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 Fig. 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 = f v = 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 = f v = 0.77, which results in a time evolution of 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 − f v parameter space, where f v ≤ 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.79the Comirnaty effectiveness on the Delta variant, υ = υ r = 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 − f v 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 f v 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 − f v 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 theregions, 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 Fig. 2 to Supplementary Fig. S2 shows how the Delta variant worsens all scenarios with respect to the Alpha variant: in all panels of Fig. 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 ( Supplementary  Fig. S3) and the Alpha variant ( Supplementary Fig. S4). 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 ( Supplementary Fig. S3a) 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 ( Supplementary  Fig. S3e). 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 (Supplementary Fig. S3c). Finally, combinations of an increased never-vaccinated fraction with other parameter changes show an interesting mix of effects. When both the nevervaccinated 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 ( Supplementary Fig. S3d). Similarly, there is a synergistic effect of the combination of the increased nevervaccinated fraction and the increased immunity waning rate ( Supplementary Fig. S3f). 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. S4).
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 analyzed 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 Fig. 2 and Supplementary Fig. S3 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 in Fig. 2 and Supplementary Figs. S2, S3 and S4, 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 f min for the whole population that would guarantee to avoid an overcritical epidemic in the long-term (for time approaching infinity, Methods): where V as 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 υ r , the fraction of nevervaccinated population d and immunity waning rate ω ( Table 1). The minimum common restrictions for the reference setup are equal to f min ¼ 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 υ r = 0.008, low d = 0.12, low ω = 0.002; Table 1 third row) we obtain V as = 0.6, and f min ¼ 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 Table 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 The studied parameters are: vaccine effectiveness a, revaccination rate υ r , fraction of nevervaccinated 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).
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 f min ¼ 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 and hospitalization numbers. For a given restriction combination f, the above analyzed asymptotic instantaneous reproduction number R Ã (Fig. 2 and Supplementary Figs. S2, S3 and S4) 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 or hospitalization numbers the system converges to. To this end, we first 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. 3) and compare it to the scenarios achieved with the Alpha variant ( Supplementary  Fig. S5). In contrast to the computation of the instantaneous reproduction number R Ã and its asymptotic values, which is 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 58 . 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. 3c). 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 Supplementary Fig. S6; the effect of simultaneous parameter variations is akin to that described earlier for the values of asymptotic R Ã and time to critical thresholds in Supplementary Fig. S3.
Again, comparison with the endemic infection numbers predicted for the Alpha variant (Supplementary Figs. S5 and S7) 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.
Besides the computation of daily infection numbers, we compute daily hospitalization numbers for the Delta (Fig. 4) and Alpha (Supplementary Fig. S8) variants in the endemic state (see the Supplementary Note 1 for details). The benefit of vaccination in reducing hospitalizations is striking: for all regions apart from the (-) region, the number of hospitalized cases in the unvaccinated population is roughly an order of magnitude larger than in the vaccinated population. In addition, the unvaccinated population is even more prone to hospitalization when the fraction of never-vaccinated is larger ( Fig. 4c and Supplementary  Fig. S8c). This negative effect of the never-vaccinated population fraction is dominant also when double changes of parameter values are considered with respect to the reference setup, both for the Delta (Supplementary Fig. S9) and the Alpha variant ( Supplementary Fig. S10). Again, the comparison between the Delta and the Alpha variant highlights the deleterious effect of increased transmissibility of the Delta variant. Clearly, to keep hospitalizations at the same level as for the Alpha variant, the Delta variant requires much stricter restrictions for the VP holders and the rest of the population.
The benefit of VPs is larger for the Alpha than for the Delta variant, and strongly depends on the (re-)vaccination rates. The above analysis demonstrates the potential risks of increased infection and hospitalization numbers. The additional freedom for VP holders can be considered beneficial for society as long as it does not lead to an uncontrolled surge of infections. Consequently, considering the relevant parameters as f v ≤ f, we estimate the benefit of VPs as a value in [0, 1] given by the fraction of the relevant f − f v parameter space where the asymptotic instantaneous reproduction number satisfies R Ã ≤ 1. For example, for the Delta variant and the reference setup, the benefit of VPs can be obtained as a fraction of the corresponding upper triangle in Fig. 2, which is colored in white or shades of blue.
The estimated benefit of VPs for the reference parameter setup and the Alpha variant is almost 0.5 (Fig. 5). Increasing the (re-) vaccination rates has the most positive and the largest impact on the VP benefit, increasing it to around 0.8 for the Alpha variant. All other considered changes to the reference setup (decreasing the effectiveness of the vaccine, preferential instead of proportional mixing, increased waning rate ω, increasing the fraction of never-vaccinated d) decrease the VP benefit. Strikingly, for half of the considered parameter setups (the reference setup and setups with single changes to it), the benefit is over two-fold reduced for the Delta variant, as compared to the Alpha variant. The maximum estimated benefit for the Delta variant is only around 0.23.

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 effectiveness. 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 consider mortality or intensive care unit utilization like some other models, albeit in the context of exploring different parameters than larger freedom for VP holders 32,34,[36][37][38][39] . In this context, the advantage of our model is that it is enriched in features such as revaccinations and waning immunity, which have been shown to be very relevant not only in the long term 22,27 . 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. Possible extensions to our model could include inter-individual variations in immunity, which would render it relevant for people with immunodeficiencies. The presented analysis has been performed assuming that without restrictions, the maximum reproduction number of the virus is R max ¼ 6 or R max ¼ 4 for Delta and Alpha variants, respectively. More transmissible variants could easily be modeled by fixing higher values of R max . 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. Clearly, the assumption of a constant vaccination rate is-although common in the literature-a simplification of reality. For instance, booster shots are usually recommended only after a minimal time span of six months after the last vaccine dose, in order to maximize the effect on affinity maturation and immunological memory. Instead, our assumptions of a fixed vaccination rate corresponds to a situation where boosting shots could be taken earlier than this minimal time span, as is implemented in some countries facing another wave. Such minimal time delay could be accounted for using a delay equation, but we do not expect that such modification would affect the conclusions from our analysis. The presented analysis assumes that the initial number of infected and the number of recovered are both negligible. Given the large heterogeneity of infection dynamics 59 , it is impossible to assign initial numbers that would be representative for all countries. We thus choose a universal reference assuming these numbers are very small or zero. The Shiny app for VAP-SIRS simulations, however, allows initializing with different values. For example, since the recovered, similarly to vaccinated, are also immunized for some time, the main conclusion when increasing the initial recovered number to some level is that it moves the dynamics forward in time in our analysis, as if the vaccinated sub-population started at this level.
Despite limitations, our model accounts for key parameters influencing infection dynamics and gives valuable insights into policies pertaining to the introduction of VPs, rendering 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. 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 Fig. 5 The estimated benefit of vaccination passes (VPs). For each parameter setting (reference setup and single or double parameter value changes to it; y axis), the benefit of VPs was evaluated (x axis) for two SARS-CoV-2 strains: Alpha (green) and Delta (red). 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. Two illustrative findings are that: (i) the minimum level of common restrictions to avoid resurgence in the reference setup has doubled from 0.29 (Alpha) to 0.62 (Delta), and that (ii) the largest VP benefit has decreased almost four-fold from around 0.8 for Alpha to 0.23 for 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, numbers of cases per million in the endemic state, numbers of hospitalizations in the endemic state, and VP benefit 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 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, unfavorable 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 utilization of tools such as the VAP-SIRS model, along with different tools available to policymakers should be explored in the context of monitoring the implementation of VPs, including the EU DCC measures, to ensure optimization of key parameters. In this manner, evidence-informed policy making would be safeguarded as would the best possible outcomes in terms of effectively combating the current pandemic.

Data availability
All source data underlying the graphs and charts presented in the main and supplementary figures are publicly available on-line in the GitHub repository: https:// github.com/eMaerthin/VAP_SIRS_Analysis as text (.csv) files.

Code 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 on-line in the GitHub repository: https://github.com/ storaged/VAP-SIRS 60