key: cord-0116764-h4s41yew
authors: Parino, Francesco; Zino, Lorenzo; Porfiri, Maurizio; Rizzo, Alessandro
title: Modelling and predicting the effect of social distancing and travel restrictions on COVID-19 spreading
date: 2020-10-12
journal: nan
DOI: nan
sha: 949b4abeba094cd51692ca724b1339d14133a7a9
doc_id: 116764
cord_uid: h4s41yew

To date, the only effective means to respond to the spreading of COVID-19 pandemic are non-pharmaceutical interventions (NPIs), which entail policies to reduce social activity and mobility restrictions. Quantifying their effect is difficult, but it is key to reduce their social and economical consequences. Here, we introduce a meta-population model based on temporal networks, calibrated on the COVID-19 outbreak data in Italy and apt to evaluate the outcomes of these two types of NPIs. Our approach combines the advantages of granular spatial modelling of meta-population models with the ability to realistically describe social contacts via activity-driven networks. We provide a valuable framework to assess the viability of different NPIs, varying with respect to their timing and severity. Results suggest that the effects of mobility restrictions largely depend on the possibility to implement timely NPIs in the early phases of the outbreak, whereas activity reduction policies should be prioritised afterwards.

Following the first report of the novel coronavirus (SARS-CoV-2) in Wuhan, COVID-19 has risen above 43 million cases and 1, 157, 509 reported deaths as of 27th October 2020 [1] ). The ongoing pandemic quickly reached Europe during February and March 2020, forcing most of the countries to implement unprecedented non-pharmaceutical interventions (NPIs) to fight the spread [2, 3, 4, 5] . Some of these interventions promote policies to reduce human-to-human interactions, for example by enforcing social distancing, halting nonessential activities, closing schools, and banishing large gatherings [2, 5] . Others limit human mobility by means of travel restrictions and bans [6] . Due to the considerable economic and social cost associated with the implementation of both of these types of policies [7, 8, 9] ), it is crucial to assess their effectiveness. Mathematical and computational epidemic models are key to accurately evaluate a wide range of what/if scenarios, predicting the evolution of the pandemic for different choices of NPIs [10, 11, 12, 6, 13, 14, 15, 16, 17] .

One of the fundamental aspects of the spread of infectious diseases is its spatial evolution and the concurrent role of human mobility patterns [18, 19, 20] . Extensive studies on mobility within the COVID-19 pandemic revealed that population movements are among the main drivers of the spatial spreading of the outbreak [4, 21] . Network structures have emerged as a powerful framework to encapsulate such mobility patterns within mathematical models of epidemics, especially by means of meta-population models [22] . This modelling paradigm is based on the definition of a set of communities (Provinces, Counties, or Regions), connected by a network that captures daily short-range commuting and long-range mobility.

Different from typical meta-population models that assume homogeneous mixing within communities [22, 23] , we propose a network structure that accounts for the inherent, heterogeneous and time-varying nature of human interactions [24, 25] , together with behavioural changes in response to the pandemic evolution [26, 27] . To this aim, individuals interact on the basis of a mechanism inspired by activity-driven networks (ADNs) [28, 29] . Our model comprehends two key aspects of social communities, mobility patterns and temporal, heterogeneous networks of contacts. Within this meta-population model, we incorporate a variation of a susceptible-infected-removed (SIR) epidemic process [30] . Such an epidemic process allows to capture several key features of COVID-19, like the existence of latency periods and the delay in the official reporting of infections and deaths. We calibrate the model on epidemic data from the Italian COVID-19 outbreak [31] , to examine different scenarios that evaluate the spatial effects of NPIs. In particular, we explore the interplay between reduction in social activity and mobility restrictions. At the modelling level, the former mechanism acts upon the network of contacts, while the latter modifies mobility patterns between communities. Our findings reveal that the timing of the interventions is essential toward effective implementations. We conclude that mobility restrictions should be applied at the early stage of the epidemic and coupled with appropriate policies to reduce social activity. Surprisingly, the impact of mobility restrictions is spatially heterogeneous. For the Italian outbreak, this results in a greater benefit for Southern regions, that is, those located far from the initial outbreak. The overall effect of early travel restrictions in these areas leads to 12% of reduction in the total number of deaths. We also examine differential interventions among age cohorts, determining that the application of severe restrictions only to the most vulnerable age cohorts would not be sufficient to effectively reduce the deaths toll. Different phenomena are observed upon the uplifting of containment measures, with the contribution of mobility restrictions being negligible. In this phase of the fight against the epidemic, policies limiting social activity (for instance, by enforcing the use of face masks or social distancing) bestow the main benefits in mitigating resurgent outbreaks.

