The effect of COVID-19 vaccination in Italy and perspectives for living with the virus

COVID-19 vaccination is allowing a progressive release of restrictions worldwide. Using a mathematical model, we assess the impact of vaccination in Italy since December 27, 2020 and evaluate prospects for societal reopening after emergence of the Delta variant. We estimate that by June 30, 2021, COVID-19 vaccination allowed the resumption of about half of pre-pandemic social contacts. In absence of vaccination, the same number of cases is obtained by resuming only about one third of pre-pandemic contacts, with about 12,100 (95% CI: 6,600-21,000) extra deaths (+27%; 95% CI: 15–47%). Vaccination offset the effect of the Delta variant in summer 2021. The future epidemic trend is surrounded by substantial uncertainty. Should a pediatric vaccine (for ages 5 and older) be licensed and a coverage >90% be achieved in all age classes, a return to pre-pandemic society could be envisioned. Increasing vaccination coverage will allow further reopening even in absence of a pediatric vaccine.


Model for SARS-CoV-2 transmission and vaccination
We developed an age-structured stochastic model for SARS-CoV-2 transmission and vaccination, based on a susceptible-infectious-removed-susceptible scheme (SIRS) and adapted from previously published models [1,2]. Mixing patterns are assumed to be heterogeneous across ages according to an age-specific social contact matrix estimated prior to the COVID-19 pandemic [3]. We assume an age-dependent susceptibility to SARS-CoV-2 infection: lower in children under 15 years of age and higher for the elderly (65+), compared to individuals of working age [4].
We simulate a two-dose vaccination campaign. In the baseline analysis, vaccination is assumed to reduce the individuals' susceptibility to SARS-CoV-2 infection and the risk of death. Breakthrough infections (i.e., infections in vaccinated individuals) are assumed to be half as infectious as those in unvaccinated individuals [5,6]. The model accounts for waning of both natural and vaccine-induced protection, where the duration is exponentially distributed (mean: 2 years [7,8]). Specifically, before waning of natural immunity individuals are fully protected against infection; while before waning of vaccine protection, vaccine reduces the probability of infection and death (with different efficacy estimates for the two endpoints) and the infectiousness of breakthrough infections. After waning of either natural or vaccine protection, individuals are considered fully susceptible. Alternative durations for natural immunity and vaccine protection are explored as sensitivity analyses (1 year and 10 years), as well as the case of a vaccine not reducing infectiousness (see Section 2.1). Immunity gained after experiencing both vaccination and infection (independently of the order) is assumed not to wane over the time scale of our simulations.
The baseline model is described by the following system of differential equations (summarized by the schematic representation in Supplementary Fig. 1 where: • the population class {a,c} represents individuals in age group a and in underlying conditions status c ("with" or "without"); Susceptible individuals are exposed to a time and age-dependent force of infection λ * (t) which is defined as: • is a scaling factor shaping SARS-CoV-2 transmissibility in the absence of physical distancing restrictions (PDRs) and other non-pharmaceutical interventions (NPIs) such as face masks or hand hygiene precautions, computed by assuming R0=3.0, as estimated for historical lineages of SARS-CoV-2 in Italy [9,10,11]. • Alpha is a coefficient representing the transmissibility increase of the Alpha variant with respect to historical lineages. We assumed Alpha = 0.5 [12][13][14][15][16].
• (1 − ) is a coefficient representing the reduction in the per-contact transmission probability ascribable to preventive transmission measures such as face masks or hand hygiene precautions. We assumed = 20% [9,10].
1} is a scaling factor representing the proportion of pre-pandemic contacts that are active at time .

•
! is the relative susceptibility to SARS-CoV-2 infection at age a. We sample susceptibility profiles from the posterior distribution estimated in [4] • !,! 0 represents the age-group-specific contact matrix, as estimated before the SARS-CoV-2 pandemic [3], whose entries describe the mean numbers of persons in age group D encountered by an individual of age group in an average day.
• ! 0,# represents the number of individuals in the population class { D,c}. For all infectious compartments, the average duration of infectiousness (1/ ) is set equal to the average generation time (6.6 days) [9]. The model described by the system of ordinary differential equations above is implemented through a stochastic discrete-time model with a time step = $ ' . At each time t, the number of unvaccinated individuals in the population class {a,c} (red compartments in Supplementary Fig. 1) who will receive a first vaccine dose is determined as a fraction !,# ( ) of the corresponding population: where !,# ( ) represents the number of first vaccine doses administered to individuals of the population class {a,c} at time t. The value of !,# ( ) is inferred from data on the daily number of first doses administered by age group in Italy between the start of vaccination (December 27, 2020) and the end of the simulated period (June 30, 2021) [17]. Details are reported in Section 1.4.

