key: cord-0466260-rt7dsq6d
authors: Rossa, Fabio Della; Salzano, Davide; Meglio, Anna Di; Lellis, Francesco De; Coraggio, Marco; Calabrese, Carmela; Guarino, Agostino; Cardona, Ricardo; DeLellis, Pietro; Liuzza, Davide; Iudice, Francesco Lo; Russo, Giovanni; Bernardo, Mario di
title: Intermittent yet coordinated regional strategies can alleviate the COVID-19 epidemic: a network model of the Italian case
date: 2020-05-15
journal: nan
DOI: nan
sha: 34d802f246e31c660bf5a442ddbf2a8232636d96
doc_id: 466260
cord_uid: rt7dsq6d

The COVID-19 epidemic that emerged in Wuhan China at the end of 2019 hit Italy particularly hard, yielding the implementation of strict national lockdown rules (Phase 1). There is now a hot ongoing debate in Italy and abroad on what the best strategy is to restart a country to exit a national lockdown (Phase 2). Previous studies have focused on modelling possible restarting scenarios at the national level, overlooking the fact that Italy, as other nations around the world, is divided in administrative regions who can independently oversee their own share of the Italian National Health Service. In this study, we show that regionalism, and heterogeneity between regions, is essential to understand the spread of the epidemic and, more importantly, to design effective post Lock-Down strategies to control the disease. To achieve this, we model Italy as a network of regions and parameterize the model of each region on real data spanning almost two months from the initial outbreak. Using the model, we confirm the effectiveness at the regional level of the national lockdown strategy implemented so far by the Italian government to mitigate the spread of the disease and show its efficacy at the regional level. We also propose that differentiated, albeit coordinated, regional interventions can be effective in Phase 2 to restart the country and avoid future recurrence of the epidemic, while avoiding saturation of the regional health systems and mitigating impact on costs. Our study and methodology can be easily extended to other levels of granularity (provinces or counties in the same region or states in other federal countries, etc.) to support policy- and decision-makers.

Regionalism is an integral part of the Italian constitution. Each of Italy's twenty administrative regions is independent on Health and oversees its own share of the Italian National Health service. The regional presidents and their councils can independently take their own actions, strengthening or, at times, weakening national containment rules. Previous studies have modelled the spread of the epidemics and its evolution in the country at the national level [1] [2] [3] [4] [5] , and some have looked at the effects of different types of containment and mitigation strategies [6] [7] [8] [9] [10] [11] . Limited work [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] has taken into account the spatial dynamics of the epidemic but, to the best of our knowledge, no previous paper in the literature has explicitly taken into consideration the pseudo-federalist nature of the Italian Republic and its strong regional heterogeneity when it comes to health matters, hospital capacity, economic costs of a lockdown and the presence of inter-regional people's flows.

In this study, we investigate the whole of the country as a network of regions, each modelled with different parameters. The goal is to identify if and when measures taken by the Italian government had an effect at both the national, but most importantly, at the regional level.

Also, we want to uncover the effects on the epidemic spread of regional heterogeneity and inter-regional flows of people and use control theoretic tools to propose and assess differentiated interventions at the regional level to reopen the country and avoid future recurrent epidemic outbreaks.

As aggregate models of the COVID-19 epidemic cannot capture these effects, to carry out our study we derived and parameterized from real data a network model of the epidemics in the country (see Figure 1a) , where each of the 20 regions is a node and links model both proximity flows and long-distance transportation routes (ferries, train, planes). The model is first shown to possess the right level of granularity and complexity to capture the crucial elements needed to correctly predict and reproduce the available data. Then, it is used to design and test differentiated feedback interventions at the regional level to alleviate the epidemic impact.

Using the model and an ad hoc algorithm to parameterize it from real data, we evaluate the effectiveness of the national lockdown strategy implemented so far by the Italian government providing evidence of its efficacy across regions. Also, we show that inter-regional fluxes must be carefully controlled as they can have dramatic effects on recurrent epidemics waves. Finally, we convincingly show that regional feedback interventions, where each of the twenty regions strengthens or weakens local mitigating actions (social distancing, inflow/outflow control) as a function of the saturation of their hospital capacity, can be beneficial in mitigating possible outbreaks and in avoiding recurrent epidemic waves while reducing the costs of a nation-wide lockdown.

Our approach successfully uncovers the regional effects of the national lock-down measures set in place by the Italian government initially in two northern regions (Lombardy and Veneto from the February 27 th 2020), and then nationally from March 8 th 2020 till May 4 th 2020. We observe that notable parameter changes, detected automatically by our parameterization procedure (see Methods), occur as an effect of such measures with a certain degree of homogeneity across all regions (see Figures S8, S9 and Table S4 in the Supplementary   Information showing the changes in the social distancing parameter over the period of interest). This confirms the effectiveness across the country of the strict social distancing rules implemented at the national level as also noted in previous work 1, 2, 14 modeling the country as a whole.