We consider a population of n individuals partitioned into a set H = {1, . . . , K} of communities, located in bounded geographical areas (administrative divisions, such as Regions, Provinces or Municipalities), where n h is the number of individuals in the hth community. Communities are connected through a weighted graph that models travel paths between them. The weight matrix W ∈ [0, 1] K×K , called routing matrix, is a matrix with with non-negative entries, zeros on the main diagonal and row sums equal to 1, such that W hk is the fraction of members of community h that move to community k per unit-time.

Individuals interact according to a mechanism inspired by ADNs [28, 29] , which accounts for the inherent, heterogeneous propensity of humans to interact with others. Specifically, individuals are divided into P baseline activity classes 0 < a 1 < a 2 < · · · < a P ≤ 1, where the activity a i of individuals in the ith class quantifies their nominal propensity to interact with others. At each time-step and for each activity class i, a randomly selected fraction a i of individuals activates and generates interactions with others. This fraction of the population is called active. Active individuals may generate interactions within the community where they are located, or they may travel and interact in other communities. This aspect is included in a mobility parameter b ∈ [0, 1] that quantifies the baseline fraction of the active population that commutes to other communities; the commuting unfolds according to the routing matrix W (Fig. 1) . The remaining fraction 1 − b of the active population does not commute and interacts locally with individuals randomly selected within their community. We assume that the P activity classes are equally distributed in different communities. Specifically, we introduce the activity distribution vector η ∈ [0, 1] P such that the expected number of individuals with activity class a i in the hth community is equal to the product n h η i . We finally introduce a parameter m ≥ 0 that captures the individuals' interaction rate, that is, the average number of interactions generated by each active individual in a time-unit.

Two parameters α ∈ [0, 1] and β ∈ [0, 1] are introduced to model NPIs. The former, α, models individuals' self-isolation due to the awareness of the disease spreading. In the model, this corresponds to scaling down the individual activity from its baseline value a i to αa i . The latter, β, captures the effect of mobility restrictions, which are modelled by scaling down the baseline mobility parameter from b to βb. 

The disease progression is modelled according to an extension of the classical SIR model (Fig. 2) , which encapsulates a latency period between contagion and infectiousness, a limited duration of the infectious period, coinciding with the peak of the viral load, and a delay for deaths reporting [2] . Specifically, we adopt a susceptible-exposed-infectious-non-infectious-removed (SEINR) compartmental model (Fig. 2 ). After contagion, infected individuals become initially exposed (E) before spontaneously moving into the infectious (I) compartment with rate ν.

Once the infectious period terminates (with rate µ), individuals transition to the non-infectious compartment (N ), before recovering (or dying) with rate γ, which is represented by the R compartment. The compartment N captures the delay between the end of the infectious period and the reporting of a death. The number of deaths is the most reliable parameter for calibration, given the uncertainty in reporting active infectious cases. The parameters have immediate interpretation: 1/ν is the average latency period of the disease (time from contagion to infectiousness), 1/µ is the average period of communicability (in which individuals are infectious) and 1/γ captures the delay before reported deaths. Hence, 1/µ + 1/γ is the average time from infectiousness to the reported death. The contagion mechanism involves an interaction. We introduce a parameter λ ∈ [0, 1] that captures the fraction of susceptible individuals that become exposed after an interaction with an infectious one. The contagion does not depend only on λ, but also on individual properties (their activity) and network structure, as well as on the prevalence of infectious individuals. We denote by Π h i (t) the contagion function, that is, the fraction of susceptible individuals with activity a i and who belongs to community h that becomes infected at time t, whose expression is detailed in the following.

We consider the generic activity class i and community h ∈ H. Let S h i , E h i , I h i and N h i be the number of susceptible, exposed, infectious and non-infectious individuals of class i in community h, respectively. Clearly,

i is the number of removed individuals in class i and community h. In the thermodynamic limit of large populations, n → ∞, we write the following system of ordinary difference equations:

