key: cord-0976257-wjeiqqm4
authors: Bonifazi, Gianluca; Lista, Luca; Menasce, Dario; Mezzetto, Mauro; Pedrini, Daniele; Spighi, Roberto; Zoccoli, Antonio
title: Study on the effects of the restrictive measures for containment of the COVID-19 pandemic on the reproduction number [Formula: see text] in Italian regions
date: 2021-12-03
journal: Eur Phys J Plus
DOI: 10.1140/epjp/s13360-021-02076-6
sha: 3460fadef7a53badd487a7fa8f17692c4dc4ad59
doc_id: 976257
cord_uid: wjeiqqm4

Since November 6th, 2020, Italian regions have been classified according to four levels, corresponding to specific risk scenarios, for which specific restrictive measures have been foreseen. By analyzing the time evolution of the reproduction number [Formula: see text] , we estimate how much different restrictive measures affect [Formula: see text] , and we quantify the combined effect of the diffusion of virus variants and the beginning of the vaccination campaign upon the [Formula: see text] trend. We also compute the time delay between implementation of restrictive measures and the resulting effects. Three different models to describe the effects of restrictive measures are discussed and the results are cross-checked with two different algorithms for the computation of [Formula: see text] .

The main purpose of this paper is to quantify the effects on the reproduction number R t of the Non-Pharmaceutical Interventions (NPI) put in place in the Italian regions since November 6th, 2020.

Italian regions and the autonomous provinces of Trento and Bolzano have been classified into the four aforementioned colored levels: red, orange, yellow and white in decreasing levels of restrictions, corresponding to four risk scenarios, for which specific restrictive measures were foreseen, as reported in Appendix. The scenarios are described in [1] and risk indicators are defined in the Ministerial Decree of April 30th, 2020 [2] . The color assignments to the Italian regions along the period November 6th, 2020-April 26th, 2021, are displayed in Fig. 1 . We note that the change of NPI in this period was quite frequent. For instance, Lombardia, a e-mail: mauro.mezzetto@pd.infn.it (corresponding author) Three possible models to evaluate the effects of the NPI on R t will be compared, and we will address two aspects that affect its computation:

• how can R t be properly synchronized with the daily number of positive swabs, from which it is computed; • the value of the time delay between the implementation of a new level of NPI and its effects on R t .

We will compare the effects of the different NPI in the nine most populated Italian regions: Lombardia, Lazio, Campania, Sicilia, Veneto, Emilia-Romagna, Piemonte, Puglia, Toscana (in decreasing level of population). All together, these regions contain 80% of the total Italian population. We excluded regions with less than 2 million inhabitants, in order to ensure a sufficiently large number of cases and allow a precise R t computation.

In early January 2021, the vaccination campaign started in Italy, and in the same period novel virus variants spread across the country. These two important factors potentially affect R t in opposite ways and we will compute their combined effect.

We compute R t from the number of new COVID-19 daily cases reported by the Italian Dipartimento della Protezione Civile [3] . In reference [4] (CovidStat from here on) we developed a simple, yet precise, algorithm to compute R t . We publish daily real-time estimates of R t , together with other indicators of the development of the Italian outbreak in [5] .

The official R t used by the Italian Government is computed from the date of onset of symptoms, not from the date of the positive swab as has been done in this paper. We verified in [6] that the two methods converge to essentially identical values of R t (see in particular Fig. 6 of reference [6] ), once the two estimations of R t are synchronized. Furthermore, since the R t used by the Italian Government is computed with the Cori et al. algorithm [7] (Cori from here on), we verified in [4] and [6] that our Covidstat algorithm computes R t in full agreement with it. In any case most of the results of this paper will be cross-checked with the Cori algorithm, by using the EpiEstim package [8] and assuming the same modeling of the generation time distribution as measured in [9] .

We prefer to keep our Covidstat algorithm as the default for the following computations because it's much faster to compute (as discussed in [4] ) and because we have full control of its computation details in all the necessary steps, in particular when we proceed to synchronize R t with the data.

