key: cord-0150972-tlbica5n
authors: Ventura, Paulo Cesar; Aleta, Alberto; Rodrigues, Francisco Aparecido; Moreno, Yamir
title: Modeling the effects of social distancing on the large-scale spreading of diseases
date: 2021-05-20
journal: nan
DOI: nan
sha: 5496b2a90f4f106eaf8c9af044aef23bcea53c46
doc_id: 150972
cord_uid: tlbica5n

To contain the propagation of emerging diseases that are transmissible from human to human, non-pharmaceutical interventions (NPIs) aimed at reducing the interactions between humans are usually implemented. One example of the latter kind of measures is social distancing, which can be either policy-driven or can arise endogenously in the population as a consequence of the fear of infection. However, if NPIs are lifted before the population reaches herd immunity, further re-introductions of the pathogen would lead to secondary infections. Here we study the effects of different social distancing schemes on the large scale spreading of diseases. Specifically, we generalize metapopulation models to include social distancing mechanisms at the subpopulation level and model short- and long-term strategies that are fed with local or global information about the epidemics. We show that different model ingredients might lead to very diverse outcomes in different subpopulations. Our results suggest that there is not a unique answer to the question of whether contention measures are more efficient if implemented and managed locally or globally and that model outcomes depends on how the full complexity of human interactions is taken into account.

The spreading of infectious diseases is a complex process involving two main aspects. On the one hand, the spreading capabilities of the pathogen depend on its biological properties, the characteristics of the host, and the many environmental factors that can play an important role. On the other hand, the pathogen (except for vector-borne diseases) can only spread if two hosts have some kind of contact. As such, the behavior of the host is a key element in the study of epidemic spreading [1, 2] . For this reason, many researches have focused on studying how animals and humans interact, in order to inform mathematical models and produce better forecasts [3] [4] [5] [6] .

Despite having acknowledged the role that behavior plays in the spreading of infectious diseases, epidemic models usually neglect the possibility that hosts will change their behavior due to an ongoing outbreak [7] . An exception to this are awareness models, in which both an epidemic and information spread at the same time in the population and host reacts accordingly (see [8] and the references therein). These behavioral changes have been observed both in animal [9] and human societies [10] . More recently, the spreading of COVID-19 has clearly demonstrated the variety of ways in which humans react to an epidemic. For instance, without forceful government intervention, traffic in Seoul's subway declined sharply following the first deaths in South Korea [11] . Conversely, several communities in other regions of the world defied social distancing measures [12, 13] .

Public health authorities base many decisions on the forecasts produced by epidemic models. It is thus of paramount importance to include the effect of behavioral changes as something inherently attached to the spreading of the epidemic [7] . Of particular interest in this regard are metapopulation models. These models provide a simple way of incorporating the spatial heterogeneity of human societies, while keeping them mathematically tractable. In essence, in a metapopulation model hosts are grouped in different subpopulations, and the exchange of individuals between subpopulations is governed by certain rules. Within each subpopulation it is possible to assume that the spreading from host to host follows the classical homogeneous mixing approach [14] [15] [16] , or include more details on individual heterogeneities through different techniques [17, 18] .

These models have been extensively used during the emergence of new pathogens to study how an outbreak might propagate globally [19] . In the particular case of COVID-19, they have also been used to study the early spreading of the disease across countries [20] , but also within countries [21] [22] [23] . Yet, these models usually focus on the effect of varying the rules governing the flux of individuals between subpopulations, or on globally modifying the transmission due to the introduction of public health interventions. In this paper, our objective is to understand the role that social distancing can play in these models and shed some light on the effect that it arXiv:2105.09697v1 [physics.soc-ph] 20 can have on the spatial spreading of epidemics when the response is heterogeneous across regions.

We implement a discrete and stochastic SIRmetapopulation model composed by V subpopulations [24] [25] [26] . In each subpopulation, individuals can interact with each other and spread the disease following the classical SIR model. These subpopulations are connected through a certain network, so that an individual may travel to another subpopulation only if it exists a direct link between the source and target subpopulations.