We now detail the contagion mechanism and derive an explicit expression for Π h i (t). To simplify the notation, we omit time t, that is, Π h i is used to denote Π h i (t). In the thermodynamic limit of large populations n → ∞ and assuming that the epidemic prevalence is small so that we can neglect the probability of having multiple interactions with infectious individuals at the same time, the quantity Π h i can be written as the sum of four different terms. The first summand accounts for the contagions caused by the fraction αa i (1 − βb) of active individuals from S h i that remains in the hth community and interacts there with infectious individuals; the second summand accounts for the infections caused by the fraction (1 − αβa i b) of S h i that remains in community h and comes in contact there with active infected individuals; the third and the fourth summands account for the contagions of the fraction αβa i b of S h i that is active and moves to other communities interacting with infected individuals or receiving interactions from active infectious individuals in the community they move to, respectively. These four terms yield

where

are the fraction of infectious and active infectious individuals that are present in community h, respectively. The quantitỹ

is the number of individuals who are located in community h, where a := P i=1 η i a i is the average activity of the population.

We calibrate the model to reproduce the COVID-19 outbreak in Italy, setting Provinces as communities, using epidemiological parameters from the medical literature [32, 33, 34] , mobility data by the Italian National Institute of Statistics (ISTAT) [35] , and data of officially reported deaths [31] . Based on available empirical data on social contacts per age groups [36] , we partition the population in two activity classes. The population below 65 years old forms the high activity class, and the population above 65 years old constitutes the low activity class. Different mortality rates are associated with the two classes to estimate the number of deaths, base on serology-informed data [37] . The details of the model calibration can be found in the following.

The Italian territory is divided into K = 107 Provinces, which are chosen as communities, extracting the corresponding population n h from the census data [35] . Provinces are grouped in 20 Regions, gathered in 5 macro-regions: North-West, North-East, Centre, South and Islands (Supplementary material, Sec. S1 and Fig. S1 ). We partition the population into P = 2 activity classes, based on age-stratified data on social contacts [36] , aggregating age-groups that form high or low number of contacts, respectively. Specifically, the former contains people below 65 years old, while the latter gathers people above 65 years old. According to the same study, we set m = 19.77. The baseline activity classes a 1 and a 2 are determined by matching the average number of contacts of the individuals in the classes. The fraction η i of population in each class is determined from the Italian age distribution [35] .

We consider two types of mobility: commuting pattern between Provinces and long-range mobility. The former is directly obtained from the 2011 census data in the ISTAT database [35] , which has been validated and adopted to model mobility in recent works on COVID-19 [12] . Comprehensive data on long-range mobility is not available. We estimate it as follows. For each province, we consider the number of nights spent in accommodation facilities over the period from February to May 2011, which represents the destinations of travellers [35] . Origins are estimated based on the flows between macro-regions [35] . Assuming uniformity within each macro-region, we set the origins proportional to the population of each Province. Finally, W is obtained by combining the two origin-destination matrices (Fig. 3) . The mobility parameter b is estimated as the fraction of population that moves outside their Province, using data from [35] .

NPIs are implemented as follows. At t = t 0 , we set α = β = 1. Then, based on empirical data [38] , we consider a linear decrease along 15 days to reach a value α low . Such a decrease begins on 5th March (day of the enforcement of the first social distancing measures) and ended on 20th March (when a severe lockdown is enacted). Similar, β is reduced to β low . Based on [38] , we determined that this switch occurred on 1st March for macro-regions North-East, North-West and Centre, and on 7th March for South and Islands. The values of α low and β low are identified from epidemic data.

Epidemic parameters are taken from the literature on COVID-19. Specifically, the latency period 1/ν and the infectious period 1/µ are taken from [2] , based on clinical estimations from [39, 40] , respectively; γ is the inverse of the difference between the average time from infectiousness to the reported death [41] and 1/µ. The infection probability λ depends on the model of social interactions. Hence, we identify it from real-world data. Table 1 reports the parameters used in our simulations.

We calibrated our model by fitting the temporal evolution of the reported deaths, during the COVID-19 outbreak in Italy. Data at the Regional level were retrieved from the official Italian Dipartimento della Protezione Civile [31] database. This database starts on 24th February, and we had extended it backward in time for 20 days (until 4th February). We filled with zero deaths the section of the database from 4th February to 20th February, and we manually corrected the database to include seven deaths that were not reported therein in the period from 20th February to 23rd February (Supplementary material, Sec. S2). To calibrate the model we focused on the period from 4th February (denoted as t 0 ) to 18th May (denoted as t end ), namely until the first relaxation of NPIs. To enhance the reliability of the data we applied a weekly moving average.

