Impact of urban structure on infectious disease spreading

The ongoing SARS-CoV-2 pandemic has been holding the world hostage for several years now. Mobility is key to viral spreading and its restriction is the main non-pharmaceutical interventions to fight the virus expansion. Previous works have shown a connection between the structural organization of cities and the movement patterns of their residents. This puts urban centers in the focus of epidemic surveillance and interventions. Here we show that the organization of urban flows has a tremendous impact on disease spreading and on the amenability of different mitigation strategies. By studying anonymous and aggregated intra-urban flows in a variety of cities in the United States and other countries, and a combination of empirical analysis and analytical methods, we demonstrate that the response of cities to epidemic spreading can be roughly classified in two major types according to the overall organization of those flows. Hierarchical cities, where flows are concentrated primarily between mobility hotspots, are particularly vulnerable to the rapid spread of epidemics. Nevertheless, mobility restrictions in such types of cities are very effective in mitigating the spread of a virus. Conversely, in sprawled cities which present many centers of activity, the spread of an epidemic is much slower, but the response to mobility restrictions is much weaker and less effective. Investing resources on early monitoring and prompt ad-hoc interventions in more vulnerable cities may prove helpful in containing and reducing the impact of future pandemics.

Human mobility is known to play a central role in the spreading process [15][16][17][18][19][20][21][22][23][24][25] as it was the case for the first wave that hit China following the inter-city mobility from Wuhan 9,26,27 . In a connected planet, rapid world-wide spread is enabled by long-distance air-, land-and sea-transportation among countries and continents, which subsequently expand the disease to new and susceptible area of the world 9,18,20,[26][27][28][29][30][31][32] . While early travel restrictions contribute to delayed disease spread, their utility is much reduced if the disease has incubation periods longer than a day or if there is asymptomatic transmission 3,5,6,8,12,[33][34][35][36][37][38][39] . Long-range mobility is mainly driven by air transportation, but restrictions to international flights have shown limited utility in mitigating the propagation of infectious diseases unless they are applied very early and in a massive manner 6,36,38 .
The majority of work on the role of mobility in regulating the spread of pandemics has focused on long range trips between countries, states and cities. Once restrictions on land, sea and air travel are imposed, and long range transportation networks cease to play a role, the core of contagion keeps developing in the main center of human activity: the densest and most crowded cities of the globe. In urban areas, the predominant transportation modes include vehicular traffic, pedestrian, bike, metros and buses 24,40 . In particular, the mixture of different transportation modes at an urban level, coupled with the spatial organization of socioeconomic activities in residential areas, creates unique mobility fingerprints [41][42][43][44][45][46][47] in which stations, office buildings and stores are continuously packed and close contact is inevitable. Historically, epidemics have influenced planning of cities 48 . In the history of urban planning, e.g. after the black plague or cholera, ventilation of streets, indoor daylight, large public squares and parks became architectural recipes to mitigate the spread of disease, in particular respiratory ailments such as tuberculosis. Such measures helped to reduce the density of cities, created safer and less infectious residential areas and increased the opportunities for physical distancing in public spaces, while allowed for better isolation of affected individuals from the larger urban population [49][50][51] .
Currently more than half the world's population lives in urbanized areas, while forecasts see this number growing to 68% in 2050 52 . A vigorous debate on how to redesign urban spaces to respond to future challenges stemming from the post-pandemic era has emerged 50,53 . While there has been prior work on the role of mobility on viral spread in urban systems 26,27,30 , showing how dense urban areas exhibit more diffuse epidemics of influenza 54 and higher attack rate peaks of SARS-Cov2 30 , scarce attention has been devoted to the intrinsic spatial structure of intra-urban mobility and its effect on the evolution of urban outbreaks. The way in which human mobility affects an epidemic progression is by regulating the population mixing patterns. This is not only an important factor at international and national levels, but it is critical in urban environments where most close contacts happen. This is particularly important in the case of airborne diseases such as COVID-19 or influenza, given that contagions may occur at fortuitous encounters due to shared transportation media, crowded areas, and business activities. Hence, understanding how the spreading of a disease depends on the structural organization of a city mobility and its population mixing patterns is of fundamental importance to devise adequate and tailored containment strategies and surveillance protocols 55 .
Here we investigate how the spatial structure of urban mobility shapes critical features of the local spreading patterns and the efficiency of the potential mitigation measures. We showed in a previous analysis that when it comes to the internal organization of city-wide mobility patterns, urban areas can be classified in a spectrum between centralized-hierarchical and decentralized-sprawled 56 . The level of centralization can be encapsulated in a metric , called the flow-hierarchy. As a first step, we analyze official surveillance data at county scale in the US and human mobility records from the Google COVID-19 Aggregated Mobility Research Dataset. The analysis performed on 22 of the most populated cities of the US indicates that intra-urban mobility hierarchy is connected to higher incidence peaks and faster epidemic growths in cities. We confirm the empirical observations through a compartmental model of epidemic spread on cities with different levels of showing higher peaks and faster outbreaks in more hierarchical cities as compared to sprawled ones, independent of their population density or shape. Moreover, our formulation enables analysis of "what-if " scenarios to assess the effect of mitigation policies, demonstrating that hierarchical cities perform and respond better at containing the epidemics when lockdowns are issued. The model is finally performed repeatedly for higher and lower values of R 0 in order to check for the stability of results in the context of new variants or other diseases, which can show higher or lower infectiousness. Our results are of fundamental interest in the context of global and national surveillance systems, improving risk-perception of epidemics, and in order to prevent new waves of infections caused by new COVID-19 variants and new viruses in general.

