# A network model of Italy shows that intermittent regional strategies can alleviate the COVID-19 epidemic

Fabio Della Rossa, Davide Salzano, Anna Di Meglio, Francesco de Lellis, Marco Coraggio, Carmela Calabrese, Agostino Guarino, Ricardo Cardona-Rivera, Pietro de Lellis, Davide Liuzza, et al.

► To cite this version:

Fabio Della Rossa, Davide Salzano, Anna Di Meglio, Francesco de Lellis, Marco Coraggio, et al. A network model of Italy shows that intermittent regional strategies can alleviate the COVID-19 epidemic. Nature Communications, 2023, 11 (1), pp.5106. 10.1038/s41467-020-18827-5. hal-03598263

HAL Id: hal-03598263

https://hal.science/hal-03598263

Submitted on 11 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

ARTICLE

https://doi.org/10.1038/s41467-020-18827-5

OPEN

<figure><img src="image_3.png" /></figure>

Check for updates

# A network model of Italy shows that intermittent regional strategies can alleviate the COVID-19 epidemic

Fabio Della Rossa ID 1,2,5, Davide Salzano ID 2,5, Anna Di Meglio ID 2,5, Francesco De Lellis ID 2,5, Marco Coraggio ID 2, Carmela Calabrese ID 2, Agostino Guarino ID 2, Ricardo Cardona-Rivera ID 2, Pietro De Lellis ID 2, Davide Liuzza ID 3, Francesco Lo Iudice ID 2, Giovanni Russo ID 4 & Mario di Bernardo ID 2✉

The COVID-19 epidemic hit Italy particularly hard, yielding the implementation of strict national lockdown rules. Previous modelling studies at the national level overlooked the fact that Italy is divided into administrative regions which can independently oversee their own share of the Italian National Health Service. Here, we show that heterogeneity between regions is essential to understand the spread of the epidemic and to design effective strategies to control the disease. We model Italy as a network of regions and parameterize the model of each region on real data spanning over two months from the initial outbreak. We confirm the effectiveness at the regional level of the national lockdown strategy and propose coordinated regional interventions to prevent future national lockdowns, 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 to support policy- and decision-makers.