Using the SEINR epidemic progression model, we computed the deaths in Province h for activity class i as a fraction of the removed individuals R h i (t), according to the class fatality ratio f 1 = 0.045% and f 2 = 5.6%. The latter was inferred from a serology-informed estimate performed on age-stratified data from Geneva, Switzerland [37] , scaled on the Italian age-distribution using census data [35] . Since we had no access to information about the initial number of exposed E h i (t 0 ) or infected I h i (t 0 ) individuals, such initial conditions needed to be identified. For each Province h and activity class i, we initialised the number of exposed and infected as a fraction k 1 and k 2 of the total reported cases at the end of the observation time (24th June) C h (t end ) from the official database [31] :

where k 1 and k 2 were identified together with the other parameters.

The parameter identification was formulated as a minimisation problem, solved by means of a dual-annealing procedure [42] . Specifically, we defined the cost function c as the weighted sum of the squared error between the number of deaths predicted by the model and the Regional real data, normalised with respect to the maximum number of deaths in the Region. To this aim, we defined the set of Regions R and the partition of Provinces into Regions as the function π : H → R, such that π(h) = r if and only if Province h was located in Region r. For each r ∈ R, we introduced:

and the cost function as the sum of c r weighted with the total number of deaths in the Region r

Here, D r (t) indicates the reported deaths in Region r, D r M A (t) the weekly moving average andD r M A = max t D r M A (t) its maximum value. Using the fatality ratio, the model predicts f R h (t) deaths in the Province h at time t. Fig. 4 shows real and simulated time-series with the identified parameters. Table 1 of the main document. The red area denotes the time interval before the implementation of activity reduction. The green area denotes the 15-days in which α decreases linearly from 1 to α low = 0.176. In the white area, activity is reduced to α low . The black, dash-dotted vertical line denotes the implementation of mobility restrictions. Orange curves illustrate the predictions from the model, blue curves are actual deaths data from [31] .

Here, we elucidate the role of NPIs in halting the spread of COVID-19. We aim at disentangling the contribution of the two most common kinds of interventions: reduction of individuals' activity, through lockdown or social distancing, and enforcement of mobility restrictions. We take as a reference the NPIs implemented in Italy (detailed in the Supplementary material, Sec. S2) and identify the NPIrelated parameters from available data. The enforcement of lockdown and social distancing policies, gradually enacted during a time-window of two weeks (from 5th March 5 to 20th March) are modelled through a linear decrease of the α parameter from 1 to α low = 0.176. The effect of the nearly complete mobility restrictions between Provinces has been observed from 1st March in the northern macro-regions and from 7th March in the southern ones [38] . We model these restrictions by setting the mobility parameter to β low = 0 on the corresponding dates.

We start from investigating the effect of mobility restrictions in combination with activity reduction policies (Fig. 5a) . We simulate the mobility restrictions as being applied on 4th February, that is, almost one month earlier than the actual date. We compare the number of deaths over the time-window that ranges from 4th February to the date of the first relaxation of NPIs in Italy (18th May). We observe that the effect of mobility restrictions becomes significant for intermediate levels of activity reduction policies (that is, 0.3 < α < 0.7). On the other hand, a negligible effect of mobility restrictions is registered for milder levels of activity reduction policies (α > 0.7) and for extremely severe activity reductions (α < 0.3). The latter, counter-intuitive finding, is due to a balance between the increased number of deaths in some northern macro-regions (close to the initial outbreak) and the decrease of deaths in others (Supplementary material, Fig. S5 )

To detail this mechanism, we examine the number of deaths in each macroregion, using different levels of mobility restrictions and setting the activity reduction to the lockdown level, α low = 0.176 (Fig. 5b) . Our results suggest that the impact of mobility restrictions is strongly dependent on their geographical location, and it vanishes if not timely implemented. Notably, we find that the timely implementation of severe mobility restrictions would have reduced the number of deaths by more than 12% in the Islands macro-region (that is, far from where the outbreak was initially located) over the duration of severe NPIs. Such an advantage becomes smaller and smaller as the considered macro-regions are closer to the initial location of the outbreak. Paradoxically, mobility restrictions becomes even slightly detrimental if applied in the North-West macro-region, where the outbreak started. This is due to the commute of infected individuals from the In (a), we illustrate the interplay between the two NPI mechanisms, assuming an earlier application of mobility restrictions, and activity reduction applied at the original date (5th March). We consider different combinations of levels of activity reduction and mobility restrictions. The heat map codes the reduction of deaths with respect to a scenario where mobility restrictions are not applied. Panel (b) details the effect of earlier mobility restrictions at a macro-regional level, assuming a level of activity reduction as per the lockdown phase (α low = 0.176). The inset illustrates the effect corresponding to the application of the same restrictions at the actual date (1st March). Figure 6 : Effect of activity reduction and mobility restrictions. We consider the total number of deaths over a time-window of 104 days from the beginning of the simulation (4th February) to the time corresponding to the relaxation of the most severe NPIs in Italy (18th May). We investigate three different intervention scenarios. In panels (a,d), both mobility restrictions and activity reduction are applied at the actual application dates. In panels (b,e) both strategies are hypothetically implemented earlier by 15 days. In panels (c,g) mobility restrictions are further set earlier on 4th February, while keeping activity reduction applied earlier by only 15 days. In panels (a-c), the province of Sud Sardegna are illustrated as an example of the positive effect of early mobility restrictions. In panels (d-f), the province of Bergamo look unaffected by mobility restrictions. In (g), we illustrate the classification of Provinces in affected and unaffected by mobility restrictions, considering the scenarios relative to panels (c,f).