Methods
Mobility data. The mobility flows are sourced from the Google COVID-19 Aggregated Mobility Research Dataset, containing anonymized trips weekly aggregated over users who have turned on the Location History setting, which is off by default. The dataset aggregates flows of populations between S2 cells (https:// github. com/ google/ s2geo metry) of approximately 5 km 2 . To produce this dataset, machine learning is applied to logs data to automatically segment it into semantic trips 56 . To provide strong privacy guarantees, all trips were anonymized and aggregated over one week using differential privacy 57 . The research presented here is carried out on the resulting weekly flow data. No individual user data was ever manually inspected, only heavily aggregated flows of large populations were handled. To further confirm the model results, we have additionally used commuting data at zip code level from the Longitudinal Employer-Household Dynamics surveys of the US census office. In particular, the information comes from the LODES dataset concerning commuting in 2016-2017 (see Data Data availability statement for the download link).
Urban mobility and its hierarchy. As mentioned, the space of the urban areas is divided in S2 cells. The weekly trip flows between pairs of areas are incorporated in an Origin-Destination matrix T(t) , whose elements T ij (t) encode the trip flow at week t between the two spatial units i and j. The diagonal terms T ii correspond to movements within the area i. The total flow of a territory is the sum of the trips over all the S2 cells present in the city, T = i,j T ij . In general, the relative trip-flow change between two timestamps t 1 and t 2 is calculated as (T(t 1 ) − T(t 2 ))/T(t 2 ) . Note that when we estimate the reduction after lockdown, t 2 refers to the mobility prior to the restrictions and t 1 after the lockdown. To calculate , a hotspot level must be assigned to each geographical unit as explained in Ref. 56 . In a nutshell, a Lorenz curve is constructed from the trip outflow of each area; the curve contains information on the differences in out-mobility between cells by ordering them in increasing order in terms of outflow and representing the fraction of cells versus the fraction of total outflow. Taking the derivative at (1, 1) establishes a threshold. Those cells over the threshold are classified as level-1 hotspots 58 . The process is applied iteratively with the remaining cells until all of them have an assigned level. Once every cell i is assigned a level L i , trip flows between cells are aggregated in a matrix S ℓm as: The value lies in the range 0 ≤ ≤ 1 with the limiting cases corresponding to flows being distributed to levels far from the considered hotspots, or flows only being concentrated at similar levels. In hierarchical cities, trips are concentrated between mobility hotspots at higher levels (spatially clustered), with taking on values closer to one. In sprawled cities, hotspots display a homogeneous spatial distribution and the flows are distributed more evenly across levels, with taking on a lower value. Empirically measured values for global cities typically lie in the range 0.7 ≤ ≤ 0.95 . While the range appears superficially narrow, the differences track well with significant variability in a variety of urban indicators associated with different levels of public transportation usage, walkability, health indicators and levels of pollution 56 . The metric is also well correlated with considerably more complex composite measures (integrating socioeconomics, land-use, population density among others) such as so-called sprawl indices 59 , and outperforms them as a predictor of urban indicators.
Epidemic data. The epidemic data has been downloaded from two sources (USAFacts and NY Times) listed in the Data availability statement below. Despite, the fact that we have global mobility data, in what is to follow, we focus on cities in the United States. This is done primarily for two reasons: first, US cities show a wide spectrum of configurations from hierarchical to extended or sprawled, with New York city and Atlanta as archetypal examples. Second, due to inherent biases and noise in collection of epidemic data, it is problematic to compare between different countries. We also restrict our analysis to the time-period corresponding to the so-called first wave of infections. For the purposes of our analysis we collected data on the 22 largest metropolitan areas in the United States, in terms of the number of counties and with populations exceeding 2 million inhabitants. The full list of cities and their corresponding values of are shown in SI Table S1. The data resolution is at the level of counties, which in a large city are well populated units with reasonable statistics (SI Tables S2-S23 detail the list of counties considered per OECD urban area).
Metrics of epidemic severity. The epidemic evolution is studied using the incidence curve as basis, which corresponds to the number of new cases detected in a day per capita. In the case of the data, the number of detected individuals depends on the testing policy, while in the models these are all the new infected individuals in each simulated day. The incidence can then be represented by curves I u j (t) for each county j of city u at day t. These incidence curves can be aggregated at full city scale by summing the cases in the counties forming the urban area and dividing by the total population, I u (t) . The first metric will be the maximum value of the curve at city u, I u max . In the case of the empirical data and to avoid fluctuations, we use a 14 days aggregation of the incidence before the peak.
Besides the direct height of the incidence peak, we are interested in how quick is the spreading process initially in every population. To quantify this, we calculate R eff that is the effective reproduction number. This has been done using the methodology developed in 60 and the code available at 61 . This method relies on a Bayesian dynamic inference of the reproductive number based on a stochastic epidemic model informed by an estimation of the serial interval mean. The meaning of R eff is related to the exponential change of the incidence: the curve grows if R eff > 1 , and it decays if R eff < 1 . R eff can be measured at different geographical scales but it has a more practical meaning in localized regions such as cities or counties where it can accurately reflect the contagion patterns. We then define R early as the average over three weeks of the effective reproductive number measured at the early epidemic stages (after the city incidence has gone over a certain case onset: 100 for the empirical data and 10,000 for the models where there is no detection problem). For a robustness test, we applied the same methodology on other types of data from the UK comparing R eff and the R early computed on the reported cases by publication date with cases by specimen date, finding strongly consistent results (see SI Figs. S34 and S35). Finally, in order to define the counties synchronization in the epidemic evolution, after mobility restrictions are implemented in a city we define R p,response as the Pearson correlation coefficient of the consecutive five weeks R eff and mobility reduction of each area from the onset of the first 100 registered cases.

Model description. The simplest version of our model is structured in a meta-population framework with
Susceptible (S), Exposed (E), Infected (I) and Recovered (R) compartments and taking as basic spatial units the S2 cells provided by the mobility data. For a given city, the population P is distributed among the S2 cells i according to their mobility, p i = P ( j T ij )/T where the index j runs over all cells including i and T corresponds to the total trips of the city. Note that while this is not a realistic configuration the goal here for the model is to uncover the effect of mobility hierarchy, and not necessarily to forecast epidemic evolution, which would require more precise data on the population mixing. Usually, only a fraction of people, M, leave the residence area, as for instance in the case of commuting between zip codes, a value of M ≈ 0.4 was measured in US cities (see SI Figs. S6 and S7). A similar value was observed in previous works for commuting in more extensive areas across the world 14,28,62 . Unless otherwise specified, the particular parameters selected for the model are shown in Table 1. Consequently in our model, when there are no travel restrictions, a fraction M = 0.4 of the population of every cell i moves and selects their destination proportional to outflow T ij from i to the other cells j (including T ii ). Non-moving individuals are not excluded from the epidemic dynamics but contacts occur inside their cell of residence. Note that we assume that the flows represent recurrent mobility, hence individuals return home www.nature.com/scientificreports/ after a time τ set at 8 hours. The incorporation of recurrent mobility in epidemic dynamics is done by using a classical approach that encodes the information of the mobility matrices into an effective transition rate 28,62,64 , see SI section S8 for details. Cities in this initial idealized configuration maintain their empirical because the outflows are proportional to those in the empirical network regardless of the value of M.
The main outcome of the model are the epidemic curves for each geographical unit (S2 cells). These cells have to be aggregated to county or city levels to allow for a fair comparison with the results obtained with the epidemiological data. A sketch with the different geographical scales involved and the type of data available at each of them can be found in SI Fig. S37.

Inclusion of lockdowns.
To model the effect of lockdowns on mobility, we need to consider that mobility flows change in both magnitude and destination in response to restrictions, with reductions that can reach up to 80% of the flows compared to baseline levels. For the purposes of simplicity and tractability, for each city we employ two mobility networks: one considering a typical day prior to restrictions with elements T ij as used thus far, and another, after mobility restrictions, with a modified mobility matrix T ′ ij . The total trips per cell i for each of the scenarios are then as the new fraction of individuals traveling outside of their residence under mobility restrictions. Since travel flows in matrix T ′ are usually smaller than flows in T , this results in an overall reduction of mobility.
During stay-at-home orders, part of the population remained isolated or only interacted with individuals in their household; not participating in the global spreading process directly. This effect has been modeled following the literature 63,65,66 , where a fraction X S of susceptible individuals do not participate in the infection process.
To model that restrictive measures are enforced when outbreaks are detected, we assume that restrictions are imposed once the prevalence (the total number of active infected cases per day) reaches a certain threshold value π th . In Fig. 1a, we plot some representative examples for the prevalence π for a variety of mitigation scenarios governed by varying X S . Stronger lockdowns flatten the curves and can even suppress propagation as occurs for X S = 0.9.
Generating instances of a city with different . To produce instances of a single city with different values of , we reshuffle a portion of randomly selected links, for instance the one going from area i to j, T ij , by picking a new destination at random (k) while maintaining the link-weights. In case a link already existed Table 1. Epidemic model parameters obtained from Ref. 63 . www.nature.com/scientificreports/ between the original area i and k, the weights of the links i → j and i → k are interchanged. This procedure creates multiple realizations of the city in terms of the arrangement of mobility flow, while preserving the population distribution and densities, as well as the total number of trips along with the number of destinations per origin. The procedure is applied equally to prior and post confinement scenarios, selecting the same new random destination for every pair T ij and T ′ ij (for a sketch see Fig. 1b).

Results
Empirical local propagation patterns and city organization. The primary response to the early onset of the epidemic were mobility restrictions corresponding to stay-at-home measures, modifying the spatial patterns of urban mobility and disrupting most of the routes used by the virus to propagate. In Fig. 2a-c, we show the observed spatial mobility layout for three cities arranged in increasing levels of centralization in their mobility structure: Atlanta ( = 0.784 ), Chicago ( = 0.868 ) and New York City NYC ( = 0.909 ). In each case, shown as red dots are the top mobility hotspots prior to lockdown. In the case of Atlanta, being a decentralized city, we see that the mobility hotspots were spread across the urban area, whereas Chicago and NYC being more hierarchical had their hotspots clustered in a centralized fashion. The new mobility hotspots that www.nature.com/scientificreports/ emerged after mitigation measures are shown as orange dots, whereas in blue, we show those areas that ceased to be top hotspots. For hierarchical cities such as NYC or Chicago, while the original top hotspots remained in the reduced mobility phase, new areas of mobility activity emerged in an extended fashion, effectively distributing the mobility over larger parts of the city, thus reducing agglomeration. This feature was observed across several cities in the United States, as reflected by plotting in decreasing order, the fractional change in distance between hotspots due to mobility restrictions as shown in SI Fig. S1. The resolution of the epidemic data allows us to run a granular analysis of the impact of lockdowns by investigating the extent to which specific areas of the city drove the epidemic spreading to other parts. To do so, we calculated the average Transfer Entropy, TE , between the incidence curves for each subsets of the city (see SI Section 2 for details), and in Fig. 2a-f we plot the results as a matrix whose rows correspond to the administrative sub-unit, and columns correspond to weeks, starting from the onset of 100 cases in the whole city, to the first week of June. The elements in the matrix are color-coded by the value of TE . For all three cities, a few regions drove the spreading of the infection before lockdown, although the strength of the driving was stronger on average in NYC and Chicago compared to Atlanta (results for four other cities, including three outside the US are reported in SI Fig. S2). Once lockdown started (marked as vertical red lines), the driving became diffused, indicating the predominance of localized spreading in sub-regions with little influence on one another. The localization in spreading appeared in parallel to the equivalent phenomenon in mobility, with a relative increase in self-flows and a decrease of inter-area flows in the administrative units as shown in SI Fig. S3. While the influence of each administrative unit on infection-spreading dissipated in a similar fashion, there is a key difference in how they were synchronized in terms of their response to mobility-mitigation, as we show next.
In Fig. 2g-i, we plot the quantity R eff as a function of mobility reduction in each city with points corresponding to sub-units and colored by time-period starting from the first 100 registered cases. Note that neighboring counties could have a correlation in R eff due to strong mobility and imported cases flows. However, we do observe a neat dispersion in the values measured across counties in the same city ( Fig. 2g-i). Even so, while in Atlanta we see that each region more or less responded independently in terms of reducing transmission, for the case of Chicago we begin to see a pattern emerging, that becomes clear when looking at NYC where the various regions were clustered temporally and reduced transmission in a monotonically decreasing fashion with mobility-reduction. These type of diagrams allows us to obtain a Pearson correlation R p,response from the scattered plots to account for the synchronization between the counties of an urban area. For example, a city like NYC in Fig. 2i has a larger R p,response than Atlanta in Fig. 2g for which R p,response is close to zero.
Generalization to a set of cities. These patterns generalize beyond the three cities considered. In SI Fig.   S4 we plot the equivalent of Fig. 2g-i for the full set of cities, and in Fig. 3a we plot R early measured at city scale in the early stages of the pandemic (three weeks after onset of 100 cases and before the lockdown) as a function of the original . There is a clear connection between transmissibility in terms of R early and the mobility hierarchy, with hierarchical cities showing an increased spread of the disease at the onset. Indeed NYC being the most hierarchical city in the United States, had a 50% higher value of R early than sprawled ones such as Cincinnati, Charlotte or Atlanta.
This increased transmissibility is also reflected in the extent of the spread of the disease as measured by plotting I max versus the original in Fig. 3b. Cities can be separated into those that had already experienced a peak in the incidence curve (Northeastern and Midwestern cities, colored in yellow), and those that were still at the early phases of the pandemic (Southern and Western cities, in red). Restricting to cities where the pandemic was already well-established, we see a clear trend whereby hierarchical cities had a much wider outbreak as compared to sprawled ones. Stronger mitigation strategies manifested in those places where the outbreak was wide-spread 67 , and a relation between mobility reduction and emerges as can be seen in Fig. 3c where we found that NYC, San Francisco and Boston reduced mobility by around 60% , whereas on the other end of the spectrum, Cincinnati, Minneapolis, Atlanta or Charlotte decreased mobility by around 40%.
The connection between mobility reduction and transmission mitigation is also far more pronounced: plotting R p,response as a function of indicates that hierarchical cities seem significantly more responsive to mitigation measures (Fig. 3d). That is, hierarchical cities exhibit more spatially synchronized slowdowns of the transmission after mobility restrictions are applied, while sprawled cities tended to asynchronous evolution at the county level. (Figs. 2 and 3) suggest that hierarchical cities experience faster and more widespread outbreaks as compared to sprawled ones. However, mobility restrictions and lockdowns in those cities were comparatively more effective. To confirm these findings while removing potential confounding factors as differences in population sizes, densities and spatial distributions, variation in type and timing of lockdowns, or noise in the data, we implement a metapopulation SEIR model driven by empirical mobility flows before and after mitigation measures.

Model confirmation of the empirical results. The empirical results
For this, we check if simulating disease spread on multiple idealizations of a city while tuning their mobility structure confirms the empirical observations of the pandemic trends. As representative examples, we consider Atlanta, Chicago and New York and simulate spread while changing by reshuffling the trip destinations. We end up with thirty cities, ten versions for each of the three, each with a different level of hierarchy score. We begin by checking whether the model can reproduce the trend in Fig. 3a in the early stages of the epidemic without lockdowns. The equivalent plot is shown in Fig. 4a, where we extract R early from the incidence curves of the simulations. In contrast to the empirical data where a significant fraction of cases is undetected, we have complete information in our simulations. Consequently, we set the epidemic onset at 10 4 accumulated cases. The Pearson coefficient R 2 for the points of each individual city ranges between 0.95 − 0.97 , indicating a clear signal for the transmissibility increasing with hierarchy. When measuring the correlation for the three cities www.nature.com/scientificreports/ collectively, confounding factors such as differences in population density and mobility distributions emerge, and R 2 decreases to 0.63 in order of magnitude agreement with the one empirically observed in Fig. 3a obtained by aggregating 22 different cities. Similar to Fig. 3b, we reproduce the observed trend of the incidence peak I max being monotonic with respect to as shown in Fig. 4b. This is reflected in the time taken to reach the incidence peak (Fig. 4c) and the final size of the outbreak (Fig. 4d), with more centralized-hierarchical cities experiencing stronger transmission, faster spread and wider prevalence in the population.
Modeling the response of cities to mobility restrictions. On the other hand, Fig. 3d suggests that cities with higher largely achieved a better reduction in transmission with lockdown measures. To simulate this, in Fig. 5a we plot the final size of the pandemic by establishing an 80% reduction in interactions ( X S = 0.8 ) with the lockdown coming into effect at π th = 5 × 10 −3 . Remarkably, we see an inversion of the curve as compared to Fig. 4d with the trends now reversed; the final size of the pandemic is lower in cities with higher and is decreased by an order of magnitude as compared to the scenario with no mitigation. The inversion of the curve resulted from a rather strong lockdown, so in Fig. 5b we show the case for a softer lockdown with X S = 0.4 ; here we see the same trend with as in the case with no lockdown, however, the size of the outbreak is reduced by a factor of five. In addition, the time taken to institute mobility restriction measures is a crucial parameter, and in Fig. 5c, we show the case for an earlier (but still soft) lockdown with X S = 0.4 and π th = 10 −3 , finding a further decrease by a factor of three in the pandemic size. In terms of assessing the effect of different flavors of lockdowns, the connection between the size of the outbreak and its dependence on , is in general a complex function of the epidemiological parameters, the extent of mobility reduction and the distribution of flows. One can get an impression of this connection by plotting in Fig. 5d the final epidemic size as a function of for different lockdown intensities X S . While the curves for X S ≤ 0.4 do not show an inversion but a progressive flattening, lockdowns with X S of 0.6 and above produce the inversion. Although, this occurs for a set of parameters, for a given infectivity β and lockdown threshold π th . www.nature.com/scientificreports/

Robustness of the results.
To confirm the qualitative trends presented thus far, we systematically explore the effects of different parameters on the model results, as well as run our model on a different dataset. A key parameter is the infectivity β and, consequently R 0 . In SI Figs. S8-S10, the epidemic size is studied as a function of for different lockdowns X S while changing β . Strong lockdowns flatten the curves with or even invert it for all the considered cases. The same occurs if we vary the times at which lockdowns are imposed through modifying the parameter π th as shown in SI Figs. S11 and S12. Additionally, the results are also robust to the fraction of people leaving geographical units M as shown in SI Figs. S13 and S14. The mobility data considered here has information on both commuting and non-commuting flows, all of which can serve as potential infection routes. However, one can repeat the analysis using different types of mobility data, and here we do so by using the publicly available commuting data from the United States census bureau's LODES database for the year 2016 collected at the zip code level. In this case, the number of residents, commuting flows and fractions of mobile population can be directly measured. In SI Figs. S15-S18, we use Atlanta as a representative example and plot R eff at early onset, the maximum incidence I max , and the size of the epidemic (for two different flavors of lockdown) as a function of (recomputed at the zip code resolution and only from commuting data) reproducing the trends seen in Fig. 5.
Thus far, we presented a systematic analysis of the parametric landscape of a minimal model that is valid in characterizing the relation between mobility hierarchy and epidemic spreading in cities. In order to focus on the characteristics unique to SARS-CoV-2 we next test the robustness of our findings with a realistic model 63 that adds age diversification of agents, age-dependent contact patterns for the individuals (POLYMOD matrices), an expanded set of health state compartments and a set of parameters fit by the authors of 63 to mimic the real data of daily occupancy in the hospitals of Paris. We employ the same rewiring procedure to scan over multiple values of . The model is then run with data from Paris and from London using both commuting and census data (SI Figs. S19-S23) as well as the Location History data (Supplementary Figs. S24-S26). In all cases, we observe an increase of R early with , as well as an inversion or flattening of the size vs lockdown measures, consistent with the results from the minimal model.
We note that while in our simulations we considered a homogeneous confinement ( X S equal in all S2 cells or zip codes), the mobility reduction as measured directly from data is heterogeneous. Consequently we modify this by assuming a linear relation between mobility reduction by zones and the confinement X S (SI Figs. S27  . Spreading by type of city. Simulations are run in each city with different values of , obtained by the randomization procedure described in the text. Each box reflects 100 runs and displays the median, quartiles, the 5% and 95% confidence intervals. In (a), R early (obtained over three weeks after the onset of 10 4 cases) as a function of , as in Fig. 3a. The peak incidence I max is shown in (b), the time to the peak since the beginning of the simulation, t(I max ) , in (c) and the final epidemic size in (d), all as a function of . All four panels correspond to the baseline mobility before lockdown. www.nature.com/scientificreports/ and S28). Furthermore, we update the values of X S weekly according to data, therefore generating restriction measures that are heterogeneous in both space and time. This refinement contains signatures of socioeconomic diversity of neighborhoods, since residents may work in sectors that were impacted differently by the mitigation measures. All these tests corroborated our main findings: the effective infectiousness of the COVID-19 depends on the urban substrate, in particular, hierarchical cities will have a more explosive epidemic spreading. The appropriate lockdown measures can flatten the size dependence on or even invert it (SI Figs. S29 and S30). Beyond city scale, at the level of counties, the models are also able to reproduce the dependence of the synchronization of the epidemic curve on the hierarchy of cities (Figs. 1g-i and 2d). In SI Fig. S31, we provide the results aggregated at the county scale in New York. As can be observed, while mobility restrictions are kept the same among all NY counties, the way in which hotspots are connected is critical to the synchronization of local outbreaks and mobility. Reproducing the empirical features requires thus hierarchy in the mobility and heterogeneity in the mobility and in the mitigation measures (confinement).
In terms of the robustness of empirical analysis, we show in SI Figs. S31-S37 the variation with temporal period of measurement of the epidemic parameters. One could decide to let more delay in the process of checking the effect of mobility restrictions on R eff . In Fig. 3d, we used a one week delay, in SI Fig. S32 we use a two week delay, finding very similar results. It is true that the value of R eff in early stages can be easily overestimated 68 , hence in Supplementary Fig. S33 we choose to measure the R early as the average value measured in a time period shifted to one week later with respect to Fig. 3a. The results remain unaffected.
The epidemiological data used so far does not report cases by onset of symptoms, but by reporting date. This could introduce some differences 69 . We were not able to find a comprehensive data for the US at county scale with onset of symptoms, hence we rely on the UK dataset (https://coronavirus.data.gov.uk/details/download) which reports cases curves for 210 Upper Tier Local Authorities (UTLA), both for reporting and specimen date. What we see in SI Fig. S34 is that the daily estimation of R eff between cases by specimen date and those by publication date tend to differ slightly. However, in our main analysis we do not compute the R eff day by day, but we take for each weak an average value of R eff for all the seven days. This allows us to keep the R eff estimation in a reasonable range lower then R eff = 3 both for specimen date and publication date cases. This difference among the two estimations of R eff tend to vanish as the weeks go by, showing a higher and higher correlation among all the 210 UTLA. Moreover, in Fig. 3 we consider the average of R eff over the first three weeks of cases to generate R early and this leads to SI Fig. S35, where the fluctuations of the first week are mitigated by the subsequent measurements, Results of a metapopulation model using the S2 cells as basic geographical units. Simulations are run in each city with different values of , obtained by the randomization procedure described in the text. Each box reflects 100 runs and displays the median, quartiles, the 5% and 95% confidence intervals. We explore the pandemic size for different lockdown scenarios: a strong lockdown, X S = 0.8 , for π th = 5 × 10 −3 in (a); a soft lockdown, X S = 0.4 , at the same prevalence in (b); another soft but earlier lockdown X S = 0.4 and π th = 10 −3 in (c); and finally, in (d), a systematic exploration with X S for π th = 5 × 10 −3 in Atlanta. www.nature.com/scientificreports/ leading to a much better estimation, which is highly correlated to the one measured on the specimen date. The same issue regarding the reported cases could affect the measure of the incidence peaks in the various cities. Here, by using the same datasets available in UK we use the same definition of incidence peak of Fig. 3. As shown in SI Fig. S36, the sum of the cases of the week before the peak in the smoothed epidemic curve is pretty consistent among the two definition of cases reporting. Transparency and colormaps are applied proportionally to the local population of the relative UTLA. As can be seen, there are some small deviations among those areas that register low population, in these cases even below 100, 000 inhabitants, hence small fluctuations in the data of low populated areas may strongly affect the measure of the incidence peak with respect to other areas, such as in our case of metropolitan areas. Finally, different lockdowns experienced in cities with diverse hierarchies may depend even on socio-economic variables 70 , such as the amount of work force employed in essential services. In Supplementary Fig. S37 we provide a multivariate analysis with socio-economic factors available at the scale of metropolitan areas (OECD defined areas, http:// www. oecd. org/ cfe/ regio nalde velop ment/ stat. htm). We do not have a specific variable for the fraction of service sector employed workers, however we have the GDP per capita, which seems to represent a good and reliable proxy for this socio-economic variable (see https:// ourwo rldin data. org/ graph er/ gdp-vsservi ces-emplo yment). As shown in SI Fig. S37, the GDP per capita only explains the 28% of the variance in the distribution of the maximum decrease of mobility in cities, versus the 38% explained by the hierarchy ( ). For sake of completeness, we provide a further analysis on the effect of hierarchy on the heterogeneity of mobility reductions in city counties. As we can see from SI Fig. S38, more centralized-hierarchical cities tend to experience more homogeneous mobility reductions among their counties during the lockdown, whereas sprawled cities exhibit higher heterogeneity among mobility reductions, causing a harder containment of the local epidemic. The only outlier is Seattle, which presents a peculiar elongated urban shape, where a mobility reduction in a central county necessarily affects all the others.

Discussion
In summary, we studied how mobility restrictions in urban areas affect the propagation of an infectious disease. We leveraged a dataset that captures aggregate flows of populations around the world, in a consistent way from the period before the pandemic to a few months after the first wave. In a previous work, we had shown that hierarchical cities have better indicators for the use of public transportation, walking, emissions per-capita and health indicators. However, their mobility structure favors spreading of infectious diseases in terms of speed and extent of the contagions. At the same time, the empirical analysis shows that lockdown and travel restriction measures can lead to better outcomes in more hierarchical cities.
Besides empirical observations, we corroborate our findings with a family of metapopulation epidemic models operating at sub-urban scales. The results of the main manuscript are obtained with the simplest version of the models, in which we include the distribution of the population and the mobility information only. This is enough to be able to generate a dependence between the speed and virulence in terms of maximum of incidence and total disease size, and the mobility hierarchy in cities. Such relation originates from the mobility configuration regarding hotspots, where the most diverse contacts concentrate and thanks to which the contagions can spread faster across the city. Hierarchical cities are characterized by intense and better connected hotspots with respect to decentralized ones. Mitigation measures in the form of mobility reductions impact mainly the hotspot and this is why hierarchical cities respond better to lockdowns. This is the main order effect, other factors such as age structure or socio-economic differences across neighborhoods can play a role as well as observed at coarser scales [44][45][46][47] . We include these effects in more elaborate versions of the model and find in all the cases that the relation between hierarchy and spreading speed and the possibility of a better response of hierarchical cities to lockdowns still hold, which means that those factors induce changes of higher order in the propagation patterns. These results have been confirmed using a set of US cities as basis for the simulations, their randomized clones, and also data from London and Paris. Still the relation between mobility hierarchy and spreading is robust and holds too in other cities and even in theoretical settings such as lattices.
Several implications present themselves from a policy making point of view: Given the limited resources for testing and surveillance in much of the developing world (and for that matter in parts of the developed world as well), it seems prudent to classify cities that are more vulnerable to spread ahead of time and deploy those resources there first. Furthermore, the analysis provides clues on the extent to which non-pharmaceutical interventions are effective. It is effective to enforce mitigation measures as early and as thoroughly as possible, given that the time-to-response is particularly crucial in hierarchical cities within or near an epidemic outbreak. Sprawled cities, with their distributed mobility and less connected outbreaks, have a larger time window within which to enforce policy measures, yet even so, the intensity of response is important to reduce the final number of cases. Similar lessons may be employed in terms of vaccination strategies, given its even more limited availability. While well established graduated procedures for vaccinations have been followed worldwide, in the face of exponential growth as has being seen in waves in India and Brazil, it may be more effective to deploy vaccines considering mobility hotspots (e.g., service workers), motivating new research in optimal spatial vaccination strategies. We note that our findings, while presented in the context of COVID-19, are generally applicable to other respiratory infectious diseases. In terms of a post-pandemic urban policies, just as outbreaks of tuberculosis and other diseases shaped the form of the modern urban space in the past, our results also shed new light on how to protect the different type of cities from infectious disease spreading.

Code availability
The code to calculate the flow hierarchy in cities is available in the following link: https:// mygit. katol az. net/ aleix/ flow-hiera rchy.