key: cord-0955399-6zuj5ffk
authors: Fierro, A.; Romano, S.; Liccardo, A.
title: Vaccination and Variants: retrospective model for the evolution of Covid-19 in Italy
date: 2022-02-28
journal: nan
DOI: 10.1101/2022.02.27.22271593
sha: 5cd7a31b5777e0b6d76e893065e08519c451e533
doc_id: 955399
cord_uid: 6zuj5ffk

The last year of Covid-19 pandemic has been characterized by the continuous chase between the vaccination campaign and the appearance of new variants that put further obstacles to the possibility of eradicating the virus and returning to normality in a short period. In the present paper we consider a deterministic compartmental model to discuss the evolution of the Covid-19 in Italy as a combined effect of vaccination campaign, new variant spreading, waning immunity and mobility restrictions. We analyze the role that different mechanisms, such as behavioral changes due to variable risk perception, variation of the population mobility, seasonal variability of the virus infectivity, and spreading of new variants have had in shaping the epidemiological curve. The fundamental impact of vaccines in drastically reducing the total increase in infections and deaths is also estimated. This work further underlines the crucial importance of vaccination and adoption of adequate individual protective measures in containing the pandemic.

The Covid-19 pandemic dramatically impacted on all aspects of world population life, causing a global lasting damage at economic, social and educational level and an enormous loss in terms of human lives. The extraordinary effort made worldwide for the development of vaccines against Covid-19 allowed an equally extraordinary achievement, such as the approval of the first vaccines in one year since the beginning of the pandemic. Such a circumstance made concrete the possibility of finally defeating the virus, limiting the further use of emergency measures, such as lock-downs, no longer economically sustainable. The first vaccine authorized by the US Food and Drug Administration for distribution in the United States and by the European Medicine Agency (EMA) for the European Union (EU) countries were the mRNA-vaccine Comirnaty (BNT162b2), produced by BioNTech/Pfizer (December 2020) and soon after the one produced by Moderna (mRNA-1273). At the same time, the adenovirus viral vector vaccine Russian Sputnik V, produced by the Gamaleya Research Institute of Epidemiology and Microbiology, was authorized and distributed in different countries. One month later another adenovirus viral vector vaccine, the Vaxzevria (ChAdOx1-S), developed by the Oxford University and produced by Astrazeneca was authorized by EMA for distribution in EU.

vaccine effectiveness against hospitalization essentially in agreement with British data, but with a 30% absolute reduction in effectiveness against disease [26] . More optimistic values for the effectiveness of vaccination against symptomatic disease caused by the Delta variants are reported in [27] , where it is evaluated at 88.0% for mRNA vaccine and to 67.0% for viral vector vaccine. Both for the Comirnaty [28] and Moderna vaccine [29, 30] , a booster dose with the same vaccine with the original SARS-CoV-2 spike protein has been strongly recommended to enforce the protection also against the Delta variant. As stressed in [31] it is hard to differentiate the effectiveness reduction against new variant from the natural decay of immunity as time passes. According to the retrospective cohort study conducted in USA in the previously cited paper, the reduction in vaccine effectiveness against Covid-19 infections over time is probably primarily due to waning immunity with time rather than the Delta variant escaping vaccine protection. In [33] the waning effectiveness of the vaccine in England, at 20 weeks or more after vaccination, was estimated to be 44.3% with Vaxzevria and 66.3% with Comirnaty, while the vaccine efficacy against hospitalization and death was confirmed. Similar results were found in data collected from the Israeli national database [34] , and in Qatar [35] . It has been also ascertained that the immunity response decreases faster in elderly people than in young individuals [36] . This paper is devoted to the retrospective analysis of the evolution of the Covid-19 in Italy, keeping into account the vaccination rate, the variant spreading and the immunity decay. To this purpose we generalize the model developed in [37] , with the introduction of the appropriate compartments for vaccinated individuals and suitably modify the parameters in order to simulate the increase in prevalence of the Delta variant, starting from mid-May and becoming dominant by August. Other variants, different from Alpha and Delta are not considered in the present model. The effect of waning immunity as time passes since the completion of the vaccine cycle and the occurrence of breakthrough infections is also included in the model. We weigh the role that different mechanisms, such as variation in the population mobility, seasonal variations of the virus infectivity, variation in risk perception and appearance of new variants have had in forging the evolution of the epidemiological curve. We find that the most relevant mechanism is the seasonal variation in the stability of the virus, followed by the awareness mechanism, that induce individuals to increase/relax self-protective measures when the number of active cases increases/decreases. The appearance of the Delta variant and the mobility variations have had instead only marginal effects. We also estimate the effect of the vaccination campaign: in absence of vaccines the emerging scenario would have been dramatic with a percentage difference in the number of total infections and total deaths respectively equal to +63% and +55%. The model also predicts the appearance of a new variant (the Omicron variant) and its becoming dominant in January 2022. the vaccinated compartments V j (with j = 1, 2, 3 indicating the first dose, the full vaccination, i.e. the second dose except for the single dose Janseen vaccine, and the third dose, respectively), the asymptomatic infected compartments I i A (t), the symptomatic infective compartments I i S (t), the diagnosed compartments I i D (t), the dead compartments D i (t) and finally the recovered compartment R(t), that includes both unvaccinated and previously vaccinated individuals. We assume that healed individuals can loose immunity over time and return to the susceptible compartment, after an appropriate time. We do not consider demographic birth and (not related to the virus) death process. The epidemic dynamic is thus governed by the fluxes of individuals among these compartments, as shown in Fig.1 , and it is fully described by a system of fifteen coupled first order differential equations for the normalized SEV I A I S I D RDS variables:

The system is closed and positive, i.e. all the state variables take non negative values for t ≥ 0, if initialized at time 0 with non negative values, and satisfy the mass conservation lawṠ + jV j + i Ė i +İ i A +İ i S +İ i D +Ḋ i +Ṙ = 0, hence the sum of the states (the total population) is constant. In our model individuals move to the vaccinated compartments V j only when the vaccine protection becomes effective, two weeks after the inoculation. Let us briefly review the main characteristic of the parameters. During the first wave of the pandemic, the world had to face with a completely novel virus and, as a consequence, it took some time to understand its mechanism of action, effective therapies and treatment modalities. As a consequence, most of the parameters that governed the evolution of the epidemic in the first period significantly changed in time. After more than one year, the knowledge and the experience in preventing, diagnosing and treating the infection made some of these parameters to reach an almost stable value. However some epidemiological parameters, as for instance individual mobility, seasonal stability of the virus, risk perception, immunity etc., are intrinsically time dependent. In order to fix them, we follow the same approach as in [37] , i. e. we try to avoid step-functions and look for reasonable functional behaviors to describe their evolution. In the following we describe the principal time dependent parameters, while the constant ones are reported in Tables 3 and 6 . The transmission rates, α and γ The model has two different transmission rate parameters, α and γ, which govern the transmission of the virus respectively from undiagnosed and diagnosed individuals. The February 23, 2022 4/20 transmission from undiagnosed symptomatic or asymptomatic individuals is still the dominant mechanism in the spread of the epidemic. Following [38] [39] [40] , we assume the transmission rates from undiagnosed vaccinated individuals, α v , to be different (and smaller) with respect to the corresponding rates for unvaccinated people, α u . Widely proven isolation protocols allow to assume the transmission of the virus by diagnosed cases, γ, to be residual and equally effective both from infected vaccinated and unvaccinated individuals, thus we do not differentiate the parameter γ between u and v individuals.

The parameters α i can be both factorized in a pure contact term, α c , that we assume to be independent on the vaccination status 1 , describing the probability per unit time that a susceptible individual meets an infected individual, and the susceptibility term, σ i (t), which takes into account the probability that a potentially contagious contact between a susceptible and an infected individual leads to a new infection:

Both terms are in principle subject to changes. Restrictive measures, such as lock-down and limitations of access to places and services, certainly affect the contact rate, α c , as it happened during the first period of the pandemic. In [37] the global effect of such significant modifications of the mobility was encoded in an appropriate mobility function obtained as the weighted average on the mobility data from the Google 

where the function is modulated in such a way to produce a doubling of the average number of contacts among individuals between February/March 2021 and the autumn period (i.e. c 1 + c 2 = 2 and c 1 − c 2 = 1), and t c and τ c are fitting parameters (the values of parameters are shown in Tables 4 and 5 ). The susceptibility terms σ i (t) depend on different factors, partly related to the behaviors of individuals -e.g. more careful use of self-protective measures such as face masks, hand washing, due to increased risk perception during the rising phases of the epidemic -and partly related to the characteristics of the virus such as seasonal variation in the stability of the virus in airborne [42] [43] [44] and increased transmissibility due to the appearance of new Covid-19 variants. The functional form of the susceptibility function is thus chosen as

where σ aw is the term encoding the awareness mechanism due to risk perception, σ term encodes the modification of transmissibility due to seasonal variability, σ var encodes the emergency of new variants and the last factor, σ i 0 , differentiate between vaccinated and unvaccinated individuals. Fixing σ u 0 as a fitting parameter, we assume the transmissibility from vaccinated cases to be lower (σ v 0 = 0.5 · σ u 0 ), according to studies on household transmission [38, 45] . Recent literature seems to question the viral load of vaccinated infected individuals to be lower, however its decline seems to be faster than for unvaccinated individuals [46] .