most affected Provinces to the rest of the Provinces and of susceptible individuals from less impacted Provinces to the rest of the country. For comparison, we also report death count for the implementation of the same restrictions on 1st March, corresponding to the actual date of their implementation. We observe that the timing of NPIs is essential; an early application of travel restrictions by one month would have saved twice as many lives. Similar results are obtained for the peak of the epidemic incidence (Supplementary material, Figs. S6 and S7). The large geographic variability of the number of deaths is confirmed in Figs. 6a-6f, which depict the total number of deaths for two exemplary Provinces under different timing and intensity of implementation of NPIs. While the Province of Bergamo (in the North-West macro-region), one of the earliest and biggest out-breaks, seems unrelieved by mobility restrictions, the Province of Sud Sardegna (in the Islands macro-region), an area much less affected by the pandemic than the former, would have largely benefited by such an intervention. To deepen this aspect, we factor out the role of the two types of NPIs by performing a non-negative matrix factorisation (NMF) [43] on the outcome of our simulations at the Province level (details in the Supplementary material, Secs. S3 and S4). Specifically, we focus on values of α ranging over a ±50% interval with respect to the value of the lowest activity coefficient α low = 0.176, identified from real-world data during the lockdown, and we simulate the early application of mobility restrictions with different intensity levels and timing. The NMF allows to partition the Italian Provinces in two sets. The first (in green in Fig. 6g ) comprises Provinces where timely implemented mobility restrictions are effective in reducing epidemic prevalence (for example, Sud Sardegna). The second set (in brown in Fig. 6g) contains Provinces for which mobility restrictions have instead a negligible impact. Predictably, most of the Provinces in the North-West (where the outbreak started) are unaffected by mobility restrictions, while the majority of Provinces in South and Islands would benefit from an early implementation of such restrictions.

The NMF allows to detect some important exceptions. For instance, the Provinces of Varese and Monza (close to the Milan metropolitan area) would have benefited from timely mobility restrictions. We believe that this is due to the initially small number of cases in those two Provinces, and to the large number of daily commuters from those Provinces to the Milan Province and other neighbouring locations, where the Italian outbreak started. Hence, the same dynamics between North and South Italy is documented again over a much smaller spatial scale, between Northern Provinces with larger initial difference in epidemic prevalence. Similar results are observed by performing the NMF for other intervention scenarios (Supplementary material, Fig. S2 ).

Finally, we discuss the possibility of implementing targeted activity reductions that act independently on the two activity classes. This allows to study the effectiveness of differential intervention policies that could aim at strongly reducing social activity for age cohorts that are more at risk of developing severe illness, while implementing mild restrictions for younger people. Instead of a single parameter α, we thus introduce two parameters α 1 and α 2 that measure the activity reduction for the high and the low activity class, respectively. The heat-map in Fig. 7 illustrates the effect of different combinations of α 1 and α 2 on the total number of deaths; the level curves help understand the trade-off in targeting the two classes. We observe that the total number of deaths is mostly determined by the parameter α 1 , that is, the activity reduction for the high activity class. Hence, our results suggest that implementing targeted stay-at-home policies in which severe activity reductions are only enforced on the age cohorts that are more at risk (in (high activity class) (low activity class) Figure 7 : Effect of targeted lockdown strategies. We consider the total number of deaths over a time-window of 104 days from the beginning of the simulation (4th February) to the time corresponding to the relaxation of the most severe NPIs in Italy (18th May). On the two axes, α 1 and α 2 correspond to activity reduction of high and low activity classes, respectively. Level curves are shown to clarify how targeted interventions should be combined to produce the same effect on the total number of deaths. The black dot represents the identified activity reduction in model calibration. our scenario, people over 65 years old) is not sufficient to reduce the overall death toll.