The representative examples of two regions, Lombardy in the North and Campania in the South, highlighted in Figure 1 , show that the model is able to correctly capture the effect of such measures in both the regions, see Table 1 . The model also captures the effect of the flow of people that travelled from North to South when the national lockdown measures where first announced on March 8 th 2020. As shown in Table 1 , the estimated number of infected predicted by the model for the Campania region in the time window March 19 th -March 30 th 2020 is detected to suddenly increase at the beginning of the next time window. This can be explained as a possible effect of the movement of people from North to South that occurred around 15 days before. Also, the data analysis consistently shows across many regions that the mortality rate changes dramatically as a function of the saturation of the hospital beds in the region (see Figure S11 and Section S2 of the Supplementary Information for further details).

After confirming the predictive and descriptive ability of the proposed model, we investigated next the influence of regional heterogeneity on the onset of an epidemic outbreak and the occurrence of possible recurrent epidemic waves. To this aim we set the model with parameters capturing the situation in each region on May 3 rd 2020, when the effects of the national lock down were fully in place, and simulated the scenario where just one of the twenty regions (e.g., Lombardy in the North of Italy) fully relaxes its lockdown. As reported in Figure   2 , we found that a primary outbreak in that region would quickly propagate causing secondary recurrent outbreaks in other regions including Emilia-Romagna and Piedmont. At the national level this would cause the onset of a second epidemic wave that, if not contained, would end up afflicting more than 25% of the entire population. An even more dramatic scenario would emerge if inter-regional flows where concurrently restored to their pre-lockdown levels (see Figure S1 ) or all regions were to relax their current restrictions concurrently (see Figure S2 in the Supplementary information).

A crucial open problem is to support decision makers in determining what form of interventions might be beneficial to avoid the onset of future outbreaks while mitigating the cost of Draconian interventions at the national level. To this aim we compared the effects of national measures (e.g., general lockdown) against those of a regional feedback strategy where social distancing measures are put in place or relaxed independently by each region according to the ratio between hospitalized individuals and the total capacity of the health system in that region. In particular, we assume each region implements a stricter lockdown when such a ratio becomes greater or equal than 20% and relaxes the social distancing rules when it is below 5% (see Methods for further details). Figure 3 confirms the effectiveness of such a local strategy, where we see that a differentiated strategy among the regions is more effective than a national lockdown in avoiding future waves of the disease ( Figure S3 ) but, most importantly, also in guaranteeing that no region exceeds its own hospitals' capacity. Moreover, intermittent regional measures yield lower economic costs for the country, as regional economies can be restarted and remain open for a much longer time (Table 2 and Table S1 ); an effect that becomes even more apparent when regions concurrently increase their testing capacity as shown in Figure 4 and reported in Table 2 .

Following the initial COVID-19 outbreak in Northern Italy, the Italian government, as many other governments around the world, adopted increasingly stricter lockdown measures at the national level to mitigate the epidemic. Despite their success, their tolling economic costs have stirred a hot national debate on whether such measures were necessary in the first place and on how to relax them while avoiding future epidemic waves. Several attempts have been made in the literature at addressing this pressing open issues by means of aggregate models (originated from the classical SIR model) to describe the effects of different intervention strategies at the national level 6, 7 . A network model has also been recently proposed to describe the spatial dynamics of the spread of the COVID-19 epidemics among the 107 Italian provinces 14 . Other works in the literature have explored the effects of intermittent measures, either periodic or as a function of some observable quantities, as a viable alternative to long, continuous periods of national lockdown. However, the effects of these strategies have only been investigated on theoretical aggregate models at the national level 7, 22 .

An important missing aspect that we considered in our study is the effect of regional heterogeneity on the efficacy of the measures taken so far and the possibility of adopting differentiated and localized intervention strategies thanks to the pseudo-federalist administrative structure of the Italian Republic. Our results confirm the effectiveness at the regional level of the national lockdown measures taken so far. They also convincingly reveal the presence of important regional effects due, for example, to the saturation of regional healthcare systems or to the presence of notable North-South flows in the country that followed the announcement of national measures. Also, contrary to previous work, we explicitly accounted for the strongly nonlinear nature of the model and the uncertainty present in the data by performing a sensitivity analysis on the estimated parameters that further confirmed the robustness of the proposed strategies for a wide range of parameter changes Our study strongly suggests for policy-and decision-makers the potential benefits of differentiated (but coordinated) feedback regional interventions, which can be used independently or in combination with other measures, in order to avoid future epidemic waves or even to contain the outbreak of potential future epidemics. Despite having been focused on Italy, our methodology and modeling approach can be easily extended to other levels of granularity, e.g. countries in a continent or counties in a state, and adapted to any other nation where regional heterogeneity is important and cannot be neglected; notable examples are countries with a federal state organization such as Germany or the United States of America.

