In the last months, the still ongoing diffusion of the Coronavirus (COVID-19) pandemia has spurred a large body of scientific contributions, attempting to explore how compartmental models [1] [2] [3] [4] can reproduce and predict the spread of the epidemics in different countries and regions [5, 6] .
Most of this work has been focusing on models in which the mixing process between people in different states or compartments does not depend on the social or geographical space where they are embedded in. However, some previous literature has shown that the (complex) structure of networks describing the way agents can meet, and possibly get infected, may affect the dynamics of the epidemic diffusion and its long-run properties [7] [8] [9] [10] [11] [12] [13] [14] [15] . Furthermore, as the virus spreads in the network, the structure of interactions between people may change over time, due to quarantining measures and/or spatialdistancing policies, which may possibly introduce a coevolutionary effect dynamically linking disease diffusion and network properties [16] [17] [18] .
Motivated by these observations, the paper introduces a generalized spatial SEIRD model that, besides the standard four compartments (susceptible, exposed, infected, recovered, dead), also considers an additional 'quarantined' state, i.e. a SEIQRD model [19] . I explore how the properties of the spread of the epidemics depend on: (i) the structure of the social/geographic network initially connecting the agents in the population, which matches infective and susceptible agents; (ii) the evolution of the share of quarantined and recovered agents (as well as so-cial distancing policies), which dynamically destroy or re-establish social links.
More specifically, I play with a finite population of agents (i.e. nodes) initially placed on four different families of interaction structures: (a) regular 2-dimensional lattices with Moore neighborhoods; (b) small-world lattice [20] ; (c) Erdös-Renyi random graphs [21] ; (d) scalefree (preferential-attachment) networks [22] . I then investigate via Monte-Carlo simulations how the epidemic diffusion is affected by network structures, as their initial average degree increases (which in turn makes their topological properties change) and as the coupled dynamics of quarantined and recovered people deletes and restores social interaction links. Finally, I examine how alternative spatial-distancing policies, which are taking again center stage in the political and social debate as the second wave of COVID-19 rolls across Europe and elsewhere, interacts with the coevolutionary process of disease diffusion and network updating.
I begin describing a simple model where no spatial distancing policies are enforced. Consider a population P of N agents living in a city, which is initially isolated from other cities. Time is discrete and, for the only sake of convenience, I will use the terms 'time periods' or 'days' as synonyms. Agents physically interact according to a simple, undirected, binary graph without self-loops, which at time t = 0 is defined as G 0 = (P, L 0 ) , where P = {1, . . . , N } and L 0 is the initial edge list, defined as the set of pairs (i, j) such that i = j, i ∈ P , j ∈ P , and (i, j) ∈ L 0 if and only if there exists an edge between i and j at t = 0. The graph G 0 -which, as we will see below, is going to evolve through time as the epidemics spreads-can be considered as describing social or geographical links through which people normally meet friends or neighbors.
At time t = 0 all nodes are in the state S (susceptible), but a randomly-chosen share θ of them becomes exposed (i.e., θN agents become in state E, due to a random inflow of infective agents from other cities). People in state E enter in an incubation period without symptoms and are not infectious. At any t > 0, I assume that each agent i ∈ P meets all its neighbors, i.e. all j ∈ V it , where V it = {j ∈ P : (i, j) ∈ L t } and L t is the current edge list. In each time period, transitions between compartments (i.e., states) occur through a parallel updating mechanism according to the following rules:
(a) An agent in state E becomes in state I (infective) after an incubation of D time periods, where D is an i.i.d random variable with probability distribution p(D).
(b) An agent in state E becomes infected with probability π = 1 − (1 − α) k if h/she meets k infective agents in its neighborhood, where α is a parameter tuning the likelihood of becoming infected in a single direct meeting and 0 ≤ k ≤ |V it |. 1 (c) An agent in state I becomes quarantined (in state Q) with a daily quarantine rate (DQR t ). Agents in state Q cannot meet anyone, i.e. they instantaneously cut all their bilateral links with their neighbors. 2 (d) An agent in state Q dies (i.e., becomes in state D) with a daily death rate (DRR t ), recovers (in state R) with a daily recovery rate (DRR t ) or stays quarantined otherwise. Recovered agents are assumed to be immunized and re-establish connections that they used to have in G 0 (provided that neighbors are still alive and are not quarantined).
A flow-chart description of model dynamics is provided in the Supplementary Material (SM), see Figure S1 .
The initial network G 0 is assumed to belong to one out of the following graph families:
(i) Regular 2-dimensional boundary-less lattices endowed with the Chebyshev distance (LA henceforth). This defines squared Moore neighborhoods of radius r LA ≥ 1 and degrees k LA i = (2r LA +1) 2 −1 for all i.
(ii) Small-worlds lattice [20] built starting from nodes placed on a ring, with rewiring probability p SW > 0 and expected average degreek SW = 2r SW , where r SW ≥ 1 is the interaction radius on the initial ring (SW henceforth).
(iii) Erdös-Renyi random graphs [21] , with link probability p ER > 0 and expected average degreek ER = (N − 1)p ER (ER henceforth).
(iv) Scale-free networks with linear preferentialattachment [22] and entrance of m SF ≥ 1 new nodes, generating an expected average degreē
To summarize network topology, I focus, besides average degree, on three statistics that have been found to influence, in general, the spread of epidemics on graphs [10] . These are: the standard deviation of node degree distribution (s k ), global clustering coefficient (c) and average path-length ( ), computed ignoring infinite path-lengths between nodes of different components. Their expected values (with standard errors) are reported in SM, Table  S1 . To get a better feel, fixing r LA ∈ {1, 2, 3, 4}, and thus k ∈ {8, 24, 48, 80}, s k , c and approximately scale ask β , with β > 0 for s k and c and β < 0 for in all networks.
All simulations refer to a population of N = 1024 agents (chosen to build a square lattice with edge L = 32) and a number of days T sufficient to reach a steady state.
The epidemic parameters of the model are calibrated using data at the national level for Italy, made available by "Dipartimento della Protezione Civile", see https://github.com/pcm-dpc/COVID-19/tree/ master/dati-andamento-nazionale, covering the period from February, 22nd onward. In the simulations, I assume for simplicity that DRR t = DRR, DQR t = DQR and DDR t = DDR and, on the basis of empirical diffusion curves, I build three epidemic scenarios:
(i) strong-impact scenario: (DQR, DDR, DRR) = (0.20, 0.10, 0.10); (ii) mid-impact scenario: (DQR, DDR, DRR) = (0.15, 0.07, 0.15); (ii) low-impact scenario:
(DQR, DDR, DRR) = (0.10, 0.04, 0.20) -see SM, Section S2 for more details.
Since the theoretical infection probability in a single meeting cannot be directly observed, I play with values of α that are (0.20, 0.10, 0.05), respectively in the three scenarios. The percentage θ of exposed agents in day 0 is set to 5% throughout. The probability distribution p(D) of incubation days is calibrated using results from [23] , who show that D is log-normally distributed with parameters (µ, σ) = (1.621, 0.418).
As to initial network structures, I experiment with average degreesk ∈ {8, 24, 48, 80}. These values result from setting r LA ∈ {1, 2, 3, 4}. Therefore, it follows that r LA ∈ {4, 12, 24, 40}, p ER ∈ {k(N − 1) −1 ,k = 8, 24, 48, 80} and m SF ∈ {4, 12, 24, 40}. See SM, Table S2 , for a summary of parameter setups.
For each choice of model parameters, I independently run M = 1000 simulations. This Monte Carlo sample size is sufficient to get standard errors for across-simulation averages small enough to ensure that differences between averages are always statistically significant.
In order to get insights about within-simulation model behavior, I keep track of several within-simulation statistics, i.e. computed in each day of the epidemic diffusion. These include: population shares in each compartment, death and cure rates, the share of agents who become infected through meetings, and the four network metrics k, s k , c, and -which change across time as the result of the evolution population shares in each compartment. Another statistics of interest is the populationaverage of the number of neighbors that each I agent has infected daily (ρ henceforth; see SM, Section S3 for details), which can be employed as a rough estimate of the basic reproduction number (R 0 ) of the epidemics. Finally, I will also look at the spatial correlation coefficients of compartments (SCCC), calculating, for each state {S, E, I, Q, R, D} the fraction of all existing edges in the network whose endpoints end up being in the same state.
To summarize the aggregate behavior of the model (i.e. across runs), the following set of additional statistics are computed: (i) peak-time of infections (PTI), defined as the first day in which the share of infected people reach its overall maximum; (ii) the shares of agents in states {S, I, R, D} at the end of simulation (EoS) and at PTI; (iii) the sum over all compartments of SCCC at PTI; (iv) the EoS share of agents who become infected through meetings; (v) the values of network metricsk, s k , c and at PTI. Furthermore, I provide an estimate of the first day after whichρ goes below one (cf. SM, Section S3).
Monte Carlo averages of all the above summarizing statistics will then be compared across initial networks families, initial average degrees, and epidemiological setups.
I begin studying the dynamic behavior of disease spreading across the four network families, focusing on the 'Mid Impact' epidemic scenario withk = 8 (see Figure 1) . Irrespective of the initial network structure, the population converges to a similar share of deaths, but in ER and SF networks a small percentage of S people still remains (see below, Section III C).
In these two networks, epidemic diffusion reaches a higher peak of infections than in the case of LA and SW , since more agents become exposed a little earlier. This is because in ER and SF networks the average number of infections per agent grows very quickly during the outbreak of the epidemic process, and then decreases earlier and more sharply than in LA and SW networks (SM, Figure S3 ). The evolution of SCCC shows, indeed, that the shares of edges linking two E or two I agents cross near to PTI and display a more abrupt inverse-U-shaped pattern over time, illustrating how the virus spreads across neighborhoods (SM, Figure S4 ).
As the epidemic process develops over time, the share of Q agents first grows and then declines. This impacts on the network structure, because quarantined agents become isolated, constraining in turn the diffusion of the disease. Furthermore, the more the infection weakens, the more quarantined people recover and re-establish some of their initial connections.
To get a better feel about this co-evolutionary process, Figure 2 shows how network metrics, normalized to match the [0, 1] interval, change during a simulation. Bothk and decrease towards their minimum value across time in LA and SW , with a pace slowing down as R people spread in the population. The decline of is due to the growing number of small connected components and isolated nodes created by Q and D agents. In LA and SW networks, however, recovered agents that reestablish their connections are able to slightly boost average degree and reconnect isolated clusters. More marked differences across network structures emerge when looking at s k and c. In LA, SW and, particularly, in ER graphs s k first increases due to the injection of Q agents, then as R and D gradually replace Q patients, it oscillates until getting to a stable level. In SF graphs, instead, s k follows the same time pattern ofk and , as initial heterogeneity is very high and cannot be further increased by the interplay between Q, R and D shares. Therefore, populations where the epidemics diffuse in LA, SW and ER networks end up having a higher final heterogeneity of degrees, while the opposite holds for SF graphs. The final clustering level, instead, is almost completely recovered, but with opposite patterns. LA and SW networks first experience an increase in c (albeit very moderate in magnitude) because the diffusion evolves less quickly. Instead, in ER and SF graphs, some triads are rapidly destroyed by Q people and then R people re-establish them when the epidemics softens.
In the SM, Section S4, I also show that, as the share of Q agents first increases and then decreases, and that of R agents keeps growing in time, the topological properties of the network changes in very heterogeneous ways, depending on the family to which it belongs. This is due to the dynamic removal and re-establishment of links -which affects in non-trivial ways, in particular, the standard deviation of node degrees and global clustering coefficients-and ultimately impacts on the properties of the diffusion process itself.
B. The Impact of Initial Average Degree I now investigate the behavior of the model when the initial average degree increases in the range {8, 24, 40, 80}, keeping fixed the epidemic scenario to the 'Mid Impact' one. If agents initially have, on average, more neighbors they can meet more infective people. Therefore, the probability to become E increases for the population in each single day. However, a largerk does not imply the at the end of the simulation (EoS) there will be a larger fraction of deaths and/or recovered, as this is mainly affected by the epidemic parameters. What changes in the speed at which the contagion evolves and some of its dynamic properties.
For example, as shown in Figure 3 both the peak-time of infections and the estimate of the first day after which ρ goes below one, quickly decrease withk. Furthermore, as the initial average degree grows, the contagion evolves more quickly in ER and, especially, in SF networks.
Furthermore, in all networks, the fraction of infected people at PTI immediately jumps up whenk increases fromz 8 to 24, and then keeps growing withk but less quickly (cf. Figure S5 in the SM). This implies that, since the epidemic scenario is fixed, the share of agents that are quarantined in the first days of the contagion increases more than linearly. Therefore, at PTI, the shares of susceptible, recovered and dead agents actually decrease with initial average degree. 
Next, I explore what happens in the model when alternative epidemic scenarios are assumed (cf. Section II C). For the sake of comparison, I keep fixedk = 8 throughout. Simulation results show that, as expected, EoS shares of dead (respectively, recovered) agents decrease (respectively, increase) in all network setups as one moves from the bad to the good epidemic scenario (see Figure 4 ). More interestingly, within the same scenario, the model behaves differently across network setups, and these differences are amplified as the contagion is less strong. Indeed, in ER and SF networks the epidemics diffuses quicker than in LA and SW graphs -as documented in SM, Figure S6 , panels (a) and (c). Therefore, LA and SW display more S (and less I) agents at PTI than ER and SW do -see Figure 4 , panel (c)-and a significantly smaller spatial correlation of compartments (panel (b) in SM, Figure S6 ). At the end of the simulation (EoS), conversely, many more susceptible agents remain in ER and, especially, in SF networks. This is due to the higher heterogeneity of the degree distribution in such networks: the existence of many small-degree nodes at the beginning of the process prevents them to be infected, especially when the contagion becomes softer and their few neighbors are quickly quarantined. As a consequence, slightly smaller shares of deaths are observed in ER and, in particular, in SF networks at EoS.
Spatial distancing (SD) is implemented in the model in a very stylized way. I assume that the city government only tracks the evolution of Q agents and enforces SD when x t (Q) > q , where x t (Q) is the current share of agents in the Q compartment and q ∈ (0, 1). The SD policy aims at making more difficult face-to-face meetings between neighbors, and can be enforced with increasing strengths. Of course, its ex-post effectiveness also depends on how strictly people follow the rules. Here, I do not separately model the ex-ante plans of the government and the response of the agents. Therefore, more formally, I define θ ∈ (0, 1) as the ex-post effectiveness of SD policy and assume that, under SD, an agent meets each neighbor in any time period t with probability ψ = 1 − θ. This implies that, under SD, an agent in state E now becomes infected with probability:
where k is the number of infective agents the agent meets in its neighborhood. I allow for two versions of SD: (i) permanent, if SD is enforced from the first day when Figure 5 plots EoS shares of agents under SD (either permanent or temporary) minus the correspondent share without SD. In each SD setup, I target the share of people ending up in either S or D compartments, and the share of deaths (D). Results show that, as expected, a permanent SD policy is better than a temporary one independently of network structure. However, especially when a strong setup is enforced in the permanent SD policy version, networked populations that benefit the most are those where agents are located on either lattices or small-worlds. Conversely, temporary SD policies are more effective in ER and, in particular, in SF networks, provided that they are implemented more rigorously. This is due to how network structures evolve during a typical run, see Figure 2 . Indeed, when a permanent SD policy is likely to be implemented, LA and SW exhibits larger average degrees and clustering than ER and SF. This prevents the infection to be transmitted more effectively during the peak. Instead, enforcing temporary SD policies allow an even smaller probability that low-degree agents remain susceptible, which is more likely to happen in ER and SF networks, due to their higher degree variability. When such a policy is switched off, ER and SF systems display higher (and more dispersed) average degrees and larger clustering than in the LA and SW cases, but the share of infected people is now smaller. Therefore, one observes less deaths. Disaggregating S and R shares shows also that, in the permanent SD case, the improvement in LA and SW is obtained via an almost similar increase of both compartments. On the contrary, when SD is temporary, much of the improvement is due to an increase in EoS susceptible agents only.
In this paper, I studied a generalized spatial SEIRD model to explore the impact of alternative social-network structures on the diffusion of the COVID-19 disease. The introduction of quarantined agents generates a coevolving process between epidemic spreading and network structure, ultimately shaping steady-state outcomes and the speed of diffusion.
In the simplest framework, without spatial distancing policies and a given benchmark choice of initial average degrees and epidemic parameters, the initial network structure does not affect the final shares of susceptible, dead and recovered people, but it strongly impact on the timing and the speed of diffusion. In ER and, in particular, in SF networks, more agents become exposed earlier and diffusion takes place quicker and more strongly than in the LA and SW cases. This is linked with how network structure coevolves across time with the shares of Q, R and D agents. Indeed, in ER and SF networks, average degree initially decreases less sharply than it does in LA and SW. Furthermore, degree variation and clustering is higher. Therefore, the probability of becoming exposed increases, as susceptible agents face larger and more clustered neighborhoods. Increasing initial average degree, while keeping fixed epidemic parameters, thus results in a faster speed of infection, especially in ER and SF networks, both in terms of smaller PTIs and average number of neighbors that each agent has infected daily. When instead different epidemic scenarios are assumed for a fixed initial degree, network structure impacts differently model behavior, and these differences are amplified as the strength of the contagion weakens. In particular, since the epidemics initially diffuses quicker in ER and SF networks, one typically observes more S (and less I) agents at PTI in LA and SW graphs, and many more remaining S agents at EoS in ER and SF networks (with slightly smaller shares of deaths). The effect of SD policies depends in the model on the strength with which they are enforced, as well as whether they are temporary of permanent. In particular, whereas permanent SD policies allow for better results than temporary ones irrespective of network structure, permanent (and strong) SD measures are more effective in LA and SW structures, whereas temporary (and strong) SD policies should be preferred if interactions occur through ER or SF graphs. This is again due to the interplay between network structure and compartment shares in the evolution of the epidemics. Indeed, switching on and off SD policies may hit the system when the topological properties of its network structure are very different, depending on the initial graph family describing social interactions.
More generally, results suggest that, in order to predict how epidemic phenomena evolve in networked populations, it is not enough to focus on the properties of initial interaction structures. In fact, if the epidemic diffusion requires quarantining people, and possibly enforcing SD policies, the coevolution of network structures and compartment shares strongly shape the way in which the virus spreads into the population, especially in terms of its speed. On the one hand, the average and standard deviation of degree distribution, as well as clustering, of initial networks are, together with epidemic parameters, important determinants of the subsequent diffusion patterns. On the other hand, the topology of social interaction structures evolves over time, due the rise and fall of Q, R and D agents, in different and non trivial ways across alternative network families, and this in turn impacts diffusion patterns. As a result, the timing and features of SD policies may dramatically influence their effectiveness.
The foregoing analysis can be extended and improved in several directions. To begin with, alternative parametrizations for the epidemic process, more in line with evidence from the ongoing second wave, could be tested. Furthermore, it would be interesting to assess the extent to which results are robust to increasing population size, additional network structures (e.g., coreperiphery graphs), and different values for the share of agents that become initially exposed. In this last respect, one could also play with alternative assumptions as to the mechanism governing the way in which exposures initially occur, e.g. allowing for the emergence of spatiallyclustered exposed agents, instead of just supposing that a randomly-chosen share of people get infected. Finally, one can perform a deeper analysis to better understand how the topology of network structures influences epidemic diffusion, for example asking whether centrality indicators such as k-coreness measures [24, 25] can help in investigating the role of super spreaders [26] . Figure S1 . Flow-chart of model dynamics in any given time period. Table S1 . 
In order to calibrate epidemic parameters for the Italian case, I use data from "Dipartimento della Protezione Civile", available at https://github.com/pcm-dpc/COVID-19/tree/master/dati-andamento-nazionale, covering the period from February, 22nd onward. The time series of total detected cases (T C) is broken into: (i) currently positive total detected cases, i.e. individuals showing severe symptoms that are therefore quarantined (Q), either at home or in the hospital; (ii) recovered patients (R); (iii) dead patients (D).
Daily quarantined, recovered, and death rates are computed as ratios to total detected cases (T C), see Figure  S2 . Since this figure is highly dependent on the number of swabs taken and asymptomatic individuals, the daily quarantined rate (DQR), in particular, is likely to suffer from strong overestimation. Indeed, several contributions have shown that for Italy actual infected people may be 5-10 times higher [27] . Therefore, in the figure we also report the DQR for these overestimation scenarios. We do not rescale DRR and DDR accordingly as we are here interested in the fraction of people who, after being quarantined, transition in the R and Q compartments.
Simulation scenarios in Table S2 are then built using the range of DQR 5 , DQR 1 0,DRR and DDR.
In each day, agent i in state S may possibly become E if h/she meets k ≥ 1 neighbors j ∈ {j 1 , . . . , j k } ⊂ P in state I. If this happens, the burden of having infected i is divided in equal shares η j h = 1/k across all infecting neighbors. If an I agent contributes over the day to the infection of more than one neighbor, all its η shares are accordingly cumulated. At the end of each day, after all infections have been eventually occurred, we computeρ as the population average (over non-zero elements) of η shares.
Theρ statistics can be interpreted as a rough estimate of the "basic reproduction number" (R 0 ) of the epidemic process, as it measures the average number of cases that each infected agent generates.
Given a single simulation, the time-seriesρ t typically goes down with t, as the set of infective agents shrinks and the number of recovered (or dead) people increases. It is therefore worthwhile to estimate t * = min{1 ≤ t ≤ T :ρ t < 1}, which is computed at the end of each simulation as one of its summarizing statistics. Since in each simulationρ t may oscillate many times around one before decreasing persistently below one, I fit the series with four different functions: (i) polynomial of degree 1: f (t) = at + b; (ii) polynomial of degree 2: f (t) = at 2 + bt + c; (ii) power: f (t) = at b ; (iii) exponential f (t) = a exp (−bt); and then for the best fit f * (according to the associated adjusted R 2 ), I define t * = min{1 ≤ t ≤ T : f * (t) < 1} for t = 1, . . . , T .
This Section reports more evidence about how the evolution of Q and R shares co-evolve, within a simulation, with the structure of the network where agents are embedded in. I still focus on a 'Mid Impact' epidemic scenario, with k = 8 for convenience, and plot in Figure S7 the within-simulation time series of network metrics, re-scaled to match the unit interval, against the shares of Q and R agents. Average path length ( ) is not taken into consideration here as its relationship withk is monotone along the process.
A quick inspection of Figure S7 suggests the process undergoes a series of phases that depend on the family of network considered. These phases are summarized in the following tables, where for each relationship between network metrics (k, s k , c), shares of agents in Q and R states, and network structure, I identify phases (i.e., subsets of timeseries evolution) where they display a particular co-movement. For example, a pair (↑, ↓) for the entry (x, y) means a phase where the network metrics x increases and the share of agent in state y decreases. 