The proposed meta-population model enables the analysis of reopening strategies to uplift restrictions while avoiding resurgent outbreaks. This has recently emerged as a key issue in the control of COVID-19 outbreaks in the medium-to long-term period [44] . We run our calibrated model to simulate the epidemic until the date of intervention uplifting. Then, we vary the values of parameters α and β to account for the uplift of the containment measures. Similar to the previous analysis, we consider a set of different options for the parameters after intervention uplifting and different times for starting the reopening strategies. Specifically, we model the uplifting of the reduction of social activity by varying the parameter α from α low , identified during the lockdown, to a value α = 0.6. Likewise, we describe the uplifting of mobility restrictions by varying the parameter β from 0 (no mobility allowed) to 1 (nominal mobility reinstated). Our results suggest that the effect of enforcing mobility restrictions upon the relaxation of NPIs is negligible and dominated by the activity reduction (Fig. 8) . We evaluate the total number of deaths in a time-window of 60 days after the relaxation date (18th May). Both at the Province level, for which we show the examples of Sud Sargegna (Fig. 8a) and Bergamo (Fig. 8b) , and at the aggregated country level (Fig. 8c) , the contribution of mobility restrictions is little or absent. These results are confirmed by other scenarios with different relaxation dates and a NMF-based analysis (Supplementary material, Secs. S5 and S6 and Figs. S3 and S4). Overall, this evidence indicates that activity reduction in the relaxation of NPIs should be thoughtfully calibrated, trading-off the risk of resurgent outbreaks and the social and economical costs associated with such policies. On the other hand, the further enforcement of mobility restrictions within the country does not seem to be beneficial in the relaxation phase. cal framework to study NPIs and elucidate their impact on epidemic spreading. Specifically, we combined a meta-population model, capturing the spatial distribution of the population and its mobility patterns [22, 23] , with an ADN-based structure, which reflects real-world features of social activity such as heterogeneity [28, 29] and behavioural traits [26, 27] . We explicitly incorporated two types of NPIs: actions aiming at reducing individuals' activity (social distancing, forbidding gatherings and, in general, any measure that curtails the number of contacts favouring the spread of the infection) and policies to restrict individuals' mobility (for instance, through travel bans). Through the lens of our modelling framework, we disentangled the effect of these two types of policies depending on the time of their implementation. We calibrated the model with data on the ongoing COVID-19 outbreak in Italy [31] .

We leveraged the model to explore a wide range of what/if scenarios on spatiotemporal dynamics of COVID-19 spreading for different combinations of NPIs. Our analysis allows to draw interesting conclusions on when and how to apply NPIs to make the fight against the spread more effective. While the level of activity reduction is unequivocally a decisive factor, the impact of mobility restrictions has a more entangled impact. First, we observed that mobility restrictions produce benefits only if applied at the early stage of the outbreak, and only if paired with appropriate activity reduction policies. Moreover, we discovered that the effect of mobility restrictions is strongly dependent on space. In fact, through a nonnegative matrix factorisation technique, we identified two sets of Provinces that are differently affected by mobility restrictions. The first set, mostly consisting of Provinces in the North (where the outbreak initially started), has little or no benefit from mobility restrictions. The Provinces in the second set, instead, would have benefited from early implementation of mobility restrictions. Surprisingly, this set includes some of the Provinces in the north (most affected area). Then, we discussed possible implementation of targeted NPIs, with severe restrictions only for age cohorts that are more at risk of developing severe illness. Our modelling framework brought to light concerning limitations in the implementation of these targeted interventions: although economical reasons may prompt these interventions, their public health value could be limited. Finally, while mobility restrictions are useful in the early stage of the outbreak, their late implementation is ineffective. A different scenario is observed for the relaxation of NPIs, where the level of activity reduction should be carefully and gradually uplifted.

Our study outlines several avenues of future research, which can be pursued leveraging the generality of the heterogeneous meta-population framework proposed in this study. Although NPIs have been homogeneously implemented nationwide through Decrees of the Prime Minister, the model could benefit from the study of heterogeneous implementation (and relaxation) of NPIs between Provinces and even the implementation of targeted mobility restrictions between specific Provinces (through the modification of the routing matrix W ), whose analysis is envisaged for future research. The outcome of such an analysis can inform policymakers on targeted interventions that may reduce social and economical costs while effectively halting the epidemic. Also, other classes may be added to our model to explore other targeted intervention policies, such as those leading to safe schools reopening. The introduction of a community that represents the rest of the world would enable the study of the impact of closures of national borders. While we considered a simple model for the epidemic progression, additional compartments and transitions may be added to capture hospitalisation or testing [45] , and used to analyse different what/if scenarios. Finally, the simplicity of our mathematical framework may be conducive to a rigorous analytical treatment, involving for example the computation of the epidemic threshold, toward providing further insight into the effect of NPIs on the epidemic spreading.