Following [37] we assume a prevalence based mechanism of rising awareness that increases the risk perception, inducing individuals to adopt more protective behaviors, with the effect of reducing the transmissibility of the virus during the peak. Thus we assume the awareness mechanism to act on the transmission rate by reducing it with a factor inversely proportional to the number of infective detected individuals, without any temporary effect of amplification or falsification [47] :

.

The reduction of contagiousness during the warm season, due to a potential decline of the stability of the virus in warm environment [42] [43] [44] [45] , already considered in [37] , seems to be further confirmed by more recent literature [48] [49] [50] . Thus we mimic the susceptibility decrease/increase, in spring and autumn respectively, through appropriate hyperbolic tangent functions:

where t term1 , t term2 , τ term1 and τ term2 are fitting parameters 2 . Finally the term

is the one encoding the susceptibility increase due to the circulation of new and more infective variants. In particular, following the ISS data [51] concerning the prevalence of Delta variant in Italy, the exponential increase is regulated in such a way to mimic the exponential transition from the susceptibility of the Alpha variant to the enhanced susceptibility of the Delta variant, that became dominant at the end of July (91.4% on July, 26). Following [21] and [23] , ∆ σ is fixed a priori equal to 0.60. Due to the circulation of the Delta variant also the parameter γ increases:

where γ 0 is a fitting parameter. The susceptibility reduction function for vaccinated people, f j Vaccinated people receive partial immunity against the infection, with an effectiveness that increases with the number of doses and decreases with the time occurred since the last inoculation and eventually with the appearing of new variants. In particular, the first dose gives only negligible immunity (estimated around 30% for both BNT162b2 and ChAdOx1 nCoV-19 vaccines in England [27] ), whereas the third dose gives optimal immunization against the Delta variant [52] (estimated around 90% for BNT162b2 mRNA Vaccine in Israel). For what concerns the second dose, as previously discussed, different timing leads to extremely different levels of protection among the double dose vaccinated individuals. It has been nowadays established [53] that the effectiveness of vaccines against the infection decreases over a period of 6 months from the inoculation of the second dose (with an estimated decay of 40%), whereas the efficacy against severe manifestation is subject to a minor decay (less than 15%). This makes it difficult to quantify the average vaccine protection level to be included in an average field model such as the one considered in the present work. To this purpose we introduce a susceptibility reduction function for vaccinated people in Eq. (1) defined as f j = 1 − E j , where E j (with j = 1, 2, 3) is the vaccine efficacy after one dose, full vaccination and booster dose, respectively. Following literature, efficacy of first dose and booster dose are fixed equal to E 1 = 30% and E 3 = 90%, respectively. The evaluation of the mean efficacy over the fully vaccinated individuals, E 2 , is instead more complicated due to different timing of full vaccination. The waning in time of different vaccines was evaluated for the US veterans in Ref. [53] . Using their data for the Pfizer-BioNTech vaccine (which is the most used in Italy in the first phase of the vaccination campaign) we find that these data are fitted with good approximation by the following function:

where e 1 , e 2 , t E and τ E are the parameters (shown in Tables 4 and 5), obtained fitting the data of Ref. [53] , and t f v is the timing of full vaccination (in [53] , the month when full vaccination is administered; here instead t f v is equal to 14 days after full vaccination). The time-depending mean efficacy of full vaccination is then evaluated averaging over the fully vaccinated individuals (from which the number of individuals, who received the booster dose, was suitably subtracted) as follows:

where N f v (t, t f v ) are individuals that received full vaccination at (t f v − 14) days and have not received yet booster dose at (t − 14) days. Following [31] we assume the reduction in vaccine effectiveness against Covid-19 infections to be primarily due to the waning immunity mechanism, rather than the Delta variant escaping vaccine protection. Thus we do not consider a further reduction factor due to new variants. Further parameters are listed below.

• χ j (with j = 1, 2, 3), the vaccination rates corresponding to first, second (or single dose for the Janseen vaccine), and booster doses. As discussed in the introduction these parameters changed discontinuously during the last year, thus we fix them as step functions through the best fits of the experimental data. In Fig. 8 the simulated evolution and the time series of vaccinated individuals, as reported by ISS, are compared.

• δ, the inverse mean latent period assumed to be the same both for vaccinated and unvaccinated people. However, the emergency of new variants has been typically characterized by a reduction of the latent and incubation period [32] . Therefore we modulate its value according to the dominant variant:

where ∆ 1 and ∆ 2 are fitting parameters, and t var and τ var are the same parameters present in Eq. (7).