Future work needs to address further aspects as for example exploring how the structural properties of the inter-regional network can influence the dynamics of the epidemic or adopting more sophisticated cost functions to design more effective region-specific mitigation strategies in other contexts or for other purposes.

Regional and national model As a regional model of the COVID-19 epidemic spread we use the compartmental model shown in Figure 4 which we found from data analysis and identification trials to be the simplest model structure able to capture the real data. The full model equations describe the dynamics of susceptible ( ), infected ( ), quarantined ( ), hospitalized ( ), recovered ( ) and deceased ( ) and are given as:

where is the infection rate which will be assumed to be the same for all regions since COVID-19 is transmitted from person to person and there is no parasite vector or evidence of environmental parameters significantly altering its infection rate, ∈ [0,1] is a parameter modeling the effects of social distancing measures in the i -th region , is the rate of infected that are detected and quarantined, is the rate of infected that needs to be hospitalized, is the recovery rate of the infected that is assumed to be equal for all regions, is the rate of quarantined who recover, is the fraction of hospitalized who recover, is the rate of hospitalized that is transferred to home isolation, is the rate of quarantined who need to be hospitalized, and is the mortality rate that was shown from data analysis (see Section S3

of the Supplementary Information) to be a function of the ratio between and the maximum number, say , of patients that can be treated in ICU at the hospitals in i-th region . is the actual population in the i-th region, i.e. the resident population without those removed because quarantined, hospitalized, deceased or recovered.

Extending previous approaches for modelling Dengue fever in Brazil 23 , we obtain the national network model of the COVID-19 epidemic in Italy as a network of twenty regions (see Figure   1a ) interconnected by links modeling commuter flows and major transportation routes among them.

The network model of Italy we adopt in this study is, for = 1, … ,20,

where in addition to the parameters and states described above, we included the fluxes ( ) between regions;

( ): ℝ → [0,1] denoting the ratio of people from region i interacting with those in region j at time t, such that ∑ ( ) = 1.

We divide the model parameterization into two stages. Firstly, we estimate from the available data the parameters of each of the twenty regional models; then, we use publicly available mobility data in Italy to estimate the fluxes among the regions.

We follow the current consensus in the literature on COVID-19 2 (see Table S5 ) by setting = 0.4 and = 1/14 for all regions. We make the ansatz that parameters remain constant over time intervals but do not assume the duration of such intervals known a priori. Therefore, we set the problem of estimating the parameters values and when they change in each region (as a likely result of national containment measures). We start estimating parameters in a region from the first date when the number of deceased and the number of recovered is greater or equal than 10.

Noting that the available data for the COVID-19 epidemics, as collected by the Dipartimento della Protezione Civile -Presidenza del Consiglio dei Ministri (the Italian Civil Protection Agency), includes for each region the daily numbers of quarantined ( � ), hospitalized ( � ), deceased ( � ) and the recovered from those who were previously hospitalized or quarantined, say � , we discretize and rewrite model (1)-(6) for each region (i = 1,…,20) as (dropping the subscripts to the parameters for notational convenience)

where measured quantities are denoted by a tilde and estimated state variables by a hat and : = + . Here, = + + + represents the total number of cases detected in region i as daily announced by the Protezione Civile.

We notice that, exploiting the available data, the predictor can be split into two parts so that two different algorithms can then be used to estimate the parameters of each part. An ad hoc identification algorithm estimates the parameters of the nonlinear part (Equations (14)-(16)) and automatically detects the breakpoints where notable parameter changes occur, while an ordinary least squares method is then used to identify the parameters of for the linear one (Equations (17)-(20)), as described in detail in Section S2 of the Supplementary Information.

Note that, as the actual number of infected is not known 9,24 , we include the number of infected at the beginning of each time window as a parameter to be estimated by the algorithm used for the nonlinear part.

The results of the identification process also show the presence of a statistically significant correlation (p-value equal to 0.034) between the value of the mortality rate, parameter in model (1)- (6) , and the saturation of the regional health system represented by the ratio between the number of hospitalized in that region ( ) and the total number of available hospital beds ( ) (See Supplementary Information Section S2 and Figure S11 for further details and function estimation).

Validation is carried out by using the parameterized model to capture the available data for each window showing an average quadratic error less than 10% over the entire dataset. The parameters identified in each window can also be used to provide model predictions of future trends of the epidemic diseases as discussed in Section S2 of the Supplementary Information.