Epidemiological data are available at [31] ; geographic, mobility and census data at [35] .

The code used for the simulations is available at https://gitlab.com/PoliToC omplexSystemLab/adn-metapopulation-2020. a first draft of the manuscript. A.R. and M.P. supervised the research and consolidated the manuscript in its present submission. All the authors contributed to the interpretation and analysis of the results and to reviewing the current submission of the manuscript. Supplementary references 18 27 S1 Italian geographic organisation 28 The Italian territory is divided into four levels of administrative entities. 29 At the finer level, there are 7, 903 municipalities (local administrative units 

19, but they were not reported in [1] , which starts reporting data from Here, we explain the non-negative matrix factorisation (NMF) through the 84 scenario of early implementations of NPIs, whose results are shown in Fig. S2 . 85 We simulated the total number of deaths, over a 104 days time-window, for 86 p different levels of α, q different levels of β, and for each Province h ∈ H.

The results of these simulations were stored in a set of non-negative matrices 88 P h ∈ R p×q + .

The NMF algorithm approximates each matrix P h as a weighted sum of 90 an arbitrary number of basis matrices, which help understand the effect of 91 NPIs on all the Provinces. We select two basis matrices, denoted as C 1 and 92 C 2 , so that

where k h1 and k h2 are two scalar weights that are specific to each matrix P h .

The identification of both weights and basis matrices is performed through 95 the following steps. First, we normalise the entries of each matrix P h between 96 0 and 1. The application of the approach is presented in Fig. S2 . As an illustration, Figure S2 : Effect of activity reduction and mobility restrictions on the total number of deaths over a time-window of 104 days from the beginning of the simulation (4th February) to the relaxation of the more severe NPIs (18th May). We plot the two basis matrices for three different intervention scenarios. In (a,b), both mobility restrictions and activity reduction are applied as per the actual application dates. In (c,d), hypothetical 15-days early implementation of both strategies. In (e,f), mobility restrictions are enacted even earlier, on 4th February, while activity reduction is applied 15 days in advance with respect to its actual implementation. In (b,d,f), we show the Province's weights related to the two basis matrices, coloured respect to the results of the k-means algorithm. This figure extends Fig. 6 in the main document, which shows two exemplary Provinces. Bergamo (panels (d)-(f)) and at the National level (panels (g)-(i)). Figure S4 : Effect of the relaxation of NPIs, for different post-relaxation levels of activity reduction and mobility restrictions. We consider the number of deaths over a time-window of 60 days, and we show the two basis matrices for three different scenarios. In (a,b), we show a late hypothetical relaxation of the restrictions on 1st July. In (d,c), we consider as relaxation date 3rd June, which corresponds to the complete relaxation of mobility restrictions in Italy. In (e,f), we consider as relaxation date 18th May, corresponding to the first uplifting of NPIs in Italy. In (b,d,f), we show the Provinces' weights related to the two basis matrices. Figure S5 : Effect of early application of mobility restrictions. We consider the total number of deaths over a time-window of 104 days from the beginning of the simulations (4th February) to the time corresponding to the relaxation of the most severe NPIs in Italy (18th May). We plot different combinations of levels of activity reduction and mobility restrictions, implemented on 5th March (actual date) and 4th February (early implementation), respectively. For each level of activity reduction α (row of the heat-map), the effect of mobility restrictions is reported in terms of the percentage difference in the total number of deaths with respect to the corresponding case with no mobility restrictions (β = 1). We show the results for six Provinces, selected as representative examples. This figure extends the results illustrated in Fig. 5 of the main document, which reports predictions aggregated at the country level. Figure S6 : Effect of early application of mobility restrictions on the peak of the epidemic incidence. We consider the maximum number of new infected individuals from the beginning of the simulations (4th February) to the time corresponding to the relaxation of the most severe NPIs (18th May). In (a), we consider the effect of different levels of mobility restrictions enacted on 4th February (almost one month earlier than the actual dates). The figure inset shows the effect of the same interventions applied on March 1 (the actual application date in North Italy). The activity reduction parameter is fixed to the lockdown level α low = 0.176. Results are aggregated at the macro-region level and reported as the variation of the infected individuals at the peak of the epidemic incidence C peak , with respect to the case with no mobility restrictions (β = 1). In (b), we consider different combinations of levels of activity reduction and mobility restrictions, implemented on 5th March (actual date) and 4th February (early implementation), respectively. For each level of activity reduction α (row), the effect of mobility restrictions is reported in terms of the percentage difference in the number of infected individuals at the peak of the epidemic incidence with respect to the corresponding case with no mobility restrictions (β = 1). Figure S7 : Effect of early application of mobility restrictions on the peak of the epidemic incidence for six Provinces, selected as representative examples. We consider the maximum number of new infected individuals from the beginning of the simulations (4th February) to the time corresponding to the relaxation of the most severe NPIs (18th May). We consider different combinations of levels of activity reduction and mobility restrictions, implemented on 5th March (actual date) and 4th February (early implementation), respectively. For each level of activity reduction α (row), the effect of mobility restrictions is reported in terms of the percentage difference in the number of infected individuals at the peak of the epidemic incidence with respect to the corresponding case with no mobility restrictions (β = 1).