<footnote n="1">¹Department of Electronic, Information and Biomedical Engineering, Politecnico di Milano, Milan, Italy. ²Department of Electrical Engineering and Information Technology, University of Naples Federico II, Naples, Italy. ³ENEA, Fusion and Nuclear Safety Department, Frascati, Rome, Italy. ⁴Department of Information and Electrical Engineering and Applied Mathematics, University of Salerno, Fisciano, Italy. ⁵These authors contributed equally: Fabio Della Rossa, Davide Salzano, Anna Di Meglio, Francesco De Lellis. ✉email: mario.di.bernardo@unina.it</footnote>
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¹⁻⁵, and some have looked at the effects of different types of containment and mitigation strategies⁶⁻¹¹. Limited work¹²⁻²¹ 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 epidemic in the country (see Fig. 1a), where each of the 20 regions is a node and the links model both proximity flows and long-distance transportation routes (ferries, trains and 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 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 epidemic 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 nationwide lockdown.

## Results

**Model formulation and fitting procedure.** To capture the regional diversity of the response to the epidemic in Italy, we derive a network model where each node represents a different region and links capture fluxes of people traveling among the regions (see Fig. 1a). Using a data-driven compartmental modelling approach, a set of ODEs is obtained describing the dynamics of six different compartments in each region (Susceptibles, Infected, Hospitalized, Quarantined, Deceased and Recovered); data analysis being used (see Methods) to define flows among compartments. The resulting model is then parameterized using a predictor-corrector algorithm applied to both a national aggregate model and to each of the twenty regional models, identifying the time points at which parameter values present significant changes. Soft constraints are used to enforce

continuity of the trajectory between different time windows and avoid parameters changing too abruptly (see Methods and Supplementary Notes for further details). Estimating all the parameters in each region allows us to fit the available data and to describe the different regional situations and the diverse impact that regional policies had on the epidemic spread in each of the Italian regions.

As further explained in the Methods and Supplementary Information, we fit the model parameters to the official data for the COVID-19 epidemic²², as collected by the Dipartimento della Protezione Civile—Presidenza del Consiglio dei Ministri (the Italian Civil Protection Agency). Also, publicly available mobility data is used to estimate inter-regional fluxes and data on the number of ICU beds⁸,⁹ to evaluate the capacities of regional health services. To assess the economic costs of national and regional lockdowns we use official data and estimates from Italian governmental agencies²³⁻²⁶. Further details on the input data and the official repositories they were obtained from can be found at https://github.com/diBernardoGroup/Network-model-of-the-COVID-19.

**Regional effects of the national lockdown.** Our approach successfully uncovers the regional effects of the national lockdown measures set in place by the Italian government initially in two northern regions (Lombardy and Veneto from the February 27, 2020), and then nationally from March 8, 2020 till May 4, 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 Supplementary Fig. 13 and Supplementary Table 4 showing the changes in the social distancing parameter $\rho_i$ 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¹,²,¹⁴ modelling the country as a whole.

The representative examples of two regions, Lombardy in the North and Campania in the South, highlighted in Fig. 1, show that the model correctly captures 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 were first announced on March 8, 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–March 30, 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, data analysis shows that the mortality rate varies as a function of the level of occupancy of the hospital beds in each region (see Supplementary Fig. 14 and Supplementary Notes for further details).

**Regional heterogeneity counts.** After confirming the predictive and descriptive ability of the proposed model, we investigated the influence of the 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, 2020, when the effects of the national lockdown 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 Fig. 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
<figure><img src="image_4.png" /><figcaption>b</figcaption></figure>

<table><thead><tr><th>Index</th><th>Region</th><th>Index</th><th>Region</th></tr></thead><tbody><tr><td>1</td><td>Abruzzo</td><td>11</td><td>Lombardy</td></tr><tr><td>2</td><td>Aosta Valley</td><td>12</td><td>Marche</td></tr><tr><td>3</td><td>Apulia</td><td>13</td><td>Molise</td></tr><tr><td>4</td><td>Basilicata</td><td>14</td><td>Piedmont</td></tr><tr><td>5</td><td>Calabria</td><td>15</td><td>Sardinia</td></tr><tr><td>6</td><td>Campania</td><td>16</td><td>Sicily</td></tr><tr><td>7</td><td>Emilia-Romagna</td><td>17</td><td>Tuscany</td></tr><tr><td>8</td><td>Friuli-Venezia Giulia</td><td>18</td><td>Trentino-Alto Adige</td></tr><tr><td>9</td><td>Lazio</td><td>19</td><td>Umbria</td></tr><tr><td>10</td><td>Liguria</td><td>20</td><td>Veneto</td></tr></tbody></table>

**Fig. 1** Schematic diagram of the network-model structure and representative regional parameters. **a** Representative graph of the network-model structure used in the paper. Only a subset of all links is shown for the sake of clarity (the complete graphs are depicted in Supplementary Fig. 8). Solid lines represent proximity links, dashed lines long-distance transportation routes (planes and trains), dotted lines show major ferry routes between insular regions and the Italian mainland. **b** Table of the Italian region names and their labels in the graph.

<table><caption>Table 1 Estimated parameter values for Campania in the South (region no. 6) and Lombardy in the North (region no. 11), where the initial outbreak occurred.</caption><thead><tr><th>Region</th><th>Breakpoint</th><th>ρ<sub>i</sub></th><th>α<sub>i</sub></th><th>ψ<sub>i</sub></th><th>κ<sub>i</sub><sup>H</sup></th><th>κ<sub>i</sub><sup>Q</sup></th><th>η<sub>i</sub><sup>Q</sup></th><th>η<sub>i</sub><sup>H</sup></th><th>ζ<sub>i</sub></th><th>I<sub>0</sub></th><th>I<sub>f</sub></th><th>R<sub>0,i</sub></th></tr></thead><tbody><tr><td rowspan="2">Campania</td><td>19/3/20</td><td>0.467</td><td>0.014</td><td>0.064</td><td>0.000</td><td>0.100</td><td>0.018</td><td>0.000</td><td>0.022</td><td>1231</td><td>1816</td><td>1.26</td></tr><tr><td>30/3/20</td><td>0.221</td><td>0.067</td><td>0.019</td><td>0.006</td><td>0.040</td><td>0.018</td><td>0.000</td><td>0.011</td><td>2231</td><td>234</td><td>0.57</td></tr><tr><td rowspan="2">Lombardy</td><td>27/2/20</td><td>0.727</td><td>0.009</td><td>0.092</td><td>0.000</td><td>0.040</td><td>0.010</td><td>0.053</td><td>0.033</td><td>1799</td><td>28900</td><td>1.69</td></tr><tr><td>19/3/20</td><td>0.303</td><td>0.018</td><td>0.056</td><td>0.000</td><td>0.027</td><td>0.010</td><td>0.029</td><td>0.024</td><td>28900</td><td>6731</td><td>0.84</td></tr></tbody></table>

These regions are highlighted in a darker colour in Fig. 1. Here, $I_0$ is the number of infected estimated in the region at the beginning of each time window, while $I_f$ is the number of infected at the end of each time window estimated by running the model (14)–(16), given the set of identified parameters and the initial condition on the infected $I_0$. The first breakpoint is the date when 10 deaths and 10 recovered were first reported in the region and the analysis started. The second breakpoint is the end of the first window and the start of the second window (ending on May 3, 2020).

not contained, would end up afflicting more than 25% of the entire population. An even more dramatic scenario would emerge if inter-regional flows were concurrently restored to their pre-lockdown levels (see Supplementary Fig. 1) or all regions were to relax their current restrictions concurrently (see Supplementary Fig. 2).

**Feedback regional interventions can be beneficial.** 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 that each region implements a stricter lockdown when such a ratio becomes greater than or equal to 20% and relaxes the social distancing rules when it is below 10% (see Methods for further details). Figure 3 confirms that a differentiated strategy among the regions (Fig. 3b) is as effective as a national lockdown in avoiding future waves of the epidemic (Fig. 3c and Supplementary Fig. 4). At the same time, an intermittent regional strategy guarantees that no region exceeds its own hospitals’ capacity and yields a lower economic cost for the country (Table 2 and Supplementary Table 1), since regional economies can be restarted and remain open for a much longer time. This advantage becomes even more apparent when regions concurrently increase their testing capacity as shown in Fig. 4 and reported in Table 2.

## Discussion

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 high 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 these 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 epidemic 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,27}$.

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
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