The R t computation with the Covidstat algorithm requires a smoothing of the original data as a preliminary step and, subsequently, an exponential fit to the smoothed data, determining in this way the time-dependent growth rate λ(t) of the number of new daily cases. The smoothing is performed applying the Savitzky-Golay Algorithm [10] in a time window of 21 days and with a second-order polynomial. 1 The fit to an exponential curve is performed within moving windows of 14 days.

R t is determined from λ = λ(t) using Eq. 1:

where θ and κ are, respectively, the scale and shape parameters of the gamma distribution that models the distribution of the generation time, the time interval between infector and infected pair [4] . Estimates of the parameters θ and κ of the generation time distribution are taken from [9] .

To study the time-dependent effects of NPI upon the trends of R t , it is important to assess the correct synchronization of the latter with the daily number of new cases. In real-time applications, it's customary to assign both the smoothing and the resulting λ of the exponential fit to the last day of the time sequence. In this way, the value of R t can be assigned to the latest day. A more reasonable option for synchronization would be to assign the smoothing and the parameter of the exponential fit to the central day of the sequence.

The way the CovidStat algorithm is defined, provides an absolute mean to determine the optimal synchronization of R t : when the new cases are at a maximum or at a minimum, the growth rate λ becomes zero and, according to Eq. 1, R t must therefore be equal to one. This feature has the additional advantage of being independent from the assumed distribution of the generation time.

In order to verify the synchronization of R t with the number of new daily cases, we determine the positions of the three prominent peaks in the period of time considered in this paper. To model them, we perform a fit to the sum of three Gompertz derivatives. Gompertz curves are widely and commonly used to describe the development of pandemics, as discussed for instance in [11] . Each derivative is defined as:

In this parametrization, a is the integral of the curve, T p it's the time when the peak is reached and k G is a growth-rate coefficient, modeling the shape, see also reference [11] .

To illustrate the procedure, Fig. 2 shows the distribution of the number of new daily cases in Lombardia. A smoothed distribution and the fit of the sum of Gompertz functions are superimposed on the data. We note the large dispersion of the number of new daily cases to be well beyond their expected Poisson statistical fluctuation. This dispersion makes the smoothing of the data necessary and sometimes produces false small secondary peaks that in turn are reflected in the evolution of R t . Table 1 displays the day in which the fitted curve reached the first peak in the nine Italian regions along with the dates when R t crossed the critical value of 1. The curves in Fig. 2 and the data in Table 1 show a good agreement between the dates in which a maximum of the new cases distribution is reached and the dates when R t = 1. However, we had to assign the R t values computed with Cori to the initial day of the 7-days time window under which it is computed, in contrast with the indications of [7] . Table 1 Fit values of the peak position of the second wave in nine Italian regions and date when R t curve crosses the R t = 1 value, according to the Covidstat and the Cori algorithms, obtained with a straight line fit to the five R t values around R t =1

Peak of new cases

The difference of times t between the peak of the new cases and the day when R t = 1 are also reported. Days are counted since February 24th, 2020, day 260 corresponds to November 9th, 2020

By adopting this synchronization (for the Covidstat algorithm the smoothing and the exponential fits are referred to the central day of the considered time window, while for Cori R t it is referred to the first day of the time window), the Covidstat R t can only be computed up to 8 days before the last day of available data, while Cori can be computed up to 7 days before the last day.

Data in Table 1 indicate a residual mismatch of a couple of days on average between the day in which a maximum is reached and R t crosses the value of 1. We decided not to adjust the synchronization for this amount of time, because it could be partially due to a bias induced by the fit to the Gompertz functions; in any case the procedure detailed in Sect. 3.1 guarantees an overall correct synchronization.