We estimate the cost of each regional lockdown as the sum of the costs for social care and the loss of added value. The costs for social care in each region were computed as the costs for layoff support (``cassa integrazione in deroga''), estimated by multiplying the number of requests 26 by 65% of the average regional monthly income 27 , together with the non-repayableloan of 600 € given to self-employed workers by the Italian Government during the national lockdown 28 . The loss of added value per day was taken from the values estimated by SVIMEZ (the Italian Association for the development of Industry in the South) in Table 3 of their online report 29 . We then compute the daily costs of the lockdown and use it to estimate the total costs of each of the simulated scenarios.

All computational analyses and the fitting of data were performed using MATLAB and its optimization toolbox. To account for the inherent uncertainty associated to the COVID-19 epidemic, each result reported in the manuscript is the output of 10,000 numerical simulations,

where we varied the values of the model parameters using the Latin Hypercube sampling method 25 . Specifically, the regional parameters , , , , , and the estimated initial condition at May 3rd 2020 , were varied considering a maximum variation of ±20% from their nominal values (indicated in Table S4 of the Supplementary Information).

Implementation and design of national and regional feedback intervention strategies

We model the implementation of regional social distancing strategies by capturing their effects as a variation of the social distancing parameters, in (7)-(8), in each region.

Specifically, we assume each region follows the feedback control rule:

where is set equal to the minimum estimated value in that region during the national lockdown (see Table S9 ) and ̅ i increased as a worst case to min(1,3 ) so as to simulate the effect of relaxing the lock-down measures in each region. (The case where ̅ i is set to a lower value equal to 1.5 is shown for the sake of comparison in Figure S5 and Figure To model the increase in the COVID-19 testing capacity of each region the parameter in region i is multiplied by a factor 2.5 which corresponds to the average increase in the number of tests carried out nationally since the COVID-19 outbreak first started.

All simulations were carried out in MATLAB with a discretization step of 1 day to match the available data sampling interval. Initial conditions for regional compartments were set as follows. Quarantined, Hospitalized, Deceased and Recovered are initially set to the datapoints available for May 3rd 2020. The number of infected is set to the value estimated by our procedure for that date while Susceptibles are initialized to the resident population form which the other compartments are removed.

The code is available at https://github.com/diBernardoGroup/Network-model-of-the-COVID-19

Further details are given in Section S5 of the Supplementary Information.

The authors declare that the data supporting the findings of this study are available within the paper and its Supplementary 

The authors declare no competing interests. shown in Figure S1 of the Supplementary Information showing even more dramatic scenarios). 15 Panels of regions adopting a lockdown are shaded in red while those of regions relaxing social 16 containment measures are shaded in greed. Blue, magenta, red, green, and black solid lines 17 correspond to the fraction in the population of infected, quarantined, hospitalized, recovered, 18 and deceased over the entire regional (panels a) or national (panel b) population averaged year of each of the simulated scenarios are reported showing the effectiveness of the 78 intermittent regional measures shown in Figure 3 and Figure Table S4 . Intermittent regional measures (Fig. 3a,b Interregional Fluxes estimation We considered two types of inter-regional fluxes associated to (i) daily commuters traveling between neighboring regions;

(ii) long distance travels covered by high-speed trains, planes, and large ferries.

To estimate commuters' fluxes, as in previous work S6 , we use the latest official country-wide assessment of Italian mobility conducted by the Italian Institute of Statistic (ISTAT) in 2011.

Specifically, we use the origin-destination matrix available at https://www.istat.it/it/archivio/139381 describing, at the municipality level, the number of people who declared themselves daily commuters for work or studying purposes. By aggregating data on a regional basis, we obtained interregional commuters' flows among regions.

For longer distance routes covered by high-speed trains, planes and large ferries we proceeded as follows.

a. For railway connections, the fraction ! !" # of the population of region i traveling to region " is obtained as

• ! # is total number of daily high-speed customers of the main Italian carrier, Trenitalia;

• & # is the total number of high-speed trains per day;

• & !" # is the number of high-speed trains connecting region ( to region ";

• % ! is the total resident population in region i.

b. For flight connections, for all the Italian operating airports (currently there are 39 according to www.enac.gov.it) and for each pair of connected airports, say (*, +) the number of weekly direct flights -$% in normal operating conditions (i.e., without any restriction due to the emergency). Then, for each pair (*, +) we computed the average capacity ⟨/ $% ⟩ of the regional fleet of the main carrier serving the route. Grouping the airports on a regional basis we compute⟨/ $% ⟩ average daily flow due to air connections ! !" & , from region ( to " as:

where 3 ! and % ! are the number of airports and the population of region (. c. For ferry connections, we considered the five regions that act as hubs for long range national maritime travel, that is, the two main insular regions, Sardinia and Sicily, together with three mainland regions, Campania, Lazio, and Liguria. The maritime flows are then obtained as

where ! ! and % ! are the average number of maritime passengers and the population of region ( respectively, while & !" * is the total number of maritime connections between regions ( and ". Note that as it is reasonable to assume that & ! ! !" ≈ & " ! "! , and as all relevant maritime routes are either from or to the main islands, in practice, it suffices to compute ! !" * only for region ( being Sicily and Sardinia.

As explained in the main text, we assume that, for each region, the parameters of the model remain constant over & time windows, but neither their number & nor their durations 6 ) , … , 6 + are assumed to be known a priori. Therefore, the identification procedure detects at the same time the breakpoints 8 ) , … , 8 +,) when notable parameters' changes are detected and, within each time-window, estimates their values that best capture the trend of the available data.

The model used to carry out the identification of regional or national parameters is to start with the discretized version of the model of the epidemic spread in each area of interest given by model (14)- (20) in the main text.

As also noted in other previous work S1-S4 , identification of SIR and SIR-modified models is highly non-convex and hence the optimization landscape is scattered with local minima that must be avoided as not being admissible. To mitigate this problem, we identified from the literature admissible intervals for the parameter values (see Table S5 ) and fixed from the available evidence in the literature S3 the parameters 9 = 0.4 and = = 1/14.

Altogether, the unknown parameters left to be estimated are [@(0), A, B, C, D -, D . , E . , E -, F] both at the national and regional level.

As mentioned in the main text, the identification procedure is carried out in two stages by considering equations (14) We use the following recursive procedure:

1. Set the initial time 8 4 as the first day in which the first 10 deaths and 10 recovered were reported in the area (region or nation). the data over the window (8 4 , 8 ) ) whose duration is therefore 6 ) : = 8 ) − 8 4 .