**Fig. 2 Only one region relaxes its lockdown.** Double scale plots of the **a** regional and **b** national dynamics in the case where only one region (Lombardy in Northern Italy) relaxes its containment measures at time 0, while inter-regional fluxes are set to the level they reached during the lockdown. (Regional dynamics when the fluxes between regions are set to their prelockdown level are shown in Supplementary Fig. 1 showing even more dramatic scenarios.) The scale on the left vertical axis (in red) applies to the fraction of hospitalized requiring ICU (red solid line) and the ICU beds capacity threshold ($T_i^H/N_i$, dashed red line). The scale on the right vertical axis (in black) applies to the infected (blue), quarantined (magenta), recovered (green) and deceased (black). The time scale, on the horizontal axis, is given in days. Panels of regions adopting a lockdown are shaded in red while those of regions relaxing social containment measures are shaded in green. Results are averaged over 10,000 simulations with parameters sampled using a Latin Hypercube technique (see Methods) around their nominal values set as those estimated in the last time window for each region as reported in Supplementary Table 4. Shaded bands correspond to twice the standard deviation. The regions identified with a red label are those where the total hospital capacity is saturated.

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
<figure><img src="image_7.png" /><figcaption>**Fig. 3** Intermittent regional measures. **a** Each of the 20 panels shows the evolution in a different 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 (see Methods) around their nominal values set as those estimated in the last time window for each region as reported in Supplementary Table 4. Shaded bands correspond to twice the standard deviation. Dashed red lines represent the fraction of the population that can be treated in ICU ($T_i^H/N_i$). 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. **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. **c** National evolution of the fraction in the population of infected (blue), quarantined (magenta) and hospitalized requiring ICUs (red) when an intermittent national lockdown is enforced with all regions shutting down when the total number of occupied ICU beds at the national level exceed 20%, reopening when it goes back below 10%. Regional dynamics corresponding to this scenario are shown in Supplementary Fig. 4. All plots are shown with a double scale. The scale on the left vertical axis (in red) applies to the hospitalized requiring ICU and the ICU beds capacity threshold, while the right vertical axis (in black) applies to the infected and quarantined. The time scale, on the horizontal axis, is given in days.</figcaption></figure>