Reproducing the SARS-CoV-2 epidemic trajectory in Italy
The reproduction number associated to the dynamical system considered can be computed as the dominant eigenvalue of the Next Generation Matrix (NGM) [21,22,23], defined as: Each block !,! 0 ) * ,;) + describes the contribution to the transmission of age-specific interactions between susceptible individuals in compartment 8  • ;) + is the relative infectiousness of SARS-CoV-2 cases among vaccinated compared to unvaccinated (equal to for vaccinated compartments, and 1 for the unvaccinated). ) the value of ( ). In particular, the selected value of ( ) will be the one that will make the model's reproduction number (recomputed after updating the ! 0 ) + ( ) to current state variables of the model) match the corresponding value of the net reproduction number as estimated on the same day from epidemic curves collected by the national integrated surveillance system [10,24] (reported in Supplementary Fig. 2).
Supplementary Figure 2. Estimates of the net reproduction number Rt the net reproduction number as obtained from epidemic curves collected by the National Integrated Surveillance System [24]. Line: mean estimates; shaded area: 95%CI.
Results discussed in the main text and in the following sections were obtained by running 300 simulations, sampling at each run a different value from the joint distribution of the transmission coefficient , the bootstrapped contact matrices !,! 0 and the relative susceptibility by age ! .

Model initialization
The model is used to simulate SARS-CoV-2 infection and vaccination dynamics in Italy between the start of the vaccination campaign (December 27, 2020) and June 30, 2021. The population by age was initialized according to the Italian age structure in 2020 [25] and statistics on underlying conditions, among which chronic respiratory disease, cardio-cerebrovascular disease, hypertension and diabetes [26] ( Supplementary Fig. 3).
Supplementary Figure 3. Fraction of the population having at least one underlying condition by age (%) [26] Fraction of initially immune individuals. The spread of SARS-CoV-2 in Italy throughout 2020 has been characterized by marked regional heterogeneities and estimates for the fraction of immune population in Italy by the end of December 2020 are not available. We considered three different Italian regions, that are representative of different levels of SARS-CoV-2 circulation: • Lombardy, the first and hardest hit region in 2020; For each region, we adapted a previously published model [1] to obtain estimates of the age-specific fraction of individuals infected with SARS-CoV-2 by December 27, 2020. Estimates obtained for Lazio were used to initialize the fraction of recovered population at the national level in the baseline analysis ("intermediate" immunity scenario), and those obtained for the two other regions were used in sensitivity analyses ("low" immunity: Campania; "high" immunity: Lombardy). For the baseline analysis (intermediate scenario), we obtained that about 16% of the overall population was infected with SARS-CoV-2 prior to the introduction of vaccination ( Supplementary Fig. 4), while for the low and high immunity scenarios these fractions are 9% and 23%, respectively ( Supplementary Fig. 5).
Under the assumption of natural immunity lasting on average 1/ 0 , we estimated the proportion f of recovered individuals that would have lost natural immunity by December 27, 2020. To this aim, we consider the time series of daily SARS-CoV-2 cases notified in Italy throughout 2020 [27] and, for each individual case i notified on day t, we determined if his/her natural immunity has waned by sampling the duration from an exponential distribution with average 1/ 7 . Left: Overall and age-specific fraction of population with natural immunity to SARS-CoV-2 in the low immunity scenario at the beginning of the vaccination campaign [1]. Right: The same as left, but for the high immunity scenario.
Initially infectious individuals. The number of infected individuals at the beginning of simulations (December 27, 2020) was determined in such a way to match the average daily number of notified cases reported in the week of simulations start (around 15,000) [27]. The reporting rate for SARS-CoV-2 infections at the end of 2020 was estimated as = / where is the infection fatality ratio estimated for Italy [28] and CFR is the case fatality ratio among SARS-CoV-2 cases reported in the last three weeks of December, assuming a delay between diagnosis and death of about 3 weeks [29]; thus, CFR= / where C is the cumulative number of cases reported in the last 3 weeks of December 2020 and M represents the number of deaths reported in the first 3 weeks of January 2021 [30]. The resulting value for the reporting rate is 41.2%.