7. If 8 ) = N /0/ the algorithm is stopped, otherwise starting from 8 ) steps 2-7 are repeated to find the next breakpoint and the new set of parameters best estimating the data in the next window until the end of the available datapoints.

Step 2: Offline refinement of the identification process

At the end of the process we will have the set of breakpoints 8 " and the parameters set in each of the windows (8 " , 8 "1) ) best fitting the data. As the number of windows can be large given

The variability in the available data, we refine the estimation results as follows to estimate the minimal number of windows able to capture qualitatively the trend of the real data.

In particular, once the window breakpoints are obtained at the end of step 1, any two consecutive windows of duration say 6 " , 6 "1) are merged into one larger window of size 6 " + 6 "1) if one of the two following conditions is verified:

a) The window size of the first window is less than 5 days.

b) The relative variation of the sum B + A as estimated in each of the two windows is less than 5%, i.e.

where j denotes the window to which the parameter estimates refers to.

If two windows are merged, then the parameters are estimated again on the entire merged window and the procedure is iterated once more in case condition b) is still verified.

As a final refinement step, we heuristically explore the effect on the fitting of perturbing the breakpoints within ±5 days from their estimated value. As a representative example, the results of the fitting procedure for the national aggregate model are shown in Table S3 and depicted in Figure S9 . The same procedure is repeated to parametrize each of the 20 regional models.

To provide a representative validation of our estimation approach, we report in Figure S10 the time evolution of the total number of detected cases at the national level predicted by the model in each time window. It is possible to see that using just 30% of the datapoints (shown in green) from all the available data (shown in red), the model predictions (solid blue lines) fit well the rest of the data in each time window both before and after the windows are merged as a result of step 2 with a maximum prediction error of 10000 units.

Step Table S2 to compute the remaining parameters.

The comparison between the model predictions and the available data is depicted in Figure S11 and Figure S12 . Values of all estimated parameters at the end of the process are given in windows that are automatically detected by the estimation procedure we proposed. for each region. In the last column the regional net reproduction numbers are computed as p 4,! = A ! 9/(C ! + o ! + =) in each time window. Figure S13 shows the distribution of the regional social distancing parameters over time showing the effects of the national lockdown at the regional level.

Observing the data and the identified parameters in Table S4 , we noticed a significant correlation between the mortality rate in each time window and the congestion of the ICU system in that region.

Specifically, we found that

where F 4 ?@ and F ) ?@ are coefficients to be estimated, while g q ! is the estimated average congestion of the hospitals in each time-window defined as the ratio between the number of hospitalized people and the number of available beds in ICU in that region, say N ! . (available from the web-page of the Italian Ministry of Health S8 and updated daily from the news published by the main Italian newspapers in each region).