and modelling 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.

## Methods

**Regional and national model.** As a regional model of the COVID-19 epidemic spread, we use the compartmental model shown in Fig. 5, which we found from data analysis and identification trials to be the simplest model structure able to capture the real data. Specifically, we constructed the model by testing how different configurations of the links among its compartments affected the model ability to capture the available data.
<table><caption>Table 2 Comparison of each of the simulated scenarios.</caption><thead><tr><th>Simulation</th><th>Total cases</th><th>Total deaths</th><th>Maximum hospitalized</th><th>Days over hospital's capacity (nation)</th><th>Regions over hospital's capacity</th><th>Economic cost [M€]</th></tr></thead><tbody><tr><td>All regions but Lombardy are locked down (Fig. 2)</td><td>10,550,000 ± 146,084</td><td>1,196,063 ± 97,122</td><td>137,640 ± 10,249</td><td>75.8 ± 2.7</td><td>3</td><td>503,355 ± 0</td></tr><tr><td>Intermittent regional measures (Fig. 3a, b)</td><td>1,986,601 ± 76,184</td><td>173,637 ± 3911</td><td>2801 ± 170</td><td>0 ± 0</td><td>0</td><td>509,142 ± 6606</td></tr><tr><td>Intermittent national measure (Fig. 3c, S4)</td><td>2,162,539 ± 194,929</td><td>205,261 ± 10,854</td><td>4481 ± 277</td><td>0 ± 0</td><td>3</td><td>562,373 ± 12,809</td></tr><tr><td>Intermittent regional measures with increased testing (Fig. 4)</td><td>1,590,459 ± 69,118</td><td>128,644 ± 2690</td><td>2057 ± 102</td><td>0 ± 0</td><td>0</td><td>366,514 ± 12,258</td></tr></tbody></table>

Metrics to evaluate the impact over 1 year of each of the simulated scenarios are reported showing the effectiveness of the intermittent regional measures shown in Figs. 3 and 4 in avoiding any saturation of the regional health systems while mitigating the impact of the epidemic. We report the average values ±1 standard deviation calculated over 10,000 repetitions of each simulation, where the parameter values are sampled using a Latin Hypercube technique centred at the nominal parameter values reported in Supplementary Table 4.

The full model equations describing the dynamics of susceptible ($S_i$), infected ($I_i$), quarantined ($Q_i$), hospitalized ($H_i$), recovered ($R_i$) and deceased ($D_i$) are

$$ \dot{S}_i = -\rho_i \beta \frac{S_i I_i}{N_i}, \quad (1) $$