Computation of age-specific vaccination rates over time
The rollout of the two-dose vaccination campaign is modeled using detailed data on the daily age-specific number of first doses administered over the considered period [17]. In the model, second doses are assumed to be administered to all individuals who have been vaccinated with one dose after an average of 42 days since the first dose. As shown in Supplementary Fig. 6, the model well reproduces the observed scale-up of the daily vaccination capacity occurred in Italy in the first half of the year 2021. The priority order of the Italian vaccination campaign was based on the WHO SAGE roadmap [31,32] prioritizing high-risk population age segments, i.e. over 80 years of age and essential workers (e.g. health care workers and teachers) and then progressively targeting younger age groups. In the absence of data regarding the presence or absence of underlying conditions in individuals targeted by vaccination, we assumed to distribute the daily doses  proportionally to the prevalence of underlying conditions in the age group considered. The evolution of vaccination coverage, by age and overall, as observed in the first half of the year 2021 is shown in Supplementary Fig. 6.

Vaccine efficacy against infection and death
In the first half of 2021 about 51 million vaccine doses have been administered in Italy, about 71% of which were Pfizer-BioNTech, 17% AstraZeneca, 10% Moderna and 2% Janssen [17]. Current evidence suggests an efficacy against infection with the Alpha variant of about 62% after two doses of AstraZeneca [33]. Estimates available for Pfizer-BioNTech suggest an efficacy against infection with Alpha of about 89% after two doses [34]. We inferred the corresponding vaccine efficacies after one dose of vaccine by assuming a 11% relative reduction compared to the efficacy observed after two doses [35,36] Fig. 7). We assume for the Moderna vaccine (mRNA type) the same vaccine efficacy of Pfizer-BioNTech and for Janssen (viral vector type) the vaccine same vaccine efficacy of AstraZeneca (Supplementary Table 1). Obtained age-specific estimates are in good agreement with estimates obtained for Italy [37] and range between 70.6% and 78.8% after the first dose (average 75.5%) and between 79.4% and 88.7% after 2 doses (average: 84.9%). Available estimates suggest a marked reduction in the risk of death in breakthrough infections: about 80.3% in vaccinated with one dose and about 96.4% in vaccinated with two doses [37]. These values include the reduced risk of SARS-CoV-2 infection in vaccinated compared to unvaccinated. To compute the infection fatality rate in breakthrough infections at different stages of vaccine protection < , we estimate an age-specific scaling coefficient ! ) + representing the risk of death given SARS-CoV-2 breakthrough infection in each stage < relative to infection in unvaccinated individuals: Where < death is set to 0 for unprotected individuals (j=0, 1, 5), to 80.3% for partially protected individuals (j=2, 3), and 96.4 for fully vaccinated individuals (j=4).

Model outputs
The main model outcomes are the age-specific number of new infections per day in the subpopulation i *,F ) + (t) for each vaccination stage j (including the unvaccinated, j=0); and the scaling factor δ(t) tuning the proportion of pre-pandemic contacts necessary reproduce the daily official estimates of Rt. Additional model outcomes are: • the number of SARS-CoV-2 confirmed cases between December 27, 2020, and June 30, 2021 and cumulative number of COVID-19 deaths over the same period; • the immunity profile of the Italian population on June 30, 2021; • estimates of the effective reproduction number (i.e., under the assumption of δ(t)=1).
For each model outcome, we report the mean values and 95% confidence intervals across stochastic simulations.