• θ i A , θ i S and θ i D , with i = v, u, the recovery rates respectively for asymptomatic, symptomatic, diagnosed vaccinated and unvaccinated individuals. The recovery rates for asymptomatic and symptomatic undiagnosed individuals are assumed not dependent on time. Those of diagnosed individuals have a significant February 23, 2022 7/20 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 28, 2022. ;  variability over time, being affected by different factors, such as the test rate (increasing when the daily number of tests increases) and the number of active cases (decreasing when the sanitary system is overload) . We assume θ u D and θ v D to be functions of time, not depending on the vaccination status of individual, i.e. θ u D = θ v D = θ D , and following the form:

where θ 1 , θ 2 , t θ and τ θ are fitting parameters. This behavior corresponds to a decreasing in time of the number of days spent by infected diagnosed individuals in the diagnosed compartment, due to increased testing efficiency of the healthcare system through the involvement of analysis centers and pharmacies in testing operations.

• κ i D and κ i , with i = v, u, the mortality rates for diagnosed and undiagnosed infected individuals, respectively. The mortality rate for unvaccinated diagnosed individuals is assumed to follow the form:

where κ 1 , κ 2 , t mor1 , t mor2 , τ mor1 and τ mor2 are fitting parameters 3 . This behavior corresponds to a decreasing of the mortality during the summer period (from κ 1 + κ 2 to κ 1 − κ 2 ), probably due to the lower age of infected people, which returns to the same value of winter/spring in autumn.

Following the data of ISS, the mortality rate for vaccinated diagnosed individuals is assumed to be lower than for unvaccinated people. In particular we assume the κ v D = κ u D /11.5, according to the relative risk estimation among vaccinated and unvaccinated individuals reported in [54] 4 .

Finally the mortality rate, κ i , for both vaccinated and unvaccinated undiagnosed individuals, is assumed equal to zero;

• i , with i = v, u, corresponding respectively to the asymptomatic percentages of infections among vaccinated and unvaccinated individuals. Although there is evidence that vaccinated individuals are protected from severe illness, to our knowledge, a systematic study of differences in the asymptomatic fraction of infection between vaccinated and unvaccinated individuals is still lacking. For this reason, we fix the fraction of asymptomatic vaccinated and unvaccinated individuals to be equal.

• η i A and η i S , with i = v, u, corresponding respectively to the detection rates for asymptomatic and symptomatic infective individuals. These parameters are chosen to be constant in time and dependent on the vaccination status of individuals, both larger for unvaccinated than vaccinated people, reflecting the tendency of unvaccinated individuals to swab more easily than vaccinated ones. 3t is any time tmor 1 and tmor 2 , such that the two hyperbolic tangent functions assume the same asymptotic value. For simplicity we can putt = (tmor 1 + tmor 2 )/2). 4 According to [54], the relative mortality risk between vaccinated and unvaccinated individuals is 9.0, 11.5 and 30.3, respectively, for individuals fully vaccinated more than 120 days before, for those fully vaccinated less then 120 days before and for those with booster dose. During the period under investigation, the number of individuals with booster dose was negligible and mostly of the population recently completed the vaccination cycle. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 28, 2022. ;  • φ, corresponding to the rate at which recovered individuals loose immunity. On the time scales here considered, the presence of healed individuals, who return to being susceptible with a rate between 180 −1 and 270 −1 day −1 (as estimated in literature), turns out to be totally irrelevant, so we choose to set it equal to zero.

For the time dependence of the parameters in Eq.s (1) we preferred the use of hyperbolic tangents, except for Eq. (7), in which we adopted the exponential law -used also in Eq. (11) -obtained by fitting the real data for the prevalence of Delta variant in Italy. However, while in Eq. (9), the hyperbolic tangent comes up as a fit of the vaccine efficacy data, published in Ref. [53] , in all other cases, the hyperbolic tangent was preferred because it allows to describe the crossover between two different values of the parameter by means of a continuous function, with derivatives of any order that are continuous in turn.

The ODE Eqs. (1) have been solved using the SciPy libraries with initial conditions reported in Table 2 of Appendix. The best fit parameters, obtained minimizing the χ-square with respect to the experimental data, are listed in Tables 3, 4, 5 of Appendix. Besides the detected cases, our model predicts a relevant number of undetected symptomatic and asymptomatic cases, as shown in Fig. 2 (lower right) . According to our prevision, the percentage of undetected cases evolves from a minimum value of 3.6% at the beginning of summer, when the outbreak slowed down, to a value larger than 43% at the end of autumn. This circumstance can be related to the tendency within the large pool of vaccinated individuals to avoid testing when asymptomatic. This mechanism is also useful to understand the result on the relative incidence among unvaccinated and vaccinated individuals. As shown in Fig. 3 , our model predicts a different impact of the epidemic on vaccinated and unvaccinated individuals, in good agreement with data published by the ISS. In the inset of Fig. 3 , the simulated relative incidence among unvaccinated detected cases and vaccinated detected ones (blue line in figure) is compared with the same quantity evaluated on the total cases (both detected and undetected) (green line), and with the experimental incidence (red dots), clearly evaluated only on detected cases. The blue line is typically slightly higher than the red dots, and systematically higher than the green line. This last circumstance is consistent with the hypothesis that, differently from unvaccinated individuals that are frequently required to test, vaccinated individuals, specially if asymptomatic, may remain undetected more frequently than unvaccinated ones.