To include social distancing effect, either policy-driven or as a consequence of fear of infection, we consider an additional coefficient that modifies the transmissibility of individuals. This coefficient mimics the behavioral responses of the population, and evolves according to the densities of infectious and recovered individuals [7] . In section III, we study this model in an homogeneously mixed population. For the remaining sections, we apply the model to a metapopulation system.

For the SIR compartmental model, individuals are assigned compartments according to their infectious status: susceptible (S) if they do not have the disease and can catch it; infectious (I) when they have the disease and can transmit it to susceptibles, and removed (R) when they no longer transmit the disease after being infectious (either by recovery or deceasing). At each model's epidemic update, the following rules determine the transitions between compartments:

• S → I: a susceptible individual in subpopulation i can become infectious with probability P i (S →

where β is the individual transmission probability, N i is the number of individuals in subpopulation i and I i is the number of infectious individuals in the same subpopulation. The term a i is called the coefficient of contact reduction, and depends on the considered strategy as we describe later on.

• I → R: an infectious individual is moved to the R (removed) compartment with probability µ, which is the inverse of the average infectious period. After this event, the individual no longer participates on the epidemic dynamics.

Throughout this paper, we set µ to 1/4 and the basic reproductive number R 0 = β/µ to 1.5, which is compatible with the parameters of an influenza-like disease. We define N S = V i=1 S i as the total number of susceptible individuals in the whole population and, analogously, N I and N R for the infectious and removed compartments. The coefficient of contact reduction a i is what determines the behavioral responses to the epidemics, and we consider different scenarios for such a response, each one with a different definition of a i :

This scenario is based on the one proposed by Eksin et al. in [7] , which here we extend to metapopulations. The coefficient of contact reduction emulates the social response to an increase in the number of infections, and it is a function of the total (global) number of infected and recovered individuals, given by:

Since a i only depends on global quantities, its value is the same for all subpopulations. The parameter k is the response strength, an adjustable exponent that calibrates the intensity of the response against the number of cases. The coefficient l is called coefficient of long-term, and determines the importance of the R compartment to the behavioral response. Following [7] , l = 0 represents the short-term strategy (ST), where the response weakens when the incidence drops. The case l = 1 is the long-term strategy (LT), in which the awareness is proportional to the total prevalence. These two limiting cases differ in many aspects, which we describe throughout the rest of the paper. 

In this variant of the previous scenario, we consider that each subpopulation (node) responds individually, according to its own number of cases. The expression for the coefficient of contact reduction is:

For comparison, we also implement a more traditional scenario in which the transmissibility is reduced by a constant factor a 0 once the overall number of cases N I + N R overcomes a given threshold N · ε I (with 0 ≤ ε I < 1). Here we use only the global information, thus a i is defined as:

and is the same for all subpopulations. This scenario can be considered as a baseline in which social distancing (or governments policies) is constant and independent of the state of the system, and it is regarded here as a benchmark that is often assumed as the behavioral response in metapopulation epidemics. [23, [27] [28] [29] .

Following a standard metapopulation framework, the mobility between subpopulations is modeled as a random walk through the links of a graph of V nodes and M links. It is controlled by a master parameter τ called the mobility coefficient, as well as the weights T ij of existing links between subpopulations (nodes) i and j. During a mobility update, each individual in subpopulation i travels to subpopulation j with probability p ij = τ · T ij /N i (0), where N i (0) is the number of individuals attributed to subpopulation i at the initial time step. Within this scheme, we have the following features:

• The average number of individuals expected to travel from i to j is τ · T ij .

• If all links are reciprocal and symmetric (i.e., T ji = T ij for all connected pairs of subpopulations i, j), then the net fluxes between the subpopulations are balanced and the population of each one remains approximately constant, fluctuating around N i (0). We use this configuration throughout the paper.

• The probability that an individual in subpopulation i travels anywhere is p i = τ j T ij /N i (0), and may vary for each subpopulation. We always set τ to be small enough so that p i < 1 for every subpopulation i.

During a single time step, we first perform the epidemic interactions in each subpopulation, then update the number of individuals that are in each state. After this, we apply the mobility rules to determine how many individuals move through each link, and update the actual numbers only after all fluxes have been calculated. This way, the results do not depend on the order at which we "visit" the subpopulations to perform the calculations. The mobility and epidemic models are schematically represented in Figure 1 .

For sufficiently low mobility between subpopulations, the local dynamics can be well described by an isolated homogeneously mixed system. Also for sufficiently high number of individuals, we can use rate equations for the expected fractions of the population in each compartment, substituting S/N , I/N and R/N by ρ S , ρ I and ρ R , respectively. The chosen expression for the behavioral response mechanism is simple enough to allow for analytical manipulation. Particularly, in the long-term (LT) strategy (i.e., with l = 1 in Eqs. 1 and 2), these rate equations can be integrated to find the trajectory of the dynamical system. In this section, we explore this tractability to extract some insights from the model.

Considering the dynamics of a single isolated subpopulation, the dynamical equations for the average fractions ρ S , ρ I and ρ R of susceptible, infectious and removed individuals can be written as:

where R 0 = β/µ is the basic reproduction number and µ is now interpreted as the recovery rate. For the longterm strategy, a(ρ I , ρ R ) = (1 − (ρ I + ρ R )) k = ρ k S , given that ρ S + ρ I + ρ R = 1. As described in [7] , we can divide equation 4 by equation 6 to obtain a separable differential equation:

For k > 0 and assuming that ρ R (0) = ρ 0 ≥ 0 and ρ I (0) = δ → 0, its solution is given by:

Alternatively, writing ρ I as a function of ρ R :

represents the trajectory of the system in a ρ I vs ρ R diagram, showing how the fraction of infectious individuals peaks as the prevalence increases in the LT strategy. Also, the final ρ R prevalence at the end of the outbreak can be found by solving equation 9 for ρ I = 0 besides the trivial solution ρ R = ρ 0 . The solid lines in Figure 2 show such trajectories for different values of the response strength k, ρ 0 = 0 and R 0 = 1.5.

For k = 0, equation 7 solves as a simple SIR model which, for ρ R (0) = ρ 0 and ρ I (0) = δ → 0, leads to a ρ I vs ρ R trajectory given by:

In the long-term strategy with k > 0, the contact reduction coefficient a(ρ I , ρ R ) is a decreasing function of ρ R + ρ I , which in turn can only increase over time. This means that a(ρ I , ρ R ) will also always increase as if the disease contention measures are held forever. This is not reasonable to assume for the long term behavior of a real situation and, for this reason, we can consider that at some point after the main outbreak, a(ρ I , ρ R ) is set to 1 again, which we call henceforth a system release (representing a release of the contention measures).

For instance, consider that a system release occurs when the fraction of infectious individuals reaches a small value δ > 0 after the main outbreak initiated. This comprises the fact that, during the time evolution given by equations 4 to 6, ρ I (t) tends to but never actually reaches zero. For the calculation of the ρ I vs ρ R trajectories though, we can consider δ → 0. After the system is released, the trajectory follows Eq. 10 with ρ 0 set to the final size of the first outbreak. In Figure 2 , the dashed lines represent the secondary outbreaks produced after the system's release.

As observed in the figure, a system release can generate a second outbreak which, for higher values of k, can be greater than the first one. This happens because, although the long-term strategy effectively reduces the size of the first outbreak (both in ρ R and ρ I ), it may leave the system under its herd immunity threshold, and thus vulnerable to new outbreaks [30] . In this homogeneous model, the threshold for herd immunity is given by ρ herd = 1 − R −1 0 , meaning that for ρ R ≥ ρ herd the fraction of infectious individuals can no longer increase. If the disease dies out before the prevalence reaches that value, the system will experience a new outbreak if the social distancing measures are relaxed and the disease is reintroduced in the population. Note that the critical value of k for reaching herd immunity at the first outbreak is k = 1, for which the final outbreak size is exactly solved as ρ * R = 1 − R −1 0 , as noticed by Eksin and others [7] . 

Besides releasing the system after the first outbreak, we propose another way to work with multiple outbreaks in the LT strategy. What makes the a(ρ I , ρ R ) coefficient to be by default strictly decreasing is its dependence on ρ R , which holds a "long-term memory". We can modify and possibly reset such memory by subtracting a constant ρ 0 from the prevalence that's considered for the contact reduction coefficient a(ρ I , ρ R ). If we do this right after the main outbreak and set ρ 0 to the recovered fraction at that time, then a(ρ I , ρ R ) is momentarily reset to 1 (i.e., no contention measures), but the system will still react in the event of another outbreak. We call this a system reset (the memory of the population is reset, but the social distancing measures still apply), in contrast with the simpler system release described in the previous section (in which the contention measures are completely removed). We are again assuming that ρ I is arbitrarily small at the end of an outbreak.

During this new round of the dynamics, equation 4 divided by equation 6 yields the following differential equation:

which is still separable, but for ρ 0 > 0 solves into a less insightful expression:

where we define:

where 2 F 1 is the Gaussian hypergeometric function. Equations 12 and 13 provide ρ R for any given ρ S . For each of these values, ρ I is calculated as ρ I = 1 − ρ S − ρ R , which finally allows the construction of a ρ I vs ρ R trajectory. Figure 3 shows the trajectories of the model with system resets, for different values of k. The first outbreak of each execution is represented by solid lines, and are traced using equation 9. For the executions that did not achieve the herd immunity threshold (which is shown as a black dot-dashed line), subsequent outbreaks are represented by dashed lines and traced using equations 12 and 13. In this scenario, stricter contention measures (that is, higher k) always cause smaller infectious peaks, but can generate more secondary outbreaks, making it more difficult to comply with the measures. The inset of Figure 3 shows in more detail the secondary outbreaks for k = 2, 4 and 8, where we can see that there are multiple outbreaks with decreasing peak size. For the short-term (ST) strategy, the coefficient of contact reduction depends only on ρ I , being written as a(ρ I ) = (1 − ρ I ) k = (ρ R + ρ S ) k . This breaks the separability of the differential equation obtained by dividing equations 4 and 6, leaving no simple method to solve it for arbitrary values of k. We can still use a classic Runge-Kutta of order 4(5) [31] to integrate the equations in time, and plot the ρ I vs ρ R trajectories.

FIG. 4. Numerical trajectories of the homogeneously mixed short-term (ST) strategy, obtained from Runge-Kutta integration of equations 4 to 6 starting with ρI (0) = 1 · 10 −5 and ρR(0) = 0. R0 is set to 1.5. Unlike the LT strategy, the system always reaches the herd immunity condition in a single outbreak, without considering any modification to the baseline model. Figure 4 shows the numerically solved trajectories of the system for different values of k. The main difference with respect to the LT strategy is that the system eventually reaches the herd immunity condition within the first and only outbreak. Higher values of k reduce the size of the peak in ρ I , at the expense of larger time-span in which the contention measures have to be sustained. This cannot be visualized in Figure 4 as the time parameter is implicit, but we further explore the interplay between peak size and time span in section V A. Note also that even for very large values of k the final fraction of removed individuals is beyond the herd immunity threshold.

The previous analysis has shown that, even if local extinction might be achievable with social distancing, it renders the system vulnerable to further reintroductions of the pathogen. Thus, if mobility between subpopula-tions is allowed, there could be spill overs from those with still ongoing outbreaks to those which contained the epidemic. To study this type of events, we now focus on the proposed model in a metapopulation using Monte Carlo simulations, with the algorithm described in section II. For the metapopulation, we use a random geometric network (RGN) of V = 50 subpopulations, constructed in a square space of length 1 and connecting subpopulations that are closer than d = 0.25. At this initial point, all links are unweighted and reciprocal. This gives an expected average degree of V πd 2 ≈ 9.82, though the specific realization we used through this paper has an average degree of 7.2. We choose such a network topology because it emulates the spatial distribution of cities in small to mid scale, where size hierarchy and long-range links are not very present.

Once this unweighted graph is constructed, we set the initial population of each subpopulation i proportional to its degree k i , according to N i (0) = N pop k i /(2M ) , where we set N pop = 10 7 . This makes the overall population to be not exactly but very close to N pop , only deviated due to truncation. Then we set the weights of existing links between subpopulations i and j as T ij = N i (0)N j (0)/N pop = T ji . Within this scheme, the local population sizes fluctuate around N (0) over time, as explained before. Moreover, the most connected subpopulations are also most populous ones, though the RGN is reasonably homogeneous in this sense. For consistency of the results, we also use a fixed subpopulation as the seed of the disease, seeding 10 infectious individuals at the beginning of the simulations, with all other subpopulations starting with susceptibles only. The seeded subpopulation was chosen to be around the center of the square space.

For the long-term (LT) strategy with global information, we can compare the outcomes of the simulations with those of a constant response after a threshold (explained in section II A). Figure 5 shows the outbreak size (total number of infected individuals) in each subpopulation for the two strategies, each one averaged over 1000 independent executions. The constant coefficient of contact reduction a 0 = 0.78 was chosen to yield a global outbreak size of ∼ 0.20, approximately equal to that obtained with the LT global strategy with k = 2. We notice that the LT strategy provides more heterogeneous outcomes between subpopulations compared to the constant response. Specially, as it can be seen in Figure 6 , subpopulations that are farther from the seed tend to have smaller outbreak sizes. This happens because, in the global LT strategy, the intensity of contention measures due to public awareness is the same for all subpopulations and increases with time, thus subpopulations that are seeded later have smaller effective reproduction numbers. The fact that our proposed model for behavioral responses with global strategy produces more heterogeneous outbreak sizes is notorious, as this is generally observed after outbreaks of real diseases. In particular, this situation is compatible with the one observed in 2020 during the COVID-19 pandemic in those countries that imposed strict global lockdowns even if only part of the country was severely affected [20, 32, 33] . In contrast, a uniform strategy would lead to more homogeneous outcomes.

We now describe a framework to compare local and global strategies in terms of their efficiency, characterized by the costs and benefits of each strategy. Direct comparison of simulations with global and local strategies using the same response strength k is inappropriate, as global strategies require greater values of k to yield similar effects. We therefore define two metrics, one for the cost and another for the effectiveness, and compare local and global strategies in parametric plots of such metrics (with k as an implicit parameter). As long-term and short-term strategies are qualitatively different, we apply different metrics to characterize each one.

The short-term strategy, in which the contact reduction coefficient only responds to the (either local or global) density of infectious individuals, is characterized by a slow progression of the system towards its herd immunity. A higher value of the response strength k represents a more effective response, which reflects into a smaller prevalence peak (i.e., the maximum ρ I ) yet longer outbreak duration. The outbreak size (i.e., maximum of ρ R ) however does not vary much with k, as it essentially depends on the herd immunity limit of the system. Therefore, the infectious peak size ρ max I is a reasonable measure of effectiveness in this case, with smaller peaks attributed to more effective strategies.

Higher values of k also imply more time in which contention measures are active to control disease propagation. This translates into a i being smaller than 1 during longer times. We can quantify "how much and for how long" the contention measures are applied with the quantity:

which accounts for the intensity and time span of the contact reduction. This provides, in a broad sense, a metric for the cost of the short term strategy, considering that contention measures typically bring costs to the population. Notice that these metrics must be calculated for each subpopulation.

For the same population and RGN structure described in section IV, we compare local and global strategies with respect to ρ max I and J for simulations with different values of k. Figure 7 shows the plots of these metrics against each other, grouped by sets of subpopulations according to their distance to the seeded subpopulation. The metric values are averaged over 1000 independent executions.

From the plots, it is clear that local strategies outperform global ones for all subpopulations, as the former produces much smaller peaks sizes for similar values of J in the considered range. This means that, for the ST strategy, the use of local information about contagions is more efficient than using a unified, global strategy. This happens because using global information may not be well suited for the epidemic situation of a given subpopulation. For instance, subpopulations that are far from the seed effectively apply the contention measures before the epidemic arrives, but as these measures are relaxed at some moment (when the global prevalence decreases), these subpopulations will still undergo an intense outbreak, given that they were only subject to minor outbreaks and the majority of the individuals in these subpopulations are susceptible.

Although we only display the results for a single value of the mobility coefficient τ = 10 −3 , the same conclusions are obtained for τ = 10 −2 and 10 −4 , but the advantage of local strategies over global ones is more pronounced with lower mobility rates. This is expected, because with lower τ values, the dynamics of each subpopulation is less coupled by mobility, making local strategies more practical.

In the long-term strategy, the overall intensity of contention measures increases with time (though it can locally decrease due to migration), as they are proportional to both I and R densities. As shown in section III for homogeneous mixing, for sufficiently high k, the epidemic spreading is halted before herd immunity is achieved, leaving the system vulnerable to secondary waves in case of a system release without other immunization policies. We can still characterize the effectiveness of the strategy in this first wave, either by its outbreak size (ρ max R ) or by the peak size (ρ max I ). For consistency with section V A, we chose to work with ρ max I as well. For the cost of the strategy, J (given by equation 14) does not provide a reliable measure, as in this case a i does not approach 1 again at the end of the main outbreak, making J overly sensitive to the time taken until the epidemic vanishes. We choose instead to simply work with a max = 1 − min a(t), which is the maximum level of contact reduction adopted by the subpopulation during the simulation. This metric disregards the temporal evolution of a i , but still provides a number that is proportional to the level of contention measures adopted by the subpopulations without the sensibility issue of J.

Using the same RGN setup, we find that the cost vs effectiveness curves are more complex for the LT strategies than for the ST ones. Figure 8 shows the ρ max R x a max parametric curves for local and global strategies, and for three values of the mobility master parameter τ = 10 −2 , 10 −3 and 10 −4 . In this case, there is not a clear advantage of local strategies over global ones, and this depends both on mobility levels and the distance to the seed. For instance, when τ = 10 −2 (Figure 8.a) ), a local strategy is better for the seed, while it is generally worse for more distant subpopulations. For immediate neighbors of the seed (one step away), the curves cross each other, making the optimal strategy to depend on the value of k.

For lower mobility values (Figure 8 .b) and c)), another interesting feature is evident: the curves for the local strategies are not monotonic with the cost, meaning that two strategies (given by two different values of k) may have the same cost but notably different effectiveness. This happens because the LT strategy can affect the invasion threshold and attack rate of the system, preventing some subpopulations from being reached by the disease for sufficiently high k. In this case, as the strategy is based on local prevalences, these subpopulations will not have to implement contention measures, which explains the decrease in the strategy cost.