$$ \dot{I}_i = \rho_i \beta \frac{S_i I_i}{N_i} - \alpha_i I_i - \psi_i I_i - \gamma I_i, \quad (2) $$

$$ \dot{Q}_i = \alpha_i I_i - \kappa_i^H Q_i - \eta_i^Q Q_i + \kappa_i^Q H_i \quad (3) $$

$$ \dot{H}_i = \kappa_i^H Q_i + \psi_i I_i - \eta_i^H H_i - \zeta_i H_i - \kappa_i^Q H_i \quad (4) $$

$$ \dot{D}_i = \zeta_i H_i, \quad (5) $$

$$ \dot{R}_i = \gamma I_i + \eta_i^Q Q_i + \eta_i^H H_i \quad (6) $$

where $\beta$ and $\gamma$ are the infection and recovery rate, respectively, which are assumed to be the same for all regions as COVID-19 is transmitted from person to person and there is no parasite vector or evidence of environmental parameters significantly altering its infection rate, $\rho_i \in [0, 1]$ is a parameter modelling the effects of social distancing measures in the $i$-th region, $\alpha_i$ is the rate of infected that are detected and quarantined, $\psi_i$ is the rate of infected that needs to be hospitalized, $\eta_i^Q$ is the rate of quarantined who recover, $\eta_i^H$ is the fraction of hospitalized who recover, $\kappa_i^Q$ is the rate of hospitalized that is transferred to home isolation, $\kappa_i^H$ is the rate of quarantined who need to be hospitalized, and $\zeta_i$ is the mortality rate that was shown from data analysis (see Supplementary Notes) to be a function of the ratio between $H_i$ and the maximum number, say $T_i^H$, of patients that can be treated in ICU at the hospitals in $i$-th region. $N_i$ is the actual population in region $i$, i.e., the resident population without those removed because they had been quarantined, hospitalized, deceased or had recovered.

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

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

$$ \dot{S}_i = - \sum_{j=1}^{M} \sum_{k=1}^{M} \rho_j \beta \phi_{ij}(t) S_i \frac{\phi_{kj}(t) I_k}{N_j^p}, \quad (7) $$

$$ \dot{I}_i = \sum_{j=1}^{M} \sum_{k=1}^{M} \rho_j \beta \phi_{ij}(t) S_i \frac{\phi_{kj}(t) I_k}{N_j^p} - \alpha_i I_i - \psi_i I_i - \gamma I_i, \quad (8) $$

$$ \dot{Q}_i = \alpha_i I_i - \kappa_i^H Q_i - \eta_i^Q Q_i + \kappa_i^Q H_i, \quad (9) $$

$$ \dot{H}_i = \kappa_i^H Q_i + \psi_i I_i - \eta_i^H H_i - \kappa_i^Q H_i - \zeta_i (H_i / T_i^H) H_i, \quad (10) $$

$$ \dot{D}_i = \zeta_i (H_i / T_i^H) H_i, \quad (11) $$

$$ \dot{R}_i = \gamma I_i + \eta_i^Q Q_i + \eta_i^H H_i \quad (12) $$

$$ N_i^p = \sum_{k=1}^{M} \phi_{ki}(t) (S_k + I_k + R_k) \quad (13) $$

where in addition to the parameters and states described above, we included the fluxes $\phi_{ij}(t)$ between regions; $\phi_{ij}(t) : \mathbb{R} \to [0, 1]$ denoting the ratio of people from

region $i$ interacting with those in region $j$ at time $t$, such that $\sum_j \phi_{ij}(t) = 1$. Note that, as a result of the identification procedure illustrated in Supplementary Notes, in Eqs. (10) and (11) the mortality rate $\zeta$ is expressed as a function of the saturation of the regional health systems whose expression is given in Supplementary Notes.

**Model parameterization from real data and model validation.** 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.

