The first 2 years of COVID-19 in Italy: Incidence, lethality, and health policies

Background The novel coronavirus disease 2019 (COVID-19) is an ongoing pandemic that was first recognized in China in December 2019. This paper aims to provide a detailed overview of the first 2 years of the pandemic in Italy. Design and methods Using the negative binomial distribution, the daily incidence of infections was estimated through the virus's lethality and the moving-averaged deaths. The lethality of the original strain (estimated through national sero-surveys) was adjusted daily for age of infections, hazard ratios of virus variants, and the cumulative distribution of vaccinated individuals. Results From February 24, 2020, to February 28, 2022, there were 20,833,018 (20,728,924–20,937,375) cases distributed over five waves. The overall lethality rate was 0.73%, but daily it ranged from 2.78% (in the first wave) to 0.15% (in the last wave). The first two waves had the highest number of daily deaths (about 710) and the last wave showed the highest peak of daily infections (220,487). Restriction measures of population mobility strongly slowed the viral spread. During the 2nd year of the pandemic, vaccines prevented 10,000,000 infections and 115,000 deaths. Conclusion Almost 40% of COVID-19 infections have gone undetected and they were mostly concentrated in the first year of the pandemic. From the second year, a massive test campaign made it possible to detect more asymptomatic cases, especially among the youngest. Mobility restriction measures were an effective suppression strategy while distance learning and smart working were effective mitigation strategies. Despite the variants of concern, vaccines strongly reduced the pandemic impact on the healthcare system avoiding strong restriction measures.


Introduction
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a new virus identified in Wuhan (Hubei, China) in late 2019 (1). SARS-CoV-2 causes the coronavirus disease 2019 (COVID- 19), an illness that ranges from mild flu symptoms to bilateral interstitial pneumonia (2). The virus spread so quickly around the world that the World Health Organization (WHO) declared the COVID-19 a Public Health Emergency of International Concern on January 30, 2020 and a pandemic on March 11, 2020, (3). Unlike other coronaviruses, the SARS-CoV-2 is able to spread through pre-and asymptomatic infections that are difficult to detect and isolate, requiring health authorities to test all contacts of confirmed cases to lower the risk of spread (4,5). The lethality of the original strain was estimated using infection fatality ratios (IFR) assessed through several national sero-surveys (6,7). While relatively low in the whole population (<1 death per 100 infections in developed countries), the risk of death is shown to increase with age (up to 10-15 deaths per 100 infections in people aged more than 75 years) and in patients who are immunosuppressed or have concomitant comorbidities (8,9). Furthermore, since the prognosis of severe cases depends on the availability of intensive care beds, lethality increases when critical care capacity is saturated (10). To address the pandemic, a global vaccination campaign was launched, and pharmaceutical industry developed candidate COVID-19 vaccines at an unprecedented speed. By the end of 2020, global Medicines Agencies had conditionally approved several vaccines based on different technologies, with others close behind (11,12). During the first 2 years of the pandemic (since December 29, 2019, to February 28, 2022), National Health Institutions detected 444,900,763 confirmed cases and reported 6,020,752 deaths worldwide . The highest number of infections favored mutations in the viral genome sequence and led to generation and spread of many viral variants (13). WHO coordinates national and subnational research aimed at sequencing RNA viral genomes detected in infected people and classifying variants of concern (VOC) that may pose a greater risk to global public health . From May 2020 to February 2022, WHO identified five consecutive VOCs: Alpha, Beta, Gamma, Delta, and Omicron. Each variant showed an increased capacity to spread (even within vaccinated people) and although the debate on virulence is still open, it would appear that all the VOCs except Omicron caused a disease with higher severity and mortality (14,15). In February 2020, Italy was the European epicenter of the SARS-CoV-2 spreading. The unexpectedly high speed of transmissions quickly resulted in hospital saturation and forced the Italian government to establish a national lockdown. Restriction measures blocked the first wave and were gradually removed in parallel with the development of a robust COVID-19 contact tracing system. To avoid lockdown during the second wave, the national government has applied a standard set of restriction measures (from soft to hard) at the regional level based on the risk of spread evaluated on a weekly basis. The risk level by geographic area (represented by a colored map: white = low, yellow = moderate, orange = high, red = highest) was evaluated by determining weekly estimates of incidence and reproductive number (R t ). During the second wave (December https://covid .who.int (accessed September , ).
27th, 2020), a national vaccination campaign was launched using two messenger RNA (Pfizer-BioNTech, Moderna) and two vector vaccines (Janssen, Vaxzevria) (14). Given the high percentage of vaccinated people in the third wave, the hospital saturation levels replaced the R t in the risk evaluation. Although vaccine protection declined over time (especially against virus variants), protection returned following administration of the booster dose especially against the development of severe infections (16)(17)(18)(19). Health institution recommended a booster shot after 4 months from the standard cycle in September 2021, and included children aged 5-11 years in the vaccine campaign in December 2021 (19,20). This study aims to provide a detailed overview of the first 2 years of the pandemic in Italy, where 13,000,000 of confirmed cases and 155,000 deaths were reported from February 2020 to 2022. The current paper is part of a larger project aimed at describing the epidemiology of Italian COVID-19 pandemic and follows an initial article introducing the method used to describe the pandemic in its 1st year (21).