When comparing model predictions and epidemiological data, it must be taken into account that both are not error-free. The error on real epidemiological data is difficult to assess, due to the stochasticity inherent in the epidemic spreading (an effect February 23, 2022 9/20 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 28, 2022. ; Fig 3. Main: incidence of total active cases (I tot = I A + I S + I D ) among vaccinated and unvaccinated individuals. Inset: the simulated relative incidence between unvaccinated detected cases and vaccinated detected ones, I u D /I v D (blue line), is compared with the same quantity evaluated on both detected and undetected cases, I u tot /I v tot (green line), and with the experimental data (red dots), provided by the ISS.

completely neglected in a deterministic compartmental model in a finite population with short-range interactions), but also due to the method of data acquisition that can introduce random, but also systematic, errors (see for instance the discontinuity present in the experimental data in mid-June due to an incorrect communication of the healed individuals by some Italian regions). A source of error is certainly linked to the diagnostic capacity, which has changed significantly during the epidemic and which is mainly linked to the amount of swabs that can be done daily.

Concerning the simulations, the main source of error is the estimation of epidemic control parameters, but also of the initial conditions. In order to reproduce the experimental curves, we had to fix many different parameters, concerning the efficacy of vaccine, the waning immunity mechanism, the spreading of new variants, and other parameters such as those appearing in Eq.s (7) and (9) . It is reasonable to wonder how much the results obtained depend on the initial conditions and on these parameters set a priori. To explore these aspects, we carried out a sensitivity analysis focusing on one crucial number, such as the total number of cases diagnosed on the last day of the simulation, and evaluating how this number changes under modifications of the external parameters and initial conditions. We find that the results obtained are very stable for a fairly wide variation of the external parameters. In particular, Fig. 4 shows that the total number of detected cases at the final day of simulation (December 16, 2021) varies less than 1% by changing the vaccine efficacy parameters in Eq. (9) over a reasonable range of variability. Fig. 4 represents the total detected cases as a function of the maximum efficacy of full vaccination, e 1 + e 2 , and the medium time of antibody decay, t E . It shows that the same total incidence is obtained when a decline of efficacy is counterbalanced by a suitable growth of the waning time. A similar variation less than 1% is observed by fixing e 1 + e 2 and changing the decay interval, τ E , in Eq. (9), together with t E (data not shown).

Similarly it is interesting to explore what happens by varying the parameters regulating the insurgence of the Delta variant and the corresponding augmented transmissibility in Eq. (7) . In particular by varying the increase in transmissibility, ∆ σ , in the range [0. 4, 1] , the time of appearance of the new variant, t var , in the range [65, 105] and the time regulating the speed in the spreading of the new variant, τ var , in the range [10, 50] (Fig.s 5) , the total number of detected cases at the final day of simulation varies less than 1%. Furthermore, it is observed that for values of ∆ σ small enough, it is completely irrelevant to vary t var and τ var in the established ranges. Similar findings are obtained varying the percentages of asymptomatic infected individuals, v and u , or the initial conditions, in a reasonable range (data not shown).

In conclusion, the errors on the simulated data due to the estimation of initial conditions and of the parameters regulating efficacy of vaccine, waning immunity, spreading of the Delta variant and fraction of asymptomatic individuals, in realistic ranges of variation, do not seem to be greater than a few percent. As previously seen, there are many different mechanisms that contribute to the February 23, 2022 10/20 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 28, 2022. ; evolution of the epidemic, such as the mobility of the population, the perceived risk, the seasonal variation in the virus stability, the appearance of new variants and, obviously, the vaccination campaign. All these mechanism are encoded in the model, through appropriate terms. It is thus interesting to analyze the contribution of each term to the epidemic evolution. Fig. 6 shows the trend of diagnosed cases in different scenarios obtained turning off one term at a time in Eq. (2) and (4) (and only for the σ var also in Eq. (8)), or assuming the absence of vaccine. In particular, in Scenario I we disregard the effect of seasonal reduction of virus stability and infectivity during the summer period, in Scenario II we 'freeze' the awareness of individuals to the initial value, when the risk perception was quite high and the attention of individuals to respect sanitary measures, such as social distancing, frequent hand washing and wearing sanitary masks, was high as well. In Scenario III we consider what would have happened in absence of Delta variant. In Scenario IV we consider the effect of freezing the mobility to a low value, as the one during the lock-down period. Finally, Scenario V is devoted to understand what would have happened if there were no vaccines available.