Once R t is synchronized with the number of new daily cases, we can proceed to plot the R t trends for the nine Italian regions considered. Data are displayed in Fig. 3 , and the colors of the dots correspond to the three NPI levels: yellow, orange and red. None of the considered regions has ever been classified as white, so we can't compute the effects of the white NPI level. We observe that all regions have the same descending trend at the beginning of the analyzed time range, because they were all subjects, at the time, to the same common NPI. Subsequently, the trends began to differ from region to region, because of the uneven NPI levels and their enforcement date.

We have to consider several factors that can influence the R t behavior beyond what is induced by the three NPI levels. Noticeably:

(a) the contagion tracking and screening procedures differ from region to region; (b) individual regions may have adopted specific additional restrictive provisions to those foreseen nationwide; in some occasions, limited parts of the regions, as single municipalities and provinces, were subject to stricter NPI; (c) some degree of freedom about the organization of schools of any degree along with Universities had been deferred to local level, more information about these points is reported in Appendix. Schools closures and partial openings, therefore, did not follow the NPI levels in the same homogeneous way. Protocols for tracking students and schools staff and the actions to be taken in case of positive cases during the school openings were left to decisions taken at regional level. (d) mobility inside each region could significantly differ because of different work activities and different network of mobility connections; (e) the vaccination campaign started in January 2021 with a potential impact on R t development; (f) virus variants appeared during the same period of time and became dominant towards its end with a potential impact on R t development.

For all these reasons, we believe that the effects induced by the NPI alone are not sufficient by themselves to fully account for the trends of R t . On the other hand, the purpose of this work is not to introduce a model capable of describing the development of the pandemic in all its aspects; we are rather interested in extracting the effects of the individual NPI.

We introduce three different models to describe the effects of the NPI on the time evolution of R t . They compute a predicted value for R t ,R t , for all the time periods in which the NPI remained unchanged. At the beginning of each new time period we setR t = R t . The time periods do not begin on the calendar day in which the NPI changed, t, but at the later day t + t D , where t D is the time required for the effects of R t to become apparent. t D is an important parameter of this analysis and will be discussed in paragraph 3.1. The models we take into consideration are:

with i=yellow, orange, red. In this model, the effect of the NPIs consists in a multiplicative factor α on R t . This method has been adopted for instance in [12, 13] . Despite its popularity in literature, this model has severe limitations. The most important one is that the values of R t can only asymptotically tend to infinite (when α is greater than one) or to zero (when α is smaller than 1).

To overcome the limitations of Model-a, we introduce a model where R t can asymptotically converge to any finite value R * t,i :

The data of Fig. 3 do not show discontinuities on the first derivative of R t , R t , in the regional trends. To guarantee the continuity of R t we additionally refine Model-b in the following way:

where the effects of the NPIs on the trend of R t are mediated by the time derivative R t (t + t D ); the weight β is an additional degree of freedom in the fits.

We determine estimates of the parameters α i , R * t , t D , N, β by minimizing a χ 2 defined globally over all considered regions, as described in the next section.

The delay time t D necessary to NPI changes to produce modifications in R t trends, is determined by a scan on the overall χ 2 computed for all the regions assuming Model-c in the period of time November 6th, 2020-January 15th, 2021. The χ 2 is defined as

where t indicates the calendar days, T is the total number of days in the considered period of time, σ R t is the error on R t . We have chosen Model-c because, as will be discussed in the following, it is the algorithm that best reproduces the data. Moreover, we limited the period of time to January 15th, 2021, because, as also discussed in the following, it is the most stable period of time before the introduction of the vaccination campaign and the spread of virus variants.

The χ 2 scan reported in Fig. 4 , obtained by varying the value of the time delay t D in the range from 0 to 22 days, indicates that an absolute minimum is reached for a delay of t D = 8 days. The χ 2 distribution is quite broad and irregular; on the other hand, the time delay is Eur. Phys. J. Plus (2021) 136:1208 Since we assigned the Covidstat R t to the center of the time interval of the fit (8 days before the last day) and the Cori R t to the start of the time interval (7 days before the last day), the value of t D = 8 days implies that changes on NPI can be detected not earlier than 16 days (Covidstat) and 15 days (Cori) after their enforcement.

