The network model of Italy we adopt in this study is, for i = 1, …, 20,7 \documentclass[12pt]{minimal} 				\usepackage{amsmath} 				\usepackage{wasysym} 				\usepackage{amsfonts} 				\usepackage{amssymb} 				\usepackage{amsbsy} 				\usepackage{mathrsfs} 				\usepackage{upgreek} 				\setlength{\oddsidemargin}{-69pt} 				\begin{document}$$\dot S_i = - \mathop {\sum}\limits_{j = 1}^M {\mathop {\sum}\limits_{k = 1}^M {\rho _j} } \beta \phi _{ij}\left( t \right)S_i\frac{{\phi _{kj}\left( t \right)I_k}}{{N_j^p}},$$\end{document}S˙i=−∑j=1M∑k=1MρjβϕijtSiϕkjtIkNjp,8 \documentclass[12pt]{minimal} 				\usepackage{amsmath} 				\usepackage{wasysym} 				\usepackage{amsfonts} 				\usepackage{amssymb} 				\usepackage{amsbsy} 				\usepackage{mathrsfs} 				\usepackage{upgreek} 				\setlength{\oddsidemargin}{-69pt} 				\begin{document}$$\dot I_i = \mathop {\sum}\limits_{j = 1}^M {\mathop {\sum}\limits_{k = 1}^M {\rho _j} } \beta \phi _{ij}\left( t \right)S_i\frac{{\phi _{kj}\left( t \right)I_k}}{{N_j^p}} - \alpha _iI_i - \psi _iI_i - \gamma I_i,$$\end{document}I˙i=∑j=1M∑k=1MρjβϕijtSiϕkjtIkNjp−αiIi−ψiIi−γIi,9 \documentclass[12pt]{minimal} 				\usepackage{amsmath} 				\usepackage{wasysym} 				\usepackage{amsfonts} 				\usepackage{amssymb} 				\usepackage{amsbsy} 				\usepackage{mathrsfs} 				\usepackage{upgreek} 				\setlength{\oddsidemargin}{-69pt} 				\begin{document}$$\dot Q_i = \alpha _iI_i - \kappa _i^HQ_i - \eta _i^QQ_i + \kappa _i^QH_i,$$\end{document}Q˙i=αiIi−κiHQi−ηiQQi+κiQHi,10 \documentclass[12pt]{minimal} 				\usepackage{amsmath} 				\usepackage{wasysym} 				\usepackage{amsfonts} 				\usepackage{amssymb} 				\usepackage{amsbsy} 				\usepackage{mathrsfs} 				\usepackage{upgreek} 				\setlength{\oddsidemargin}{-69pt} 				\begin{document}$$\dot H_i = \kappa _i^HQ_i + \psi _iI_i - \eta _i^HH_i - \kappa _i^QH_i - \zeta \left( {H_i{\mathrm{/}}T_i^H} \right)H_i,$$\end{document}H˙i=κiHQi+ψiIi−ηiHHi−κiQHi−ζHi/TiHHi,11 \documentclass[12pt]{minimal} 				\usepackage{amsmath} 				\usepackage{wasysym} 				\usepackage{amsfonts} 				\usepackage{amssymb} 				\usepackage{amsbsy} 				\usepackage{mathrsfs} 				\usepackage{upgreek} 				\setlength{\oddsidemargin}{-69pt} 				\begin{document}$$\dot D_i = \zeta \left( {H_i{\mathrm{/}}T_i^H} \right)H_i,$$\end{document}D˙i=ζHi/TiHHi,12 \documentclass[12pt]{minimal} 				\usepackage{amsmath} 				\usepackage{wasysym} 				\usepackage{amsfonts} 				\usepackage{amssymb} 				\usepackage{amsbsy} 				\usepackage{mathrsfs} 				\usepackage{upgreek} 				\setlength{\oddsidemargin}{-69pt} 				\begin{document}$$\dot R_i = \gamma I_i + \eta _i^QQ_i + \eta _i^HH_i$$\end{document}R˙i=γIi+ηiQQi+ηiHHi13 \documentclass[12pt]{minimal} 				\usepackage{amsmath} 				\usepackage{wasysym} 				\usepackage{amsfonts} 				\usepackage{amssymb} 				\usepackage{amsbsy} 				\usepackage{mathrsfs} 				\usepackage{upgreek} 				\setlength{\oddsidemargin}{-69pt} 				\begin{document}$$N_i^p = \mathop {\sum}\limits_{k = 1}^M {\phi _{ki}} \left( t \right)\left( {S_k + I_k + R_k} \right)$$\end{document}Nip=∑k=1MϕkitSk+Ik+Rkwhere in addition to the parameters and states described above, we included the fluxes ϕij(t) between regions; \documentclass[12pt]{minimal} 				\usepackage{amsmath} 				\usepackage{wasysym} 				\usepackage{amsfonts} 				\usepackage{amssymb} 				\usepackage{amsbsy} 				\usepackage{mathrsfs} 				\usepackage{upgreek} 				\setlength{\oddsidemargin}{-69pt} 				\begin{document}$$\phi _{ij}\left( t \right):{\Bbb R} \to \left[ {0,1} \right]$$\end{document}ϕijt:R→0,1 denoting the ratio of people from region i interacting with those in region j at time t, such that \documentclass[12pt]{minimal} 				\usepackage{amsmath} 				\usepackage{wasysym} 				\usepackage{amsfonts} 				\usepackage{amssymb} 				\usepackage{amsbsy} 				\usepackage{mathrsfs} 				\usepackage{upgreek} 				\setlength{\oddsidemargin}{-69pt} 				\begin{document}$$\mathop {\sum}\nolimits_j {\phi _{ij}} \left( t \right) = 1.$$\end{document}∑jϕijt=1. Note that, as a result of the identification procedure illustrated in Supplementary Notes, in Eqs. (10) and (11) the mortality rate ζ is expressed as a function of the saturation of the regional health systems whose expression is given in Supplementary Notes.