For each scenario we evaluate the increase in the total infected detected cases and in total deaths with respect to the experimental initial values of ISS at time t 0 (∆I scen = I scen tot (t f ) − I ISS tot (t 0 ) and ∆D scen = D scen (t f ) − D ISS (t 0 )) and in Table 1 we report the percentage variations of these quantities in the specific scenario with respect to those of the reference simulation ((∆I scen − ∆I ref )/∆I ref and analogously for the death term). As already observed in [37] , we find that the seasonal modulation of virus transmissibility is an essential ingredient in order to explain the evolution of the epidemic curve during the summer period. Indeed by assuming this term to be constant (violet line in Fig. 6) , and equal to its value in winter time, an agreement between data and simulation would have been obtained only in the first peak, then the curve for diagnosed cases, after the achievement of a minimum at the beginning of June, would start to rise again, reaching a second and significantly higher peak in summer, followed by a decrease to a long and still high plateau. In such a circumstance at the end of the simulation period (December 16, 2021) the number of total infections and total deaths would have been significantly higher (175.1% and 95.7% respectively) than the observed ones.

Analogously, if the awareness mechanism hadn't been at work, a larger and higher peak would have been reached (blue line in Fig. 6 ), followed by a rapid decrease to zero of diagnosed individuals already around the beginning of September, so the subsequent increase in both the contact and seasonal terms would not have been sufficient to produce an increase in diagnosed cases. Let us note that this behavior corresponds to the situation in which the epidemic had spread always with the same high level of awareness that population had at the end of February, when Italy was in lock-down. But, in correspondence of the peak, the awareness is lower than in the reference simulation (due to the form of Eq. (5) and to the higher value of I D in the peak), while it is larger in the summer period. As a consequence, this Scenario, is characterized by a negative percentage variation in ∆I tot with respect to the reference simulation, while February 23, 2022 11/20 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted February 28, 2022. ;  the percentage variation in ∆D is positive due to the enlargement and rise of the first peak, in a phase of epidemic with higher mortality. Both contributions of the contact term, Eq. (3), and of the Delta variant, Eq. (7), have minor and very similar impacts on the epidemic evolution and their shutdown determines only the reduction of the summer peak and a slower autumn growth (brown and green lines in Fig. 6, respectively) .

A separate analysis deserves the Scenario V with no vaccines. What would have happened in absence of vaccines? Fig. 6 shows the time evolution of active cases in absence of vaccinations (orange curve). During the first peak the curve does not differ much from the experimental one, consistently with the fact that the number of vaccinated individuals was still contained in this time window. However, from the summer period, and specially in autumn, when the vaccination campaign reaches a significant percentage of the population, scenario V foresees a significant increase compared to the case of the reference simulation/experimental data. In particular according to our simulation, if vaccines were not available, at the end of simulation, one would have obtained a percentage increase of 63% in the total detected infections and of 55% in the total deaths. So far we have presented the model previsions up to mid-December. If we let the simulation run beyond this moment, the agreement between model predictions and real data becomes gradually worse. Actually this is not surprising because of the presence of the Omicron variant that became more and more important at the end of the year. On the other hand, the disagreement between the model and reality can give us a quantitative estimate of the presence of Omicron variant in Italy. In Fig. 7 the discrepancy between the diagnosed active cases provided by the ISS and the simulated data, I D , are plotted as a function of time. If we attribute the difference between observed and expected data to the presence of the Omicron variant, we can conclude that it appeared in Italy around the beginning of December, grew rapidly and became dominant in mid-January. This picture is in substantial agreement with the data provided by ISS [55], which gives prevalence estimates at the national level on January 3, 2022 for the Delta variant 19.22% (range: 0.0% − 66.7%) and for the Omicron variant 80.75% (range: 33.3% − 100%), despite the prevalence percentages are measured in a non-reproducible way within our model (the virus is sequenced by sample on positive swabs). The data in Fig. 7 may be indeed influenced by a different permanence of the individuals affected by the two variants of the virus in the compartment of the diagnosed individuals, for example if individuals with the Delta variant had on average a more severe disease, they would remain longer than the individuals affected by Omicron variant in the I D compartment, systematically affecting the prevalence of the two variants estimated in this way.