As a compromise between the estimates reported in the literature on COVID-19$^{2,14}$ in Italy (see Supplementary Table 5), we set $\beta = 0.4$ and $\gamma = 1/14$ [days$^{-1}$] for all regions. We make the ansatz that parameters remain constant over time intervals $T_k$ but do not assume the number or 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 the parameters in each region from the first date when the number of deceased and the number of recovered is greater than or equal to 10.

Note that the official data for the COVID-19 epidemic$^{22}$, 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 ($\tilde{Q}_i$), hospitalized ($\tilde{H}_i$), deceased ($\tilde{D}_i$) and the daily number of individuals that recovered from those who were previously hospitalized or quarantined, say $\tilde{R}_i^Q$. To fit the model to these data, we discretize and rewrite Eqs. (1)–(6) for each region ($i = 1, ..., 20$) as (dropping the subscripts to the parameters for notational convenience)

$$ \dot{S}_i(t+1) = \dot{S}_i(t) - \rho \beta \frac{\dot{S}_i(t) \dot{I}_i(t)}{N_i(0) - \tilde{Q}_i(t) - \tilde{H}_i(t) - \tilde{D}_i(t)} \quad (14) $$

$$ \dot{I}_i(t+1) = \dot{I}_i(t) + \rho \beta \frac{\dot{S}_i(t) \dot{I}_i(t)}{N_i(0) - \tilde{Q}_i(t) - \tilde{H}_i(t) - \tilde{D}_i(t)} - \gamma \dot{I}_i(t) - \tau \dot{I}_i(t) \quad (15) $$

$$ \dot{C}_i(t+1) = \tilde{C}_i(t) + \tau \dot{I}_i(t) \quad (16) $$

$$ \dot{Q}_i(t+1) = \tilde{Q}_i(t) + \alpha \dot{I}_i(t) - \eta^Q \tilde{Q}_i(t) - \kappa^H \tilde{Q}_i(t) + \kappa^Q \tilde{H}_i(t) \quad (17) $$

$$ \dot{H}_i(t+1) = \tilde{H}_i(t) + \psi \dot{I}_i(t) - \eta^H \tilde{H}_i(t) + \kappa^H \tilde{Q}_i(t) - \kappa^Q \tilde{H}_i(t) - \zeta \tilde{H}_i(t) \quad (18) $$

$$ \dot{R}_i^Q(t+1) = \dot{R}_i^Q(t) + \eta^Q \tilde{Q}_i(t) + \eta^H \tilde{H}_i(t) \quad (19) $$

$$ \dot{D}_i(t+1) = \tilde{D}_i(t) + \zeta \tilde{H}_i(t) \quad (20) $$

where measured quantities are denoted by a tilde and estimated state variables by a hat and $\tau = \alpha + \psi$. Here, $C_i = Q_i + H_i + D_i + R_i^Q$ 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 Eqs. (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 Eqs. (17)–(20), as described in detail in Supplementary Notes. Note that, as the actual number of infected is not known$^{9,29}$, 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.071) between the value of the mortality rate, parameter $\zeta_i$ in model (1)–(6), and the saturation of the regional health system
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

**Fig. 4** Intermittent regional measures with increased COVID-19 testing capacity. **a** Each of the 20 panels shows in a double scale plot the evolution in the region named above the panel 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 (see Methods) around their nominal values set as those estimated in the last time window for each region as reported in Supplementary Table 4. Shaded bands correspond to twice the standard deviation. Dashed red lines represent the fraction of the population that can be treated in ICU ($T_{N}/N_{i}$). 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. Regional COVID-19 testing capacities are assumed to be increased by a factor 2.5 (see Methods) with respect to their current values. **b** National evolution of the fraction of infected (blue), quarantined (magenta) and hospitalized requiring ICUs (red) obtained by summing those in each of the 20 regions adopting intermittent regional measures. All plots are shown with a double scale. The scale on the left vertical axis (in red) applies to the hospitalized requiring ICU and the ICU beds capacity threshold, while the right vertical axis (in black) applies to the infected and quarantined subjects.