Study design
This study analyzed public data of COVID-19 in Italy collected in the national registry by the Civil Protection (CP) and the National Health Institute (ISS).

Settings
The Italian Government declared a health emergency status on January 3, 2020 and extended it to March 31, 2022. The CP was delegated to manage the process and established a system to collect COVID-19 data in a national registry (managed by the ISS). Aggregate data on incidence and vaccination are published daily. The ISS reviews and updates the registry data to account for data reporting delays and regional recounts and releases an updated report with details including the age distribution of detected cases. The European Center for Disease Prevention and Control (ECDC) collects VOCs continental data through the European Surveillance System (TESSy).

Participants
All confirmed cases of COVID-19 in Italy.

Outcomes
The primary outcomes were: (1) the number N k of persons who became infected on the kth day of pandemic; (2) the number D k of persons who died (over time) among N k (i.e., the

Data sources/measurement
Aggregate data from the national COVID-19 registry and the vaccine campaign are stored in public repositories and updated daily. The data include daily counts of performed tests, of diagnosed cases and fatalities who tested positive using the polymerase chain reaction or the rapid antigen test (beginning on January 8, 2021), and of persons who received vaccine shots by region . The ISS provides a weekly report that includes the median age of detected cases, estimates of vaccine protection and (beginning on December 7, 2020) the distribution of detected cases by 10-year age class (22). The ECDC releases European data on VOCs , the National Institute of Statistics releases data from the sero-survey (May 25-July 15, 2020) and on Italian population .

Statistical analysis
As already highlighted by De Natale et al. at the onset of the pandemic, the high number of asymptomatic infections makes deaths more suitable than detected cases for estimating incidence (23). Given the probability p k of dying after having caught the infection on the kth pandemic day (k ∈ Z + ), we used the negative binomial distribution to estimate the daily number of infections (N k ) from the resulting deaths over time [D k ; Section Modeling the Incidence of Infections (Negative Binomial Distribution)]. First, we estimated D k by applying the weighted moving average to deaths (that are recorded by the occurrence date, Section Estimating D k and υ k : Weighted Moving Average). Second, we modeled the probability p k accounting for the age at infection, VOCs prevalence and population vaccination level (Section Modeling the Daily Probability to Die p k ). Using other simple assumptions, we evaluated excess death (for health system saturation) and lives saved by vaccines (Sections Excess Death and Vaccine Effect). Finally, we used the number of detected cases υ k among N k to check the admissibility (N k > υ k ) of estimates (Sections Estimating D k and υ k : Weighted Moving Average and Checking Estimates). In the following, we will proceed with the mathematical formulation, which will be progressively upgraded, in the next sections, to consider the more complex probabilities involving age classes, different strains, and vaccination level. Once the main formulas are established, the estimated variables used to determine the solutions will be given in the Section Estimating the Daily Lethality.

Modeling the incidence of infections (negative binomial distribution)
Let X (j) k be the binary random variable representing the outcome (1 = dead; 0 = recovered) of the jth person infected on the kth day of the pandemic be the random variable representing the rank of the daily infection resulting in the D k -th death, the probability of N (D k ) k follows a negative binomial distribution with parameters D k and p k P N with k, D k ∈ Z + , n ≥ D k . We estimated the number of daily infections (with the related 95% CI) as the mean of Equation (1) Estimating D k and υ k : Weighted moving average Let d k,k+j and υ k,k+j be the number of persons infected on the kth pandemic day who died or were diagnosed j days after the infection, the number of deaths (D k ), and detected cases (V k ) among infections on the kth pandemic day can be evaluated as Since only the corresponding number of events by the occurrence date (of death or diagnosis) is available D k and V k were estimated as Let T dead and T diagn represent the time from infection to death and diagnosis, respectively, and α k and β k be the binary variables representing the events to die (α k = 1) or be alive (α k = 0) and to be diagnosed (β k = 1) or undetected (β k = 0) on the kth pandemic day, π (k+j) j and θ (k+j) j can be expressed as the conditional probability to die or be diagnosed j days after the infection The ISS provided estimated quartiles (Q1, Q2, and Q3) of the time distributions from symptoms to death and diagnosis during three different periods (March-May/2020, June-September/2020, and October/2020-December/2020). The ISS estimates for time to death are admissible under symmetric distributions except during the summer period [where there was strong bias from clusters of vacationers (24)]. These biased estimates were not considered and the remaining, which are equivalent [ Table 1 in (21)], were extended to the whole studied period. We added 5 days [the mean time from infection to symptoms (25)] to ISS estimates to obtain the corresponding parameters of the probability density function of the time from infection to death and diagnosis If necessary, we adjusted for symmetry by replacing the median with the center of first and third quartile and assumed that the functions in Equation (6) follow the truncated normal distribution where µ and σ are the mean and standard deviation of the parent general normal probability with µ = Q 3 +Q 1 2 and σ = Q 3 −Q 1 1.34896 . Of note, the Equation (4) with probabilities Equation (5) derived from Equation (7) can be also interpreted as a weighted moving average of period 2µ + 1 on time series d .,k+j and υ .,k+j in (3) Modeling the daily probability to die p k Let X j,ξ ,V and Y k,j,ξ ,V be the binary random variables representing the events "to die after the infection" and "to be infected on the kth pandemic day, " respectively, by 10-year age class (j: 0-9, 10-19, . . . , 80-89, 90+ years), VOC (ξ : 0 = original strain; 1 = Alpha; 2 = Beta; 3 = Gamma;4 = Delta; 5 = Omicron), and vaccination level (V: 0 = unvaccinated; 1 = uncompleted basic cycle; 2 = completed basic cycle more than 4 months ago; 3 = completed basic cycle in the last 4 months; 4 = received a booster shot). By assuming that the conditional probability p k,j,ξ ,V to die after having caught the infection on the kth day does not depend on k, we have that Let N k,j,ξ ,V be the number of infected people on the kth pandemic day by age class, VOC, and vaccination level and N k,.,.,. be the total number of infections on the same day, the overall probability p k to die among infections on the k-th day is equal to Now, let RR j,ξ ,V be the risk ratio to die of people with the vaccination level V (= 0, 1, 2, 3, 4) vs. unvaccinated (V = 0) by age class and VOC and let N k,j,ξ ,. and N k,j,.,. be the number of infections on the kth pandemic, respectively, by age class and VOC (with any vaccination level) and by age class (with any VOC and vaccination level), the Equation (9) can be rewritten through the Equation (10) as Finally, let S j,ξ ,0 (t) and S j,0,0 (t) be the distributions of survival time of unvaccinated people in the jth age class, respectively, for the VOC ξ (= 0, 1, 2, 3, 4, 5) and original virus strain (ξ = 0), under the assumption of proportional hazards we have that where hR j,ξ ,0 are the hazard ratios by VOC (ξ = 0, 1, 2, 3, 4, 5 vs. ξ = 0) by age class for unvaccinated people (V = 0). By integrating the Equation (12) over the whole pandemic period, we obtain the following identity Frontiers in Public Health frontiersin.org Ferrante . /fpubh. . and since S j,ξ ,0 = 1 − p j,ξ ,0 and S j,0,0 = 1 − p j,0,0 the Equation (11) can be expressed as The Equation (14) is the lethality equation I introduced to compute the pandemic parameters of interest. We can notice that it depends on k only through the daily distribution of infection by age, VOC, and vaccination level and that we can derive the lethality by variant as is the age distribution of infections due to the variant ξ . All still unknown quantities used to univocally determine the result will be specified in the paragraph Estimating the daily lethality.

Excess death
Let P and P j be, respectively, the Italian population and its subgroup in the j-th age class and Y (j) i be the binary random variable indicating that the ith person in the jth age class has been infected. If the virus spreads randomly within the population, we would have that the distribution of cases by age class equals that of the whole population and the related probability of dying can be obtained from Equation (14) by replacing the proportion of infected people by age class ( ) with the corresponding proportion in the population ( Since COVID-19 transmission began among younger people and eventually spread within the elderly (26), it was assumed that the spread was out-of-control if the distribution of detected cases by age class followed the age structure of the population (15). Through the deaths that resulted from the product between the Equation (2) and the Equation (16), the excess death was defined as the following difference:

Vaccine e ect
Avoided infections By rewriting the conditional probability in Equation (8) as ratio of probabilities, the relative risk in Equation (10) can be expressed as Under the assumption that the vaccines had no impact on the risk of catching the infection (P Y k,j,ξ ,V = 1 = P Y k,j,ξ ,0 = 1 ), the relative risk (17) reduces to: Finally, let Pop k,j,ξ ,V be the population at risk on the kth pandemic by age class, VOC, and vaccination level and D k,j,ξ ,V be the number of deaths in each group, through the relationship D k,j,ξ ,V = P X j,ξ ,V = 1, Y k,j,ξ ,V = 1 Pop k,j,ξ ,V , the relative risk (18) can be rewritten as By replacing in Equation (14)  (20)

Saved lives
If the vaccines have no effect against death, the relative risks RR j,ξ ,V in Equation (14) would be equal to 1 and the lethality would reduce to  (21), we obtain an estimate of the number of deaths D * * k that would have occurred without vaccines By multiplying the Equation (2) for the Equation (21), we obtain an estimate of the number of deaths D * k that would have occurred without vaccines among the infected people

Checking estimates
We studied the ratios r k of detected cases υ k among estimated infections N k on kth day (Figure 2) If r k (i) >1 (i.e., υ k > N k ), then the estimatedp k overestimates the actual p k on the kth pandemic day p k <p k .

Estimating the daily lethality
Available data to estimate quantities in Equation (14) were used as follows: 1) The probability p j,0,0 of dying among unvaccinated (V = 0) people in the jth age class who were infected with the original strain (ξ = 0) was estimated using the IFR by age class ( IFR j,0,0 ) in (6) p j,0,0 = IFR j,0,0 .
2) The daily unvaccinated population by age (N k,j,.,0 ) was estimated as the difference between the ISTAT population 7 and the vaccinated people 4 . 3) As estimates of hazard ratios hR j,ξ ,0 were used those from (27,28) and since those for Alpha, Beta, Gamma, and Delta are not determined by age class, we considered them constant by age. 4) Let D j,ξ ,V , C j,ξ ,V , and Pop j,ξ ,V be the number of deaths, of detected infections, and of population by age class, VOC, and vaccination level, respectively. The ISS provided estimates of the relative rate (vaccinated/unvaccinated) of deaths (RD j,ξ ,V ) and of infections (RC j,ξ ,V )  (22). We used those estimates to assess the relative risk in Equation (10) and in Equation (19) as follows 5) As estimates of the fraction of vaccinated infections by age class and VOC ( ), the daily fraction of vaccinated population by age class 4 were used.
Estimating f Let P j and P (Med k ) k,j be the people of age j in the Italian and in fictitious populations, respectively, on the kth pandemic day, using the definition of "median" we have that where 100+ indicate people aged 100 years and more. By assuming that the ratios between infected people at ages greater than the median and the oldest (100+ years) are equal to those in the Italian population, and that the ratios between those at ages smaller than or equal to the median Frontiers in Public Health frontiersin.org Ferrante . /fpubh. . and the youngest (0 years) are also equal to those in the Italian population and P (Med k ) k,100+ are determined from the Equation (24) and can be used to derive all the remaining fractions (f Since the method returns one probability estimate per week (ISS median age refers to a week), each pair of values was linearly connected. By assuming RR j,ξ ,V = hR j,ξ ,0 = 1 in Equation (14), the resulting death probability is equal to By replacing the Equation (26) in the Equation (2), we estimated the number of infections ( N) from the beginning to the middle day of the ISTAT seroruvey (June 19, 2020; after 120 gg from the beginning) as

Waves
In epidemiology, an internationally accepted definition of "wave" does not still exist; the term refers to the appearance of a plot of cases over time. In this paper, the word "wave" is used to indicate the part of the plot that lies between two local minima.

Health policy evaluation
The effects of applied health policies were determined by comparing the weekly and bi-weekly incidence rates before and after the day (k) they entered into force with k ∈ Z + and j = 7, 14.

Data cleaning methods
Dates of the ISS cumulative distribution of detected cases by age correlate with the most recent update. Those distributions were reordered to have non-decreasing functions. Since the VOC prevalence data in TESSEy cover a week, the data were linearly fitted to obtain daily prevalence. In addition, frequencies recorded before the official date of the first available samples were set to zero.

Incidence curve
The incidence curve in Figure 1 shows five waves. The   (Table 1). From the third to the fifth wave, the proportion of infections among young individuals (0-19 years) increased by up to 30% (Supplementary Figure 1).

Lethality
Estimates of infection-related deaths during the 1st months of the pandemic were provided by a fictitious population (26) and adjusted using a correction factor (27) of 1.02%. Daily lethality ranged from 2.8% (first wave: April 9, 2020) to 0.15% (last wave: December 30, 2021), causing a total of 152,358 deaths with a peak of 723 deaths on November 5, 2020 (Figures 2, 3). After the peak, lethality in the first wave decreased to 1.0% on June 11, 2020, generating 32,739 deaths. Initially, lethality of the second wave, which caused 62,595 deaths, decreased to 0.8% on August 15, 2020, and then increased to a peak of 1.6% on December 8, 2020, then remained stable. The lethality in the third wave peaked at 1.4% on February 8, 2021, then decreased to 0.45% by June 8, 2021, and caused 28,596 deaths. The lethality in the fourth wave, which caused 3,620 deaths, initially decreased to 0.38% by July 6, 2021, and then increased to a peak of 0.91% by September 30, 2021. The lethality of the fifth wave, which caused 24,808 deaths, decreased continuously from a peak of 0.91% on October 2, 2021 to 0.15% on December 30, 2021, and then slowly increased to 0.18% by the end of February 2022. Two periods had the lethality higher than the threshold: the first 4 months of the first wave (March-June, 2020), with an excess of 14,134 deaths; and the last 3 months of the second wave (November, 2020-January, 2021), with an excess of 6,478 deaths ( Figure 3).

Impact of variants of concerns
Of five VOCs, three became prevalent (Alpha, Delta, and Omicron), each replacing prior variants at a faster  pace. Omicron was responsible for the higher number of infections (7,557,368), while Alpha was associated with the largest number of deaths (29,167). Without vaccines, Delta would be the most virulent VOC with > 70,000 deaths (Supplementary Table 1).

Vaccine e ect
Vaccines reduced infections by 38% (from 25,045,987 to 15,604,551) and deaths by 62% (from 185,850 to 71,760). Of 114,090 lives saved, 62,902 (55%) would have resulted from infections prevented and 51,188 (45%) from the infections that occurred. Without vaccines, the expected number of infections would have been 30,274,455 (30,160,115-30,389,006), and the expected number of deaths 266,448 for a lethality of 0.88% (Figures 1-3 and Table 1).

Health policies e ects on estimated curves
The strongest restriction measures affecting all the population (initial stay at home, the November 2020 introduction of standardized prevention measures based geographic risk, and the Christmas 2020 and Easter 2021 restrictions) strongly reduced the curve rates (up to 1,000%) within a 1st week of their introduction. Industrial lockdown and specific restrictions (including 75% of high-school students who received distance learning) implemented in October 2020 are associated with smaller (from −37 to −227%) and slower (concentrated in the 2nd week) rate reductions. Curve rates increased after school openings (except after the last one of January 2022) and reductions in smart working, especially during the 2nd week. Rate increments after school openings reduced over time. In-shop Christmas 2020 incentives increased the incidence rate by 90-100% and the death rate by 45-60%. The introduction of a compulsory green pass reduced the rate of infection and death curve by 150-200 and 40-70%, respectively. The mandatory use of the FFp2 mask in closed places reduced the curve rates by 45-85% ( Table 2). The introduction of rapid tests (from January 2021) increased the percentage of infections detected among children, particularly when schools were open (Supplementary Figure 1).

Discussion
This paper provides a comprehensive picture for the first 2 years (February 2020-2022) of the COVID-19 pandemic in Italy, including the impact of VOCs, the vaccine campaign (until the third shot), and an evaluation of government health policies using only public data.

Virus spread
Almost 40% of COVID-19 infections have gone undetected, likely because they were asymptomatic or paucisymptomatic. During the first wave, the virus primarily spread in the north of the country and was highly concentrated in the Lombardy region. The virus likely arrived in Italy through the airport system of Milan (the largest city in Lombardy), which includes one intercontinental and two international airports (one of which in Bergamo, the most hard-hit Italian city). Like other respiratory viruses, SARS-CoV-2 spreads directly or indirectly through person-to-person contact (especially in indoor environments) (29). Lombardy is the Italian region with the highest level of daily commuting for work or school . Indeed, a study highlighted the association between the regional patterns of viral spread during the first wave and the origindestination matrix of goods and food transportation and for the https://www.istat.it/it/archivio/ population (30). Another component responsible for the rapid spread of COVID-19 was unpreparedness. During the onset phase of the pandemic, hospitals followed WHO guidelines and tested only people with a known link to China (thus accelerating the spread of the virus). Fortunately, the quick stop to all national mobility on March 12, 2020, confined the virus to northern regions. Data collected during the second wave revealed similar patterns to those reported by other studies: the virus infects younger people first followed by those >70 years of age (26). Since retired people have fewer daily contacts than students and workers, who often use public transport and share indoor environments with others, it is likely that school and work transmission impacted the onset of the familial transmission chain among the elderly. Infected students and workers carried the virus home, transmitted the infection to other family members (31) and increased the probability that older relatives (including grandparents) would become infected, especially through presymptomatic or asymptomatic infections. During the summer months, those of 20-29 years of age were the  hardest hit, presumably because of increased nightlife and other social activities. Without public health policies, it is likely that two waves would have occurred per year (similar to the fourth and fifth waves as shown in Figure 1): a winter wave (resulting from a higher number of transmissions from indirect contacts) and a summer wave (largely resulting from direct contacts). The former is longer (October to June), peaks in January, and is more virulent because it mainly involves families; the latter is shorter (July to September), peaks in August, and is less virulent because it primarily involves single people.

Lethality
The estimates of lethality obtained by the fictitious populations (26) provided a cases count by June 20, 2020, that was very close to that estimated by the national ISTAT sero-survive 6 , with an error of 2%. This indicates that with a high median age of detected cases (with respect to that of the national population), the fictitious population provides reasonable estimates when the age distribution of cases is unknown. The virus's lethality was extremely high in the first 3 months of the pandemic, when the median age of detected cases was much higher than that of the population. This is the result of two serious errors: a lack of screening tests to reduce transmission from adult/young to elderly and the use of nursing homes to support hospitals with saturated capacity in the hardest hit regions (which increased infections among the most at risk population) (32). Lethality was higher than expected even during the second wave. The introduction of reliable rapid tests allowed a massive test campaign that kept lethality under the threshold from the right tail of the second wave. Lethality would have been under the threshold even without vaccines (Figure 2 dotted curves). Of the three prevalent VOCs in Italy, the first two were more virulent than the original strain, the third was not. To ensure that vaccines remain updated and new and potentially harmful variants are identified early, ongoing research on the evolution of virus genomes is crucial.

Health policies
The stronger the restriction measure, the higher its efficacy and the more quickly it took effect. Supporting the assumption that the incidence curve was largely underestimated at the onset of pandemic, the initial lockdown impacted the death rate more than the infections rate. After the second wave, the large-scale screening helped to monitor the actual size  of the outbreaks, especially among young students (often asymptomatic). Quarantining infected grandsons (tested at school) and parents (tested at work) likely protected the grandparents. This is supported by increased rates of infection and death curves after school openings (except the last one) and reductions in smart working. However, in schools, those increments declined over time until disappearing, presumably because school protocols became more and more effective. Even if the mandatory use of the ffp2 in the public transport reduced the curves rate by up to 70% during 2021 Christmas holidays, their true effect is shown after the schools opening, where rates drastically decreased of 1,500-2,800%. Although protective effect of the vaccines was reduced by the emergence of new variants, vaccination saved more than 110,000 lives and avoided the saturation of the health system without a need for stronger restriction measures even during periods of high virus circulation.

Advice
It is necessary to monitor the evolution of SARS-CoV-2 in greater depth and to develop mathematical models that can predict future changes in its genome. A flexible pandemic plan able to adapt to the evidence of the data (collected through a digital and multi-connected surveillance system) should be developed through a multidisciplinary approach and shared with international health authorities. It should contain measures that are tailored for different combinations of virus transmissibility and virulence (from low-low to high-high). Estimates of territorial origin-destination matrices can help to simulate possible spatiotemporal patterns of virus spread. Initial settings of a public health response should refer to an "average" or "worst case" scenario and updates should follow data evidence.

Data availability statement
Publicly available datasets were analyzed in this study. This data can be found at: https://github.com/floatingpurr/covid-19_ sorveglianza_integrata_italia/tree/main/data.

Ethics statement
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent from the participants' legal guardian/next of kin was not required to participate in this study in accordance with the national legislation and the institutional requirements.

Author contributions
PF: study concept and design, acquisition of data, analysis and interpretation of data, drafting of the manuscript, critical revision of the manuscript for important intellectual content, and statistical analysis.