SARS-CoV-2 confirmed cases and COVID-19 deaths.
The daily number of SARS-CoV-2 confirmed cases is obtained from the daily infections by assuming a reporting ratio ( ) of 41.2% (see Section 1.3). The daily number of deaths in each subpopulation {a,c} with vaccination status < is obtained by applying the age-specific infection fatality ratio estimated for each population class ( !,# ) to daily infections, considering a delay of 3 weeks ( G ) between symptom onset and death [29], and rescaled by a factor (1 − ) to account for improvement in COVID-19 treatment. Specifically, where ! ) + is the relative risk of death in breakthrough infections, as computed in Section 1.5. To match the cumulative number of deaths detected by the surveillance system over the considered period, we obtain a reduction of mortality by 20%.
The age-specific IFR for each population class ( !,# ) was estimated as: Age group (years) ε F denotes the overall IFR in class c, as estimated from Lombardy data [28].
• ! denotes the age-specific IFR as estimated independently of the presence of underlying conditions [28]. • the scale factor ν is determined in such a way to minimize the root mean square error between ! and ̃! = ∑ !,# ⋅ !,# # , and !,# denotes the proportions of individuals of age a in class c in the Italian demographics [26].
Effective reproduction number. For each time t, the effective reproduction number is computed as the dominant eigenvalue of the Next Generation Matrix defined in Equation (1) by setting ( ) = 1, to account for complete resumption of pre-pandemic contacts (i.e., in absence of interventions and behavioral change).

No vaccination scenario
To evaluate the impact of the vaccination campaign, we need to define a suitable counterfactual scenario. Simply removing vaccination from the model and letting the epidemics evolve in absence of interventions would be an unrealistic scenario in terms of public health. Therefore, we chose as counterfactual a scenario where governmental restrictions and individual behavior would result in the same epidemic trajectory as the one observed in the presence of vaccination. We set the daily number of first doses !,# ( ) to 0 for all t, and recalibrated at each time step the value of ( ) necessary to match the time-series of the net reproduction number estimated from surveillance data ( Supplementary Fig. 2) [10,24]. Then, we compared the obtained estimates of social contacts, deaths, and effective reproduction number with those of the main analysis. ( ) = ( * ) where * corresponds to June 30, 2021. • adding a multiplying scale factor (1 + Delta ) when assuming that the Delta variant is dominant. The baseline value assumed for Delta is 50% [38][39][40]. This estimate neglects the increase in natural immunity due to SARS-CoV-2 infection occurring during the summer. This choice is supported by the low viral circulation in Italy in the summer, with 322,191 reported cases between July 1 and September 7, 2021, corresponding to about 1.3% of the total Italian population being infected after accounting for underreporting. By comparison, about 13.3% of all Italians have been fully vaccinated over the same period. We also did not mechanistically model the transmission dynamics of SARS-CoV-2 during the summer months because of the complexities arising from the co-circulation of two distinct strains. For the purpose of our estimates, the specific dynamics of COVID-19 transmission in the summer is irrelevant. In fact, the reproduction numbers estimated here only depend on the profile of population immunity at the considered date.

Future vaccination scenarios for the Delta variant
Future vaccination scenarios. We provide estimates on the reproduction number to be expected for the Delta variant under a set of vaccination scenarios , corresponding to different improvements of the vaccination coverage with respect to September 7, 2021. Vaccination scenarios are defined as follows: let Cov * ( ) be the vaccination coverage achieved by September 7, 2021 in age group , we define the age-specific vaccination coverage in age group under scenario as: where N;? is the minimum age to which the vaccine is administered (12 years in the baseline analysis, 5 years when considering a pediatric vaccine). We explored values of Ω between 60% and 100% with incremental steps of 5%.
For each vaccination scenario , we compute the reproduction number associated to different social contact levels * (between 0% and 100% of pre-pandemic contacts, with incremental steps of 5%) as the dominant eigenvalue of the Next Generation Matrix defined as: % ) * and ! 0 ) + (Ω) represents the number of individuals of age D with vaccination status < according to vaccination scenario .
As a baseline, we assume a transmissibility increase for the Delta variant with respect to Alpha ( Delta ) equal to 50% [38][39][40]. Alternative values of Delta (i.e. 25% and 75%) are explored as sensitivity analyses.

Description
We assess the robustness of our results with respect to alternative assumptions on: • the average duration of natural immunity (1/ 7 ); • the average duration of vaccine-induced immunity (1/ ) ); • the relative infectiousness of SARS-CoV-2 breakthrough infections ( ); • the natural immunity level at model initialization (end of December 2020) Supplementary Table 2 summarizes the main model parameters and assumptions in the baseline and in alternative sensitivity analyses.

Results
Supplementary Figs. 8-10 show that model results are most sensitive to the duration of natural immunity and the initial immunity profile, but that overall the estimates are quite robust. In particular, the estimate for June 30, 2021 of the average proportion of fully susceptible individuals ranges between 34% and 38% ( Supplementary Fig. 8), that of social contacts ranges between 44% and 53% of pre-pandemic contacts ( Supplementary Fig. 9) and that of the effective reproduction number ranges between 1.7 and 2.1 ( Supplementary Fig. 10).
Results obtained for future vaccination scenarios suggest that, if the vaccine did not reduce the infectiousness of breakthrough infections, the social contact levels that could be achieved without causing an epidemic would be about 45-65% (depending on the coverage scenario) compared to 55%-70% for the baseline ( Supplementary  Fig. 11). The duration of vaccine-induced protection affects very marginally results for the prospective vaccination scenarios ( Supplementary Fig. 12), due to the fact that the majority of vaccines have been administered too recently for waning to have a major effect. On the other hand, shorter duration of natural immunity, or a lower level of natural immunity at model initialization, would result in a slightly lower level of social contacts that could be resumed without causing an epidemic (namely, 50-70% for a duration of immunity of 1 year and 50-70% for an initial immunity profile that is lower than estimated -see Supplementary Figs Supplementary Figure 14. Sensitivity analysis with respect to the level of initial immunity. Heatmap of the mean estimated reproduction number for different vaccination scenarios (x axis) and different levels of social activity (y axis); n=300 stochastic model realizations, under different initial levels of natural immunity: low immunity scenario (left) and high immunity scenario (right).