represented by the ratio between the number of hospitalized in that region ($H_i$) and the total number of available hospital beds in intensive care ($T_{N,i}^H$)
(See Supplementary Notes and Supplementary Fig. 14 for further details and function estimation).

Validation is carried out by using the parameterized model to capture the available data for each window showing a mean squared 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 disease as discussed in Supplementary Notes.

**Cost estimation.** 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
20. Van Mieghem, P., Omic, J. & Kooij, R. Virus spread in networks. *IEEE/ACM Trans. Netw.* **17**, 1–14 (2008).

21. Wang, W. & Zhao, X.-Q. An epidemic model in a patchy environment. *Math. Biosci.* **190**, 97–112 (2004).

22. Dipartimento della Protezione Civile. Dipartimento della Protezione Civile-Presidenza del Consiglio dei Ministri official COVID-19 data repository. http://github.com/pcm-dpc/COVID-19/tree/master/dati-andamento-nazionale (2020).

23. INPS. Dati al 10 maggio su Cassa integrazione ordinaria, assegno ordinario, richieste di pagamento SR41 e Cassa integrazione in deroga. https://www.insp.it/nuovoportaleinps/default.aspx?itemdir=53641 (2020).

24. ISTAT. Reddito netto regioni e tipo di comune. http://dati.istat.it/Index.aspx?QueryId=22919 (2020).

25. INPS. Osservatori statistici e altre statistiche, Indennità 600 euro. https://www.insp.it/docallegatiNP/Mig/Dati_analisi_bilanci/StatInBreve_indennita600.pdf (2020).

26. SVIMEZ Associazione per lo sviluppo dell'industria nel mezzogiorno, L'impatto economico e sociale del Covid-19: Mezzogiorno e Centro-Nord. http://lnx.svimez.info/svimez/wp-content/uploads/2020/04/svimez_impatto_coronavirus_bis.pdf (2020).

27. Meidan, D. et al. Alternating quarantine for sustainable epidemic mitigation. Preprint at https://arxiv.org/abs/2004.01453 (2020).

28. Stolerman, L. M., Coombs, D. & Boatto, S. SIR-network model and its application to dengue fever. *SIAM J. Appl. Math.* **75**, 2581–2609 (2015).

29. Li, R. et al. Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV2). *Science* (2020).

30. Helton, J. C. & Davis, F. J. Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems. *Reliab. Eng. Syst. Saf.* **81**, 23–69 (2003).

## Acknowledgements

M.d.B., M.C. and A.G. wish to acknowledge support from the European Project FET-OPEN Cosy-Bio. M.d.B., M.C., P.D.L. and F.L.I. acknowledge support from the research project PRIN 2017 "Advanced Network Control of Future Smart Grids" funded by the Italian Ministry of University and Research (2020-2023). P.D.L. and F.D.R. acknowledge support from the Compagnia di San Paolo, Istituto Banco di Napoli—Fondazione for supporting their research under the grant STAR 2018, project ACROSS.

## Author contributions

M.d.B. with support from F.D.R. and P.D.L., D.L., G.R. and F.L.I., who contributed equally to this work, designed the research; F.D.R., A.D.M, P.D.L. and D.S. carried out

the model parameterization; F.L.I., C.C. and A.G. analysed the data and estimated the inter-regional flux matrix; F.D.L., M.C. and R.C. wrote and checked the numerical code used for all simulations; D.L., M.C. and R.C. investigated feedback strategies for mitigation and containment; F.D.L., D.S., M.C., C.C., A.G. and R.C. carried out the numerical simulations with inputs from the rest of the authors; M.d.B. wrote the manuscript with inputs from F.D.R., P.D.L., A.D.M., D.L., F.L.I. and G.R.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41467-020-18827-5.

Correspondence and requests for materials should be addressed to M.d.B.

Peer review information *Nature Communications* thanks Sally Blower and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020