Study on the effects of the restrictive measures for containment of the COVID-19 pandemic on the reproduction number $R_t$ in Italian regions

Abstract

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 $R_t$, we estimate how much different restrictive measures affect $R_t$, and we quantify the combined effect of the diffusion of virus variants and the beginning of the vaccination campaign upon the $R_t$ 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 $R_t$.

1 Introduction

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 $6{\mathrm {th}}$, 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 $30{\mathrm {th}}$, 2020 [2]. The color assignments to the Italian regions along the period November $6{\mathrm {th}}$, 2020–April $26{\mathrm {th}}$, 2021, are displayed in Fig. 1. We note that the change of NPI in this period was quite frequent. For instance, Lombardia, the most populated Italian region with about 10.0 million inhabitants, changed level 15 times in about six months.

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.

Fig. 1

NPI periods in the Italian regions, characterized by the colored markers white, yellow, orange and red, in the November $6{\mathrm {th}}$, 2020–April $26{\mathrm {th}}$, 2021, period of time. Besides the nine regions previously mentioned, it is also reported Sardegna, the only Italian region to reach the white restriction in the time period considered

2 $R_t$ and its synchronization with the new COVID-19 cases

2.1 Computation of $R_t$

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 $\lambda (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.¹ The fit to an exponential curve is performed within moving windows of 14 days.

$R_t$ is determined from $\lambda =\lambda (t)$ using Eq. 1:

$$\begin{aligned} R_t = (1+\lambda \theta )^\kappa \,\,, \end{aligned}$$

(1)

where $\theta$ and $\kappa$ 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 $\theta$ and $\kappa$ of the generation time distribution are taken from [9].

2.2 Synchronization of $R_t$

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 $\lambda$ 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 $\lambda$ 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:

$$\begin{aligned} g(t;a,k_G,T_p) = a \cdot k_{G} e^{-e^{-k_{G}(t-T_p)}-k_{\tiny {G}}(t-T_p)}\,\,\,. \end{aligned}$$

(2)

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$.

Fig. 2

Distribution of the number of new daily cases in Lombardia (squares) together with their smoothing (blue line) and the fit to the sum of three Gompertz functions (orange line). Dotted curves indicate the three single Gompertz components called G1, G2, G3. The $R_t$ values, as computed with the Covidstat algorithm (green curve) and with the Cori algorithm (brown curve) are superimposed with a different scale indicated on the right of the plot. The dashed horizontal orange line is drawn to mark the critical value of $R_t$ =1

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].

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.

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

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.

2.3 $R_t$ data

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.

Fig. 3

$R_t$ values as a function of time from November $6\mathrm {th}$, 2020, to April $26\mathrm {th}$, 2021, in the nine considered regions. The colors of the dots correspond to the three levels of the NPI (yellow, orange, red) locally enforced in different periods of time

3 $R_t$ models

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$, $\hat{R}_t$, for all the time periods in which the NPI remained unchanged. At the beginning of each new time period we set $\hat{R}_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:

• Model-a

  $$\begin{aligned} \hat{R}_t (t+t_{D} +1) = R_t (t+t_{D}) \cdot \alpha _i \end{aligned}$$

  with i=yellow, orange, red. In this model, the effect of the NPIs consists in a multiplicative factor $\alpha$ 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 $\alpha$ is greater than one) or to zero (when $\alpha$ is smaller than 1).

• Model-b

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

  $$\begin{aligned} \hat{R}_t (t+t_{D} +1) = R_t (t+t_{D}) + \frac{R_t (t+t_{D})-R^*_{t,i}}{N} \end{aligned}$$

  The parameter N defines the number of steps needed to reach $R_t^*$.