Coronavirus Disease (COVID-2019) Situation Reports

The effect of control strategies to reduce social mixing on outcomes of the COVID-19 epidemic in Wuhan, China: a modelling study

Effect of non-pharmaceutical interventions to contain COVID-19 in China

The effect of human mobility and control measures on the COVID-19 epidemic in China

An investigation of transmission control measures during the first 50 days of the COVID-19 epidemic in China

The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak

The impact of COVID-19 on small business outcomes and expectations

Economic and social consequences of human mobility restrictions under COVID-19

A nationwide survey of psychological distress among Chinese people in the COVID-19 epidemic: implications and policy recommendations

Opinion: What models can and cannot tell us about COVID-19

The challenges of modeling and forecasting the spread of COVID-19

Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures

Modelling the impact of testing, contact tracing and household quarantine on second waves of COVID-19

Mathematical models to guide pandemic response

Modelling COVID-19

COVID-19 and SARS-CoV-2. Modeling the present, looking at the future

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

The role of the airline transportation network in the prediction and predictability of global epidemics

The hidden geometry of complex, network-driven contagion phenomena

Multiscale mobility networks and the spatial spreading of infectious diseases

Population flow drives spatio-temporal distribution of COVID-19 in China

Epidemic processes in complex networks

Critical regimes driven by recurrent mobility patterns of reaction-diffusion processes in networks

Epidemic thresholds in dynamic contact networks

Temporal networks

Modelling the influence of human behaviour on the spread of infectious diseases: a review

Effect of individual behavior on epidemic spreading in activity driven networks

Activity driven modeling of time varying networks

Continuous-time discrete-distribution theory for activity-driven networks

Mathematical models in population biology and epidemiology

COVID-19 Italia -Monitoraggio situazione. GitHub

Evolving epidemiology and transmission dynamics of coronavirus disease 2019 outside Hubei province, China: a descriptive and modelling study

Report of the WHO-China joint mission on coronavirus disease 2019 (COVID-19). Geneva; 2020

SARS-CoV-2 viral load in upper respiratory specimens of infected patients

Nazionale di Statistica

Social contacts and mixing patterns relevant to the spread of infectious diseases

Serology-informed estimates of SARS-CoV-2 infection fatality risk in

COVID-19 outbreak response, a dataset to assess mobility changes in Italy following national lockdown. Sci Data

Incubation period of 2019 novel coronavirus (2019-nCoV) infections among travellers from Wuhan, China

Clinical presentation and virological assessment of hospitalized cases of coronavirus disease 2019 in a travel-associated transmission cluster

Incubation Period and Other Epidemiological Characteristics of

Novel Coronavirus Infections with Right Truncation: A Statistical Analysis of Publicly Available Case Data

Efficiency of generalized simulated annealing

Learning the parts of objects by non-negative matrix factorization

Assessing the impact of coordinated COVID-19 exit strategies across

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

COVID-19 Italia -Monitoraggio situazione. GitHub

Gazzetta Ufficiale della Repubblica Italiana

Data clustering: a review

The authors are indebted to Alessandro Vespignani for precious discussion. 

L.Z. and A.R. conceived and designed the research with inputs from all the authors. F.P. performed the parameter identification and the numerical studies and wrote

The authors declare that they have no competing interests.