We adapt the metapopulation model with a long-term global strategy to consider system resets, as explained in section III B. This allows us to analyze secondary outbreaks caused by a relief in the contention measures, represented here as a reset in the "memory" of the contact reduction mechanism (that is, the value of ρ R in the a i coefficient).

The system resets are implemented as follows: the contention measures are activated when the overall fraction of infectious individuals ρ I surpasses an activation threshold ε in I = 10 −4 (activation threshold), and then deactivated (possibly after an outbreak) if ρ I goes under another threshold ε out I = 0.8 · 10 −4 < ε ∈ I . At each activation event, the "memory" of the system is reset, that is, the coefficient of contact reduction is given by a i (t) = (1 − (ρ I (t) + ρ R (t) − ρ 0 )) k , with ρ 0 set to the value of ρ R at the moment of the activation event. This mechanism is slightly different from that used in section III B for the rate equations model, but the distinction ε out I < ε in I is needed in stochastic simulations so that small fluctuations around ε I do not increase the number of outbreaks. To control the number of secondary outbreaks, we also limit the number of system reset events to 10. If this number of reset events is reached before the extinction of the epidemics, the system is then released, that is, a i is permanently set to 1.

In Figure 9 , we show the time series of the fraction of infectious individuals of each subpopulation (colored thin curves), as well as the global fraction of infectious individuals (black shaded curve) and the value of a i over time (red dashed curve) for a typical execution of the model. Using the same RGN population as in previous sections, we set the response strength as k = 5 and the mobility coefficient to τ = 10 −3 . Each vertical ascent of the a i value represents the occurrence of a system reset. We notice that, at each reset, a new set of local outbreaks occurs, raising again the overall incidence and leading the system to readopt contention measures. Differently from the homogeneously mixed case (as in Figure 3 ), the secondary outbreaks can be greater in peak size than the first one. This feature is an essential difference between homogeneously mixed populations and structured ones like a metapopulation, and occurs because some (usually many) subpopulations are not reached by the primary outbreak.