• Model-c

  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:

  $$\begin{aligned} \hat{R}_t (t+t_{D} +1) =R_t (t+t_{D}) + \frac{R^{'}_t(t+t_{D}) + \beta (R_t (t+t_{D})-R^*_{t,i})}{N} \end{aligned}$$

  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 $\beta$ is an additional degree of freedom in the fits.

We determine estimates of the parameters $\alpha _i$, $R_t^*$, $t_{D}$, N, $\beta$ by minimizing a $\chi ^2$ defined globally over all considered regions, as described in the next section.

3.1 The time delay between NPI changes and their effect upon $R_t$

The delay time $t_{D}$ necessary to NPI changes to produce modifications in $R_t$ trends, is determined by a scan on the overall $\chi ^2$ computed for all the regions assuming Model-c in the period of time November $6{\mathrm {th}}$, 2020–January $15{\mathrm {th}}$, 2021. The $\chi ^2$ is defined as

$$\begin{aligned} \chi ^2 = \sum _{t=1}^{T}{\left( \frac{R_t (t+t_{D}) - \hat{R}_t (t+t_{D})}{\sigma _{R_t (t+t_{D})}}\right) ^2} \end{aligned}$$

(3)

where t indicates the calendar days, T is the total number of days in the considered period of time, $\sigma _{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 $15{\mathrm {th}}$, 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 $\chi ^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 $\chi ^2$ distribution is quite broad and irregular; on the other hand, the time delay is expected to have a quite broad distribution. Two contributions to the dispersion of $t_{D}$ are the finite width of the generation time distribution, and the collection and processing time of molecular swabs that introduce a dispersion of 7.7 days as we estimated in [6].

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.

Fig. 4

$\chi ^2$ scan of the fit to $R_t^*$ with Model-c as a function of the delay time $t_{D}$

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 $1{\mathrm {st}}$, 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.

Fig. 5

$R_t$ trend in Sardegna. The change to the white level was established on March $1{\mathrm {st}}$, 2021. Considering the required delay of 8 days, we display the change on March $9{\mathrm {th}}$, 2021. The steep ascent of $R_t$ begins on March $14{\mathrm {st}}$, 2021. 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 after having applied the 8 days time delay

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

4 Effects of the three different NPI levels

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.

Fig. 6

Results of the fits for the parameters of the Model-a (upper left plot; what is reported in dots are the fitted estimates of the parameter $\alpha$, Model-b (upper right plot; the fitted estimates of $R_t^*$ are reported as squares) and Model-c (lower plot; the fitted estimates of $R_t^*$ are reported as squares) for all the regions and their combination (Italy columns). The colors of the fits, yellow, orange and red, correspond to the three NPI levels. Error bars are smaller than the markers and masked behind them. The horizontal blue dotted lines mark the threshold values $\alpha =1$ and $R_t^* =1$

Table 2 Fit values of the Model-a $\alpha _i$ parameters ($i=\hbox {yellow}$, orange, red) and of the $R^*_{t,i}$ parameters of Model-b and Model-c

The overall $\chi ^2/{\mathrm {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 $\chi ^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 $\chi ^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 values $\hat{R}_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 $\chi ^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$.

Fig. 7

Estimated $\hat{R}_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 value $\hat{R}_t (t)$ is set equal to $R_t$ (t) at any change of NPI level

In Fig. 8 we display the fits to $R_t$, $\hat{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 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.

Fig. 8

$R_t$ values as function of time in the nine considered regions. The colors of the dots correspond to the three NPI levels, yellow, orange, red after having applied a delay of $t_{D} =8$ days. 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. The blue lines represent the output of the fits, $\hat{R}_t$, computed with Model-c for each individual region

5 Effects of the vaccination campaign and virus variants

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 $15{\mathrm {th}}$, 2021, was 4.3% (0%) with the first (second) dose, respectively, while, as April $26{\mathrm {th}}$, 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 $18{\mathrm {th}}$, 2021, while this percentage was 17.8% at February $4{\mathrm {th}}$, 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 $6{\mathrm {th}}$, 2021–January $15{\mathrm {th}}$, 2021, and January $16{\mathrm {th}}$, 2021–April $26{\mathrm {th}}$, 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.

Table 3 Fit values of the asymptotic values $R_t^*$ for the three NPIs levels yellow, orange and red, computed with Model-c (see text) in two different periods of time

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.

6 Conclusions

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 $6{\mathrm {th}}$, 2020, to April $26{\mathrm {th}}$, 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 $\beta$) 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.

Appendix

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 $20{\mathrm {th}}$, 2020 [17], and March $13{\mathrm {th}}$, 2021 [18], nursery and elementary schools had been always open, except the period March, $3{\mathrm {rd}}$, 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 $6{\mathrm {th}}$, 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