The model has some limitations. The first one is that the system is assumed to be closed and protected from the injection of new cases from abroad. This circumstance is not fully justified, specially in the summer period, when the touristic flows increase. However, according to data published by ISTAT (National Institute of Statistics), even if the international tourist flow in 2021 was in recovery with respect to the year 2020 (+22, 3%), it was still far from the levels of 2019 (−38, 4%) [56] . The extension of the model to open system is left to future work.

Secondly, it has been shown that the vaccine efficacy to protect against severe infections is higher than the efficacy against mild or asymptomatic infection. Our model doesn't distinguish the symptomatic cases according to the severity of symptoms, being pauci, mild and severe symptomatic cases, as well as hospitalized cases, all included in the symptomatic compartment. It would be interesting to rearrange the model in order to measure differences in hospitalization and severity of symptoms between vaccinated and unvaccinated individuals. However this would involve the introduction of new compartments into the model, circumstance that we avoided in order to keep the model simpler in the present paper.

In [57] authors study the attenuation of antibody titres after the second dose, showing that the most important factor in determining the waning immunity is sex, age and smoking. Our model does not take into account the age structure of population and not even the sex groups, and thus it is not able to capture differences in infections, mortality and recovering among different groups of individuals. We leave this interesting in-depth analysis to a future work.

Finally, it will be interesting to extend the stochastic models proposed in Ref.s [58, 59] for the pandemic H1N1 to the forth wave of Covid-19 epidemic: in this case indeed the almost total absence of mobility restrictions, which were instead present during the previous waves, would allow to apply the social contact hypothesis [60] in order to reproduce the epidemic spreading. 

In conclusion, we presented a deterministic mean field model in which the vaccinated compartments with different number of doses have been suitably introduced. The model is able to reproduce the epidemic spreading in Italy during the third and in part the fourth wave of Covid-19. The analysis of the ingredients that must be taken into account in order to reproduce the epidemiological curves teaches a lot about the disease. The strong seasonal trend of the epidemic has been confirmed, together with the role of awareness mechanisms that allow to mitigate the epidemic spreading through individual protective behaviors adopted when the risk perception increases. The effects of the appearance of the Delta variant and the contact increase during the summer months were instead only marginal, causing a slight rise in the summer peak, being the mitigating effect of summer temperatures stronger. The model also predicts with remarkable accuracy the appearance of the Omicron variant and its becoming dominant in January 2022. According to our model, in absence of a vaccination campaign, the total number of infections and deaths would have been dramatically higher, confirming the fundamental role of the vaccine in containing the pandemic and saving human lives. Table 5 . Time intervals and reference days in Eq.s (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) , obtained as fitting parameters. ∆ σ 0.60 transmissibility rise due to Delta variant, Eq.(7) v = u 0.65 asymptomatic infected individual percentage f 1 0.70 susceptibility reduction function after one vaccine dose f 3 0.10 susceptibility reduction function after booster φ 0 day −1 immunity lost rate Table 6 . Further parameters, whose values are fixed a priori. vaccinated individuals is negligible at t 0 , I v D 0 , I v A0 , I v S 0 and E v 0 are put equal to zero; the initial values for I u A0 , I u S 0 and E u 0 , which are not experimentally observables, are obtained extrapolating simulation results from previous model developed in [37] ; finally S 0 is obtained as difference between 1 and the sum of the initial values of all the other compartments.

Fitting rates present in Eq.s (1), whose value is not dependent on time, are given in Table 3 . Constants present in Eq.s (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) , obtained as fitting parameters, are listed in Table 4 . Days and intervals of time present in Eq.s (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) , obtained as fitting parameters are reported in Table 5 (in the simulations the time 0 corresponds to the initial time, t 0 reported in Table 2 , i.e. February 19, 2021) . Further parameters are reported in Table 6 . Fig. 8 shows the comparison between the simulated evolution and the time series of vaccinated individuals, as reported by ISS.

As already stressed in [37] most of the epidemiological compartmental models published in literature include more than three compartments and a relevant number of parameters that regulate fluxes between these compartments, within the specific model. Few available observables (i.e. data sets) do not allow to univocally fix the parameter values in the phase space. The basic idea at the base of such modelization is to exhibit a possible set of parameters, that, within the specific model/formulation, allows to reproduce the epidemic evolution.

Ministero della salute -Direzione generale della Prevenzione Sanitaria

Ministero della salute -Direzione generale della Prevenzione Sanitaria

EMA AstraZeneca's COVID-19 vaccine: EMA finds possible link to very rare cases of unusual blood clots

SARS-CoV-2 breakthrough infections in vaccinated individuals: measurement, causes and impact

A systematic review of methodological approaches for evaluating real-world effectiveness of COVID-19 vaccines: Advising resource-constrained settings

Efficacy and safety of the mRNA-1273 SARS-CoV-2 vaccine

