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 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 prelockdown levels (see Supplementary Fig. 1) or all regions were to relax their current restrictions concurrently (see Supplementary Fig. 2). 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 (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_i^H{\mathrm{/}}N_i$$\end{document}TiH/Ni, 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.