Each point in Figure S14 is the value F ! estimated for each region in each time window against the number of hospitalized in the corresponding region averaged in the time window. As illustrated in Figure S14 , a least square linear fitting yields F 4 ?@ = 0.016, F ) ?@ = 0.0013.

We considered three different types of intervention strategies at the regional level that can be used individually or in combination.

Each region modulates its lockdown measures so as to switch them on or off according to the relative saturation level of its health system.

Namely, the social distancing parameter is modulated as follows:

where A ! is set equal to the minimum estimated value in that region during the national lockdown (corresponding to the value given in the last window for each region reported in Tab. S4 ) and, as a worst case scenario, A̅ ! is set to three times A ! (or unity if 3A ! >1) to simulate a relaxation of the containment measures in that region. Note that all people resident in ( not commuting to any other region are assumed to stay and move in region ( itself by setting

All the numerical analyses presented in the paper were performed with Matlab. The code is available at https://github.com/diBernardoGroup/Network-model-of-the-COVID- 19 and was designed to 1) Load the estimated parameter values and the inter-regional fluxes and to iterate the discretized model dynamics of each region considering the presence of interregional fluxes (MATLAB script "siqhrd_network_main.m"). The script also computes the regional and national net reproduction numbers.

2) Implement the differentiated regional intervention strategies described in Section S4. Further details can be found in the accompanying README.TXT file included with the code in the software repository above. Figure S1 [Phase 2, only one region relaxes its lockdown] a. Regional and b. national dynamics in the case where only one region (Lombardy in Northern Italy) relaxes its containment measures at time 0 and the fluxes between regions are set to their pre-lockdown level. Panels of regions adopting a lockdown are shaded in red while those of regions relaxing social containment measures are shaded in green. blue, magenta, red, green, and black solid lines correspond to the fraction in the population of infected, quarantined, hospitalized, recovered, and deceased averaged over 10,000 simulations with parameters sampled using a Latin Hypercube technique around their nominal values set as those estimated in the last time window for each region as reported in Table S4 . Shaded bands correspond to twice the standard deviation. The black dashed line identifies the total fraction of the population that can be treated in ICU (N ! . /3 ! ). The regions identified with a red label are those where the total hospital capacity is saturated. a b Figure S2 [Phase 2, all regions relax their lockdown measures] ] a. Regional and b. national dynamics in the case where all regions relax their current restrictions restoring fluxes to their pre-lockdown level. Blue, magenta, red, green, and black solid lines correspond to the fraction in the population of infected, quarantined, hospitalized, recovered, and deceased in the population averaged over 10,000 simulations with parameters sampled using a Latin Hypercube technique+ around their nominal values set as those estimated in the last time window for each region as reported in Table S4 . Shaded bands correspond to twice the standard deviation. The black dashed line identifies the fraction of the population that can be treated in ICU (N ! . /3 ! ). The regions identified with a red label are those where the total hospital capacity is saturated. a b Figure S3 [National lockdown] ] a. Regional and b. national dynamics in the case where no region relaxes its containment measures, while all regions restore the interregional fluxes to their pre-lockdown level. Blue, magenta, red, green, and black solid lines correspond to the fraction in the population of infected, quarantined, hospitalized, recovered, and deceased averaged over 10,000 simulations with parameters sampled using a Latin Hypercube technique around their nominal values set as those estimated in the last time window for each region as reported in Table S4 . Shaded bands correspond to twice the standard deviation. The black dashed line identifies the fraction of the population that can be treated in ICU (N ! . /3 ! ). The regions identified with a red label are those where the total hospital capacity is saturated.