Safety and Efficacy of the BNT162b2 mRNA Covid-19 Vaccine

Impact and effectiveness of mRNA BNT162b2 vaccine against SARS-CoV-2 infections and COVID-19 cases, hospitalisations, and deaths following a nationwide vaccination campaign in Israel: an observational study using national surveillance data

Prevention and Attenuation of Covid-19 with the BNT162b2 and mRNA-1273 Vaccines

Safety and Efficacy of the BNT162b2 mRNA Covid-19 Vaccine through 6 Months

Safety and efficacy of the ChAdOx1 nCoV-19 vaccine (AZD1222) against SARS-CoV-2: an interim analysis of four randomised controlled trials in Brazil, South Africa, and the UK

Oxford-AstraZeneca COVID-19 vaccine efficacy

Effectiveness of heterologous ChAdOx1 nCoV-19 and mRNA prime-boost vaccination against symptomatic Covid-19 infection in Sweden: A nationwide cohort study. The Lancet Regional Health

Understanding COVID-19 vaccine efficacy

Influence of age on the effectiveness and duration of protection of Vaxzevria and CoronaVac vaccines: A population-based study. The Lancet Regional Health

Associations of the BNT162b2 COVID-19 vaccine effectiveness with patient age and comorbidities at daily resolution

Africa needs more genome sequencing to tackle new variants of SARS-CoV-2

Coronavirus variants are spreading in India -what scientists know so far

Dynamics of antibody response to BNT162b2 vaccine after six months: a longitudinal prospective study. The Lancet Regional Health

COVID vaccine immunity is waning -how much does that matter?

Epidemic scenarios of Delta variant in France in the summer 2021

Effectiveness of Covid-19 Vaccines against the B.1.617.2 (Delta) Variant

Evaluation of a Booster Dose (Third Dose). Vaccines and Related Biological Products Advisory Committee Briefing Document FDA.gov

Safety and immunogenicity of SARS-CoV-2 variant mRNA vaccine boosters in healthy adults: an interim analysis

SARS-CoV-2 neutralization with BNT162b2 vaccine dose 3

Effectiveness of mRNA BNT162b2 COVID-19 vaccine up to 6 months in a large integrated health system in the USA: a retrospective cohort study

Transmission, viral kinetics and clinical characteristics of the emergent SARS-CoV-2 Delta VOC in Guangzhou

Duration of Protection against Mild and Severe Disease by Covid-19 Vaccines

Waning Immunity after the BNT162b2 Vaccine in Israel

Waning of BNT162b2 vaccine protection against SARS-CoV-2 infection in Qatar

Waning immune humoral response to BNT162b2 covid-19 vaccine over 6 months

Beyond the peak: A deterministic compartment model for exploring the Covid-19 evolution in Italy

Effect of Vaccination on Household Transmission of SARS-CoV−2 in England

Impact of BNT162b2 vaccination and isolation on SARS-CoV-2 transmission in Israeli households: an observational study

Vaccination with BNT162b2 reduces transmission of SARS-CoV-2 to household contacts in Israel

Water, climate change, and Covid-19: prioritising those in water-stressed settings

Evidence that high temperatures and intermediate relative humidity might favor the spread of Covid-19 in tropical climate: a case study for the most affected Brazilian cities

Effective transmission across the globe: the role of climate in Covid-19 mitigation strategies. Lancet Planet Health

Initial report of decreased SARS-CoV-2 viral load after inoculation with the BNT162b2 vaccine

Community transmission and viral load kinetics of the SARS-CoV-2 delta (B.1.617.2) variant in vaccinated and unvaccinated individuals in the UK: a prospective, longitudinal, cohort study

Lattice model for influenza spreading with spontaneous behavioral changes

Impact of temperature on the dynamics of the COVID-19 outbreak in China

Temperature dependence of COVID-19 transmission

Association between meteorological factors and daily new cases of COVID−19 in 188 countries: A time series analysis

Effectiveness of a Third Dose of BNT162b2 mRNA Vaccine

SARS-CoV-2 vaccine protection and deaths among US veterans during 2021

Comuniato Stampa Movimento Turistico In Italia -gennaio/settembre

Attenuation of antibody titres during 3-6 months after the second dose of the BNT162b2 vaccine depends on sex, with age and smoking as risk factors for lower antibody titres at 6 months

A Lattice Model for Influenza Spreading

Multiple Lattice Model for Influenza Spreading

A.L. and A.F. acknowledge financial support of the MIUR PRIN 2017WZFTZP "Stochastic forecasting in complex systems". The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

The values of R 0 , D 0 , V 10 , V 20 , V 30 , reported in Table 2 , correspond to the current values at the beginning of the period under examination, t 0 ; since the number of