We can better understand this feature by looking at the local incidence. Figure 10 shows a network map of the greatest peak size (given by ρ max dividuals for some selected subpopulations (lateral panels). The curve of panel (a) corresponds to the seeded subpopulation, displaying a single and pronounced peak in the early stage of the process. Other subpopulations may also present a single main outbreak (as in panel (e)), two (separate or close) outbreaks ((d) and (c)) or even more outbreaks ((b) and (f)). This spatial heterogeneity driven by contention measures has been observed in several locations during the COVID-19 pandemic [21, [34] [35] [36] .

The features of a simulation with system resets are better presented through a single, typical execution of the simulation. We show that the observed pattern of multiple outbreaks is a solid feature of our model by counting the number of outbreaks of each time series, then plotting histograms for 1000 independent executions. We use a simple algorithm to determine the number of outbreaks of each local time series, described as follows: every time the fraction of infectious individuals crosses a given threshold φ in = 1 · 10 −3 from bellow and, posteriorly, another threshold φ out = 0.5φ in < φ in from above, an outbreak is accounted, and the time interval between these two crossings is regarded as a single outbreak. The difference between the two thresholds reduces spurious detection of outbreaks due to stochastic fluctuations, though this may still marginally occur. Figure 11 shows the statistics of accounted local outbreaks. On panel (a), the number of outbreaks for each subpopulation and each of the 1000 executions is put into a single histogram, showing that most subpopulations present a single outbreak, but often higher number of outbreaks also occur. On panel (b), we take the average number of outbreaks over all subpopulations in a single execution, then plot a histogram for the executions. From it we see that the typical average number of outbreaks is between 1 and 1.4, meaning that most executions present at least some subpopulations that undergo multiple outbreaks.

We have presented a model that incorporates dynamical behavioral responses to epidemics, which could be driven by governmental policies or by the endogenous response of individuals (e.g. fear of infection), in the context of metapopulations. The emerging features of the model are very rich, showing a complex landscape of outcomes depending on the implemented strategy.

First, we have shown that for isolated populations, strong social distancing measures (LT) are able to ablate the epidemic, but render the system vulnerable. Indeed, as long as the prevalence is below the herd immunity threshold, a reintroduction of the pathogen after the measures are lifted will inevitably lead to a second wave. On the other hand, soft social distancing measures (ST) contain the exponential growth of the epidemic and leave the system in a state in which the prevalence is over the herd immunity threshold, preventing further outbreaks. Note, however, that for a severe disease in which the infection fatality rate is substantial, larger prevalence implies a larger number of deaths. Thus, stronger social distancing policies will reduce the number of deaths at the expense of leaving the system vulnerable, while softer measures will control the spread but increase the number of deaths. These observations are particularly important in the context of metapopulations. Since subpopulations are not isolated, in the absence of additional control measures, the disease might be seeded again in disease-free areas by individuals who where infected in other regions. Thus, besides long and short term strategies, we also need to incorporate whether the contention measures use global or local information.

We compared the LT and ST strategies of our model, in an attempt to address the question of whether the contention measures are more efficient if implemented and managed locally (at each subpopulation) or globally (equally to the whole population). For the ST strategy, we measured its cost by the J metric that quantifies the intensity and duration of the contention measures, and used the maximum number of simultaneous infections (i.e., the peak size) to quantify the strategy's performance. We showed that the local strategy always outperformed the global one, which is valid for any value of the mobility parameter. For the LT strategy, we quantified the cost by the maximum intensity reached by the contention measures, and the performance by the peak size. In this case, the cost/benefit relationship was more complex, depending on the mobility rate and the distance to the seeded subpopulation. Global strategies are generally better for distant subpopulations (with respect to the seed). This situation was precisely the one that emerged in the first wave of infections in 2020 during the COVID-19 pandemic in countries that imposed countrywide lockdowns even if only one region was severely af- fected [20, 32, 33 ]. Yet, this is reversed if the mobility is low enough so that the local strategy can stop the epidemics before these nodes are reached. For the seeded subpopulation and its neighbors, local strategies are typically preferred. For this type of strategy, therefore, the choice between global and local strategies is not trivial and should be addressed appropriately.

However, it is important to notice that long-term strategies assume a permanent adoption of contention measures, which leaves the system vulnerable to secondary outbreaks if these measures are lifted. We address this by implementing memory resets which, after outbreaks, instantly lifts the intensity of contention measures but leaves its mechanism active. Due to the residual presence of infectious individuals, the system faces secondary outbreaks. Thus, subpopulations that experienced a milder first wave, might have worse outcomes in subsequent waves. Besides, regardless of the contention measures, once they are lifted the system keeps progressing towards the herd immunity threshold, but each population might experience waves of different intensity at different times. The spatial heterogeneity driven by contention measures depicted in this paper is compatible with the one observed during the COVID-19 pandemic in several regions of the world [21, [34] [35] [36] .

In summary, we have shown than even in a simple metapopulation model a very complex scenario emerges. Depending on the mobility and distance to the seed, local vs global strategies yield different results and introduce different heterogeneities. As such, the full complexity of human gatherings and behaviors should be accounted for to effectively deal with emerging diseases. 

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

Close encounters of the infectious kind: methods to measure social mixing behaviour

High-resolution contact networks of free-ranging domestic dogs canis familiaris and implications for transmission of infection

Social Contacts and Mixing Patterns Relevant to the Spread of Infectious Diseases

Projecting social contact matrices to different demographic structures

A high-resolution human contact network for infectious disease transmission

Systematic biases in disease forecasting-the role of behavior change

Epidemic spreading with awareness and different timescales in multiplex networks

Social network plasticity decreases disease transmission in a eusocial insect

The Public's Response to the 2009 H1N1 Influenza Pandemic

Changes in Subway Ridership in Response to COVID-19 in Seoul

Israel's response to the COVID-19 pandemic: tailoring measures for vulnerable cultural minority populations

Partisan differences in physical distancing are linked to health outcomes during the COVID-19 pandemic

Spatial Heterogeneity in Epidemic Models

Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations

Human mobility networks and persistence of rapidly mutating pathogens

Comparing large-scale computational approaches to epidemic modeling: Agent-based versus structured metapopulation models

Metapopulation epidemic models with heterogeneous mixing and travel behaviour

Human Mobility Networks, Travel Restrictions, and the Global Spread of 2009 H1N1 Pandemic

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

Ferreira. Outbreak diversity in epidemic waves propagating through distinct geographical scales

A data-driven assessment of early travel restrictions related to the spreading of the novel COVID-19 within mainland China

Evaluation of the potential incidence of COVID-19 and effectiveness of containment measures in Spain: a data-driven approach

Grenfell. Metapopulation Dynamics of Infectious Diseases

Seven challenges for metapopulation models of epidemics, including households models

Spatial epidemiology of networked metapopulation: an overview

Towards a Characterization of Behavior-Disease Models

Metapopulation Network Models for Understanding, Predicting, and Managing the Coronavirus Disease COVID-19

Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases

Data-driven estimate of SARS-CoV-2 herd immunity threshold in populations with individual contact pattern variations. medRxiv

A family of embedded runge-kutta formulae

Villafáfila Gomila, Lluis Carbo Saladrigas, Adoración Hurtado Fernández

Nieves Ascunce Elizaga, María Ederra Sanz, Carmen Ezpeleta Baquedano, Ana Bustinduy Bascaran, Susana Iglesias Tamayo, Luis Elorduy Otazua, Rebeca Benarroch Benarroch, Jesús Lopera Flores, and Antonia Vázquez de la Villa. Prevalence of SARS-CoV-2 in Spain (ENE-COVID): a nationwide, population-based seroepidemiological study

Stefano Merler, Patrizio Pezzotti, and the COVID-19 working group. Epidemiological characteristics of COVID-19 cases and estimates of the reproductive numbers 1 month into the epidemic

A spatial analysis of the COVID-19 period prevalence in U.S. counties through June 28, 2020: where geography matters?

An interactive web-based dashboard to track COVID-19 in real time

Impact of the accuracy of case-based surveillance data on the estimation of time-varying reproduction numbers. medRxiv