As a cross-check, we display the trend of R t in Sardegna, which is not one of our selected regions, but it is the only one that achieved a white color assignment on March 1st, 2021, a decision which resulted in a pronounced rise of R t . The R t trend for Sardegna is shown in Fig. 5 and is in good agreement with our evaluation of t D = 8 days.

With an analogous χ 2 scan, the number of steps N for Model-b and Model-c is fixed to N = 9.

We estimate the parameters of each model for all the individual Italian regions and for their combination (we will call it Italy from here on). As a cross-check, we perform all the calculations by using both the Covidstat and Cori algorithms. The results are displayed in Fig. 6 . The values for the Italian sample are reported in Table 2 .

The overall χ 2 /ndof values are 33.4, 25.8, 20.0 for Model-a, Model-b, Model-c, respectively, as computed with the CovidStat R t algorithm for Italy and 36.2, 28.4, 24.3 for the same models, using the Cori R t algorithm. None of the χ 2 is fully satisfactory, but, as already discussed in Sect. 2.3, these models cannot claim to completely describe the course of the pandemic.

It's clear from the comparison of the χ 2 values that Model-a does not provide a fit quality to the data comparable to Model-b and Model-c. The latter provides the best fit to the data and in the following we will consider it as our reference. In all the Italian regions the red set of NPI turns out to provide better containment than the orange set and this latter is better than the yellow set, well beyond the errors of the fits. Absolute numbers in the Italian regions can differ (as discussed in Sect. 2.3, different regions implemented additional and localized containment measures beyond the NPI considered here), but the average values in Italy are in good agreement with the regional trends.

The results indicate that the yellow NPI level has been, in general, inadequate to guarantee the containment of the pandemic, since its fitted R * t value is bigger than 1 for all cases (with the only exception of Lazio with an estimated R * t = 1.0). The orange NPI level produces trends dangerously close to R * t = 1, while the red NPI level guarantees a good containment of the pandemics. We recall that, as reported in Appendix, the main additional restrictions set by the red level were: further restrictions to mobility, closure of all the stores with few exceptions and suspension of all sport competitions. Parameters computed with the CovidStat and Cori algorithms are in general in good agreement.

We display in Fig. 7 the fitted valuesR t computed with the three models superimposed to the R t data in the case of Lombardia. A visual inspection of Fig. 7 confirms the outcome of the χ 2 of the fits: Model-c is the model that better follows R t trends. It isn't fully satisfactory, but we have discussed in Sect. 2.2 that some features of the R t trends are due to fluctuations of the data on positive swabs and we discussed in Sect. 2.3 that the NPI we are modeling here are not the only factors that can influence the trends of R t .

In Fig. 8 we display the fits to R t ,R t computed with Model-c for all the regions. The colors of the R t dots in Fig. 8 are different from those of Fig. 3 because here we have applied the delay time t D . The R t trends differ significantly from region to region, it should be reminded Fig. 7 EstimatedR t values in Lombardia computed with the Model-a (grey dotted curve), Model-b (pink dotted curve) and Model-c (blue curve) compared with the measured values of R t (colored dots). The colors of the R t dots correspond to the three NPI levels yellow, orange, red. They are delayed by 8 days as computed in the text. For this reason, the first 8 values of R t are colored in grey since they don't correspond to any of the three NPI levels here considered. We remind the fitted valueR t (t) is set equal to R t (t) at any change of NPI level that following the color assignments of Fig. 1 the regions had different NPI assignments at different times. The resulting phenomenology results quite complicated. It's significant that Model-c on the average can follow the trends of R t in all the Italian regions.