a Figure S4 [Phase 2, regional dynamics when intermittent national measures are enforced as shown in Fig. S3c of the main text] . Each of the 20 panels shows the evolution in the corresponding region of the fraction in the population of infected (blue), quarantined (magenta) and hospitalized requiring ICUs (red) averaged over 10,000 simulations with parameters sampled using a Latin Hypercube technique around their nominal values set as those estimated in the last time window for each region as reported in Table S9 . Shaded bands correspond to twice the standard deviation. The black dashed line identifies the fraction of the population that can be treated in ICU (N ! . /3 ! ). National lockdown measures are enforced with all regions shutting down when the total number of occupied ICU beds at the national level exceed 20% (windows shaded in red, green when relaxed). The regions identified with a red label are those where the total hospital capacity is saturated. a b Figure S5 [Phase 2, intermittent regional measures, q D = Ä. Å D ]. a. Each of the 20 panels shows the evolution in the region named above the panel of the fraction in the population in the population of infected (blue), quarantined (magenta) and hospitalized requiring ICUs (red) averaged over 10,000 simulations with parameters sampled using a Latin Hypercube technique around their nominal values set as those estimated in the last time window for each region as reported in Table S4 of the Supplementary Information. Shaded bands correspond to twice the standard deviation. The black dashed line identifies the fraction of the population that can be treated in ICU (N ! . /3 ! ). Regions adopt lockdown measures in the time windows shaded in red while relax them in those shaded in green. Differently from Figure 3 , when lockdown measures are relaxed A̅ ! is set to 1.5 times A ! . During a regional lockdown, fluxes in/out of the region are set to their minimum level. b. National evolution of the fraction in the population of infected (blue), quarantined (magenta) and hospitalized requiring ICUs (red) obtained by summing those in each of the 20 regions adopting intermittent regional measures. Table S4 . Shaded bands correspond to twice the standard deviation. The black dashed line identifies the fraction of the population that can be treated in ICU (N ! . /3 ! ). National lockdown measures are enforced with all regions shutting down when the total number of occupied ICU beds at the national level exceed 20%. The regions identified with a red label are those where the total hospital capacity is saturated. Regions adopt lockdown measures in the time windows shaded in red while relax them in those shaded in green. Differently from Figure  S4 , when lockdown measures are relaxed A̅ ! is set to 1.5 times A ! . During a regional lockdown, fluxes in/out of the region are set to their minimum level. b. National evolution of the fraction in the population of infected (blue), quarantined (magenta) and hospitalized requiring ICUs Table S4 . Shaded bands correspond to twice the standard deviation. The black dashed line identifies the fraction of the population that can be treated in ICU (N ! . /3 ! ). Regions adopt lockdown measures in the time windows shaded in red while relax them in those shaded in green. During a regional lockdown, fluxes in/out of the region are set to their minimum level. Regions COVID-19 testing capacities are assumed to be increased by a factor 2.5 (see Methods) with respect to their current values. Differently from Figure 4 , when lockdown measures are relaxed A̅ ! is set to 1.5 times A ! . During a regional lockdown, fluxes in/out of the region are set to their minimum level. b. National evolution of the fraction in the population of infected (blue), quarantined (magenta) and hospitalized requiring ICUs (red) obtained by summing those in each of the 20 regions. Figure S11 [Identification of the regional models -Steps 1,2] Comparison of each of the regional model predictions for the total number (expressed in thousands of people) of detected cases in each region (solid magenta line) against the available data points. Parameters are set to the values estimated at the end of Steps 2 carried out for each region. Vertical black lines denote the breakpoints from one estimation window to the next. Figure S12 [Identification of the regional models -Step 3] Comparison of each of the regional model predictions of the total number (expressed in thousands of people) of recovered (green), quarantined (magenta), hospitalized (red), deceased (black) and recovered (green) in each region against the available data points (plotted as circles of the same colour). Parameters are set to the values estimated at the end of Step 3 carried out for each region. Vertical black lines denote the breakpoints from one estimation window to the next. Table S4 Figure S14. F ! fitting of g q ! as a function of p 2 = 0.280 % < 0.001 and % < 0.001. Normality of the residuals has been tested with the Lilliefors test (p-value 0.12). Intermittent national measure, Figure 3c , Figure S4 2 

Continuity constraint on the number of infected at the national level from window j to window j+1.

Continuity constraint on the number of infected at the regional level. The constraints are relaxed by 10% of the national estimated infected to account for the fact that in estimating the region parameter we are neglecting the influx of infected from other regions. Ê-≤ 0.1, Ê. ≤ 0.1;

We assume that the daily number of people hospitalized from quarantine and discharged but still positive (and vice versa) is no higher than 10% of the total.

We assume B! = C â ! + o O ! does not differ from the national estimate B̂ of more than 30%. D -,"1) = D -,"

We assume the recovery rate of those quarantined at home remains the same from a time window to the next as this parameter is likely to be time-invariant. In any case, removing this constraint, we observed no significant change of this parameter from a time window to the next. Step 1