During the period of time considered in this study, two new important facts happened that could potentially impact on the trends of R t . The vaccination campaign began in Italy early in January 2021. According to the data reported in the COVID-19 Opendata Vaccini repository [14] , the fraction of the Italian population vaccinated by January 15th, 2021, was 4.3% (0%) with the first (second) dose, respectively, while, as April 26th, these fractions were increased to 20.1% (8.3%).

On the other hand, virus variants (variants of concerns or VOC), specifically lineage B.1.1.7 (English variant), began to spread out in Italy in the same period: according to [15] , 86.7% of the new cases in Italy were due to this variant at March 18th, 2021, while this percentage was 17.8% at February 4th, 2021 [16] .

Both VOC and vaccinations have a potentially important impact on the spread of the contagion; therefore we performed a fit with Model-c in two separate period of times: November Eur. Phys. J. Plus (2021) 136:1208 6th, 2021-January 15th, 2021, and January 16th, 2021-April 26th, 2021, to measure a possible overall effect of these two new factors in the progress of the pandemics. The results of the fits are reported in Table 3 . The combined effect of virus variants and the vaccination campaign in the considered period of time, resulted in a worsening of NPIs effectiveness. The yellow NPI level, for instance, worsened its effectiveness in containing the spread of contagion by 15%. After this period of time, the effect of the variant B.1.1.7 reached its maximum, while the vaccination campaign continues. So, hopefully, if no further variants will influence the R t trends, the effect of the considered NPI will improve.

The goal of this paper is the computation of the effectiveness of the three NPI called yellow, orange and red introduced in the Italian regions from November 6th, 2020, to April 26th, 2021.

The task to single out the effects of the three NPI is particularly complicated because, as we discussed in the text, several other pharmaceutical and nonpharmaceutical measures influenced the trends of R t in the same period. We had been facilitated by the fact that the changes of NPI happened several times in the period taken into consideration, often at different times in different regions. Furthermore, we introduced several original procedures to improve the analysis capability of separating the effects of the NPI from the other sources.

We firstly synchronized the R t values with the dates in which the NPI changed in the Italian regions. The synchronization is a necessary step for any procedure that tries to compute changes in R t trends due to changes of NPI at given calendar dates. To this purpose we introduced two original methods, the first of which exploits the fact that our original Covidstat algorithm directly connects R t with the growth rate of an exponential fit.

We identified critical issues in the model often used in literature to quantify the effects of NPI measures. We introduced two methods to overcome them and the fit to the data shows that the two new models clearly outperform the original one. We selected the so-called Model-c as a reference for all the results of the paper.

We demonstrated that our model is able to describe with just 4 parameters (3 asymptotic values R * t and a weight parameter β) the development of the pandemics in the 9 Italian regions considered in this paper. The development of the pandemics in the Italian regions is quite differentiated because they underwent different NPI restrictions at different times. The agreement between the model and the data is not complete because several local restrictive measures had been taken and because the poor quality of the positive swabs data introduces important noise in the R t trends. We discussed these features in the paper.

The model provides a coherent picture wherein all the Italian regions the red set of NPI results to have smaller R * t values of the orange set and this latter has always R * t smaller than the yellow set, well beyond the statistical errors. We find, that averaged over the nine most populated Italian regions, only the red NPI level can produce a significant reduction of R t , below the threshold value of 1. The orange level keeps R * t around 1, while the yellow level targets R t values above one.

We measured a worsening of the containment effects in the period of time when the variant of concern lineage B.1.1.7 spread out in the country, even if, in the same period, the vaccination campaign started.

All the computations have been performed with two R t algorithms: the Covidstat algorithm that we developed [4] and the widely used standard Cori et al. algorithm [7] . We found a good agreement between the results obtained using the two algorithms. Table 4 continued

All the restrictions of Yellow, plus:

• Prohibited movement in and out of regions and municipalities (with some exceptions) • Coffee shops and restaurants closed 7 days a week.

Take-away allowed until 10:00PM. Home deliveries without restrictions

All the restrictions of Orange, plus:

• It is forbidden any movement, even within its own municipality, except for reasons of work, needs and health. • Closed all stores, except supermarkets, grocery stores, newsstands, tobacconists, pharmacies, laundries, hairdressers. • Suspended all sports competitions except those recognized as being of national interest. Closed sports centers.

Prevention and response to COVID-19: evolution of strategy and planning in the transition phase for the autumn-winter season, Ministero della Salute and Istituto Superiore Sanità

Dati COVID-19 Italia

A simplified estimate of the Effective Reproduction Number Rt using its relation with the doubling time and application to Italian COVID-19 data

A study on the possible merits of using symptomatic cases to trace the development of the COVID-19 pandemic

A new framework and software to estimate timevarying reproduction numbers during epidemics

EpiEstim: Estimate Time Varying Reproduction Numbers from Epidemic Curves

The early phase of the COVID-19 outbreak in Lombardy, Italy

Smoothing and differentiation of data by simplified least squares procedures

The use of Gompertz models in growth analyses, and new Gompertz-model approach: An addition to the Unified-Richards family

Imperial College COVID-19 Response Team, Estimating the effects of nonpharmaceutical interventions on COVID-19 in Europe

Variant Of Concern) in Italia: lineage B.1.1.7, P.1, P.2,lineage B.1.351, lineage B.1.525

Prevalenza della variante VOC 202012/01, lineage B.1.1.7 Italia

We acknowledge the effort of ISS of making public the data of the symptomatic cases of COVID-19. The present work has been done in the context of the INFN CovidStat project that produces an analysis of the public Italian COVID-19 data. The results of the analysis are published and updated daily on the website https://covid19.infn.it/. The project has been supported in various ways by a number of people from different INFN Units. In particular, we wish to thank, in alphabetic order: Stefano Antonelli (CNAF), Fabio Bredo (Padova Unit), Luca Carbone (Milano-Bicocca Unit), Francesca Cuicchio (Communication Office), Mauro Dinardo (Milano-Bicocca Unit), Paolo Dini (Milano-Bicocca Unit), Rosario Esposito (Naples Unit), Stefano Longo (CNAF), and Stefano Zani (CNAF). We also wish to thank Prof. Domenico Ursino (Università Politecnica delle Marche) for his supportive contribution.

The datasets used for this analysis are quoted in bibliography. Some of the article plots are available under https://covid19.infn.it/ and can be downlowded in csv format. All the remaining material is available from the corresponding author on reasonable request.]

A brief description of the rules in place for the four NPI levels denominated as White, Yellow, Orange, Red is reported in Table 4 .Following the Ministerial Decrees of October 20th, 2020 [17] , and March 13th, 2021 [18] , nursery and elementary schools had been always open, except the period March, 3rd, 2021-April 2nd, 2021, when they closed. High schools and Universities were suggested to organize distance teaching, with some degree of freedom at local level. On April 6th, 2021, it has been decided to bring the distance teaching at 50% and on April 26th, 2021, this percentage will become 70% to 100%. middle schools (three years, 12 to 14 years old students) had different implementations in the different regions, the most common was to follow for the first year the rules of elementary schools and for the following two years the rules of high-schools. Table 4 Brief description of the rules in place for the four NPI levels denominated as White, Yellow, Orange, Red

• Requirement to wear masks outdoors• Social distancing of 1 meter• Curfew: decision taken at the regional level, Sardegna chose a curfew from 11:30PM to 5:00AM.

All the restrictions of White, plus:• Curfew from 10:00PM to 5:00AM.• Shopping centers closed on holidays.• Museums and exhibitions closed.• Distance education for high schools; in-person teaching for preschools, elementary schools, and middle schools. Universities closed, except for laboratories • Public transportation must travel with max 50% occupancy • Coffee shops and restaurants open until 6:00PM. Takeaway allowed until 10:00PM. Home deliveries without restrictions • Suspension of activities of amusement arcades, betting shops. • Pools, gyms, theaters, and movie theaters remain closed. Open sports centers.