Step 2 If  ηQ  ηH  ζ  α  ψ  κH  κQ  ti  R0  Abruzzo 0,485 944  1083 0,010 0,003 0,026 0,029 0,051 0,005 0,099 21-mar-20 1,29  0,321 1083 747  0,010 0,000 0,022 0,025 0,049 0,000 0,087 25-mar-20 0,89  0,194 847  210  0,010 0,019 0,014 0,078 0,003 0,005 0,000 14-apr-20  0,51  Aosta  0,283 425  262  0,010 0,096 0,036 0,062 0,016 0,028 0,000 30-mar-20 0,77  0,122 262  41  0,010 0,260 0,011 0,062 0,000 0,079 0,000 14-apr- Values of estimated parameters for each region at the end of the identification process. Dates are given corresponding to breakpoints between estimation windows that are automatically detected by the estimation procedure we proposed. Regional net reproduction number are reported in the last column clearly showing the increasing effect of the national lockdown measures taken by the government on March 8 th 2020. The parameter values used to simulate the exit from the Lock-Down considered in this paper are those shaded in grey for each region corresponding to the estimated parameter in the last window ending on May 3 rd 2020. Table S5 [Comparison with the parameters estimated in other papers in the literature for the national aggregate models]. Comparison between the parameter's values we used in our work and those used in other papers proposing national models for the COVID-19 epidemic in Italy that recently appeared in the literature. Notice that, unfortunately, often it is not possible to pin down a specific parameter in a different model that clearly corresponds to one of ours, and vice-versa. This is because of the different meaning that compartments have in the models and the dynamics of people between the compartments which do not always overlap in an unambiguous manner. When we had to use a time constant, say τ, to determine the value of a parameter, say k, that describes the rate of change form compartment A to B, we proceeded as follows. We set k = p / τ, where p ∈ [0, 1] is the percentage of people that has been estimated to actually move from compartment A to B .   nom  min  max  nom min  max  nom min  max  nom min  max  nom min  max  nom min  max  nom min  max  nom min  max  nom min  max  nom min  max  nom min  max  nom min  max   our model  SIQHRD  variable 0,122 1 -

SEII [S17]

Source [S6] [S3]

[S10]

[S11]

[S12]

[S13]

[S14]

[S15]

[S2]

--β pre-lockdown β post-lockdown ---

A Modified SIR Model for the COVID-19

Contagion in Italy

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

The Mathematics of Infectious Diseases

Parameter Identification for a Stochastic SEIRS Epidemic Model: Case Study Influenza

Epidemic Analysis of COVID-19 in China by Dynamical Modeling

Can the COVID-19 Epidemic Be Controlled on the Basis of Daily Test Reports?

On Fast Multi-Shot Epidemic Interventions for Post Lock-Down Mitigation: Implications for Simple Covid-19 Models

Impact of Non-Pharmaceutical Interventions (NPIs) to Reduce COVID-19 Mortality and Healthcare Demand

Estimating the Number of Infections and the Impact of Non-Pharmaceutical Interventions on COVID-19 in 11

Optimal Resource Allocation for Control of Networked Epidemic Models

A Distributed Dynamics for Virus-Spread Control

The Effect of Network Topology on the Spread of Epidemics

Spread and Dynamics of the COVID-19 Epidemic in Italy: Effects of Emergency Containment Measures

An Immunization Model for a Heterogeneous Population

A Deterministic Model for Gonorrhea in a Nonhomogeneous Population

On the Dynamics of Deterministic Epidemic Propagation over Networks

A Structured Epidemic Model Incorporating Geographic Mobility among Regions

On the Predictability of Infectious Disease Outbreaks

Virus Spread in Networks

An Epidemic Model in a Patchy Environment

An Alternating Lock-down Strategy for Sustainable Mitigation of COVID-19

SIR-Network Model and its Application to Dengue Fever

Substantial Undocumented Infection Facilitates the Rapid Dissemination of Novel Coronavirus (SARS-CoV2)

Latin Hypercube Sampling and the Propagation of Uncertainty in Analyses of Complex Systems

Dati al 10 maggio su Cassa integrazione ordinaria, assegno ordinario, richieste di pagamento SR41 e Cassa integrazione in deroga

Stipendi, l'Italia è uno 'scivolo': Regioni e province dove si guadagna più

Osservatori statistici e altre statistiche

SVIMEZ Associazione per lo sviluppo dell'industria nel mezzogiorno, L'impatto economico e sociale del Covid-19: Mezzogiorno e Centro-Nord

Parameter Identification for a Stochastic SEIRS Epidemic Model: Case Study Influenza

The Reproductive Number of COVID-19 Is Higher Compared to SARS Coronavirus

Modelling the COVID-19 Epidemic and Implementation of Populationwide Interventions in Italy

A Modified SIR Model for the COVID-19 Contagion in Italy

On the Definition and the Computation of the Basic Reproduction Ratio R 0 in Models for Infectious Diseases in Heterogeneous Populations

Spread and dynamics of the COVID-19 Epidemic in Italy: Effects of Emergency Containment Measures

Epidemic Calculator

Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand

Epidemic analysis of COVID-19 in China by dynamical modeling

Can the COVID-19 epidemic be controlled on the basis of daily test reports

On Fast Multi-Shot Epidemic Interventions for Post Lock-Down Mitigation

Estimating the number of infections and the impact of nonpharmaceutical interventions on COVID-19 in 11 European countries

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV2)

Estimates of the severity of coronavirus disease 2019: a model-based analysis" The Lancet Infectious Diseases