]Av. Trabalhador Sancarlense, 400, 13566590, São Carlos - SP, Brazil

# Modeling the effects of social distancing on the large-scale spreading of
diseases

Paulo Cesar Ventura1<EMAIL_ADDRESS>[  Alberto Aleta2 Francisco A.
Rodrigues3 Yamir Moreno2,4 1 Instituto de Física de São Carlos, Universidade
de São Paulo, São Carlos, SP, Brazil.
2 ISI Foundation, Turin, Italy
3 Instituto de Ciências Matemáticas e de Computação, Universidade de São
Paulo, São Carlos, SP, Brazil.
4Institute for Biocomputation and Physics of Complex Systems (BIFI), 50018
Zaragoza, Spain
5Department of Theoretical Physics, University of Zaragoza, 50018 Zaragoza,
Spain

###### Abstract

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.

††preprint: APS/123-QED

##  I Introduction

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 Funk2010Sep ; Read2012Dec . For this reason,
many researches have focused on studying how animals and humans interact, in
order to inform mathematical models and produce better forecasts Wilson2019 ;
Mossong2008Mar ; Arregui2018Dec ; Salathe2010Dec .

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
eksin2019systematic . 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 daSilva2019Sep  and the references therein). These
behavioral changes have been observed both in animal Stroeymeyt2018Nov  and
human societies SteelFisher2010Jun . 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
Park2020Apr . Conversely, several communities in other regions of the world
defied social distancing measures Waitzberg2020Dec ; Gollwitzer2020Nov .

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 eksin2019systematic . 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 Lloyd1996Mar ; Colizza2008Apr ; Aleta2017Mar , or
include more details on individual heterogeneities through different
techniques Ajelli2010Dec ; Apolloni2014Dec .

These models have been extensively used during the emergence of new pathogens
to study how an outbreak might propagate globally Bajardi2011Jan . In the
particular case of COVID-19, they have also been used to study the early
spreading of the disease across countries Chinazzi2020Apr , but also within
countries Costa2020Dec ; Aleta2020Oct ; Aleta2020Dec . 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 can have on the spatial
spreading of epidemics when the response is heterogeneous across regions.

##  II Model description

We implement a discrete and stochastic SIR-metapopulation model composed by
$V$ subpopulations Keeling2004Jan ; Ball2015Mar ; Wang2014Oct . 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 eksin2019systematic . In section III, we
study this model in an homogeneously mixed population. For the remaining
sections, we apply the model to a metapopulation system.

Description | Symbol
---|---
Individual infection probability | $\beta$
Recovery probability/rate | $\mu$
Basic reproduction number | $R_{0}=\beta/\mu$
Contact reduction coefficient | $a_{i}{}$
Response strength exponent | $k$
Long-term coefficient | $l$
Mobility coefficient | $\tau$
Threshold for infectious fraction | $\varepsilon_{I}$
Individual counts in each state in subpopulation $i$ | $S_{i},I_{i},R_{i}$
Global individual counts in each state | $N_{S},N_{I},N_{R}$
Global fractions in each state | $\rho_{S},\rho_{I},\rho_{R}$
Number of nodes (subpopulations) on network | $V{}$
Number of links on network | $M$
Number of individuals in subpopulation $i$ | $N_{i}$
Total number of individuals in population | $N_{pop}$
Maximum distance for connection on RGN | $d$
Link weight between subpopulations $i$ and $j$ | $T_{ij}$
Long-term strategy | LT
Short-term strategy | ST
network | RGN
Table 1: List of symbols and acronyms.

###  II.1 Epidemic spreading

Figure 1: Scheme of the epidemic and mobility models in a metapopulation. The
SIR epidemic dynamics occur inside each subpopulation, where homogeneous
mixing is assumed. Also at each time step, individuals move between
neighboring subpopulations $i$ and $j$ according to a mobility matrix
$T_{ij}$.

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 $\rightarrow$ I: a susceptible individual in subpopulation $i$ can become
infectious with probability $P_{i}(S\rightarrow
I)=1-(1-a_{i}\beta/N_{i})^{I_{i}}$, where $\beta$ 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 $\rightarrow$ R: an infectious individual is moved to the R (removed)
compartment with probability $\mu$, 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 $\mu$ to $1/4$ and the basic reproductive number
$R_{0}=\beta/\mu$ to $1.5$, which is compatible with the parameters of an
influenza-like disease. We define $N_{S}=\sum_{i=1}^{V}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}$:

  1. 1.

Social distancing based on global information

This scenario is based on the one proposed by Eksin _et al._ in
eksin2019systematic , 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:

$a_{i}=\left(1-\frac{(N_{I}+l\cdot N_{R})}{N}\right)^{k}$ (1)

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 eksin2019systematic , $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.

  2. 2.

Social distancing based on local information

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:

$a_{i}=\left(1-\frac{(I_{i}+l\cdot R_{i})}{N_{i}}\right)^{k}$ (2)

  3. 3.

Constant response after threshold

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\cdot\varepsilon_{I}{}$ (with $0\leq\varepsilon_{I}{}<1$). Here we use only
the global information, thus $a_{i}{}$ is defined as:

$a_{i}=\begin{cases}1&N_{I}+N_{R}<N\varepsilon_{I}{}\\\ a_{0}&N_{I}+N_{R}\geq
N\varepsilon_{I}{}\end{cases}$ (3)

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.
Perra2011Aug ; Aleta2020Dec ; Calvetti2020Jun ; Manfredi2013 .

###  II.2 Mobility model

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 $\tau$
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}=\tau\cdot 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
$\tau\cdot 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}=\tau\sum_{j}T_{ij}/N_{i}(0)$, and may vary for each subpopulation. We
always set $\tau$ 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.

##  III Analytical insights for homogeneously mixed populations

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 $\rho_{S}$, $\rho_{I}$ and $\rho_{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.

###  III.1 Long-term strategy

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

$\displaystyle\dot{\rho_{S}}$ $\displaystyle=-\mu
R_{0}\,\rho_{S}\rho_{I}\cdot\,a(\rho_{I},\rho_{R})$ (4)
$\displaystyle\dot{\rho_{I}}$ $\displaystyle=\mu
R_{0}\,\rho_{S}\rho_{I}\cdot\,a(\rho_{I},\rho_{R})-\mu\rho_{I}$ (5)
$\displaystyle\dot{\rho_{R}}$ $\displaystyle=\mu\rho_{I},$ (6)

where $R_{0}=\beta/\mu$ is the basic reproduction number and $\mu$ is now
interpreted as the recovery rate. For the long-term strategy,
$a(\rho_{I},\rho_{R})=(1-(\rho_{I}+\rho_{R}))^{k}=\rho_{S}^{k}$, given that
$\rho_{S}+\rho_{I}+\rho_{R}=1$. As described in eksin2019systematic , we can
divide equation 4 by equation 6 to obtain a separable differential equation:

$\frac{\text{d}\rho_{S}}{\text{d}\rho_{R}}=-R_{0}\rho_{S}^{k+1}$ (7)

For $k>0$ and assuming that $\rho_{R}(0)=\rho_{0}\geq 0$ and
$\rho_{I}(0)=\delta\rightarrow 0$, its solution is given by:

$\rho_{S}=\left[\frac{1}{(1-\rho_{0})^{k}}+kR_{0}(\rho_{R}-\rho_{0})\right]^{-\frac{1}{k}}$
(8)

Alternatively, writing $\rho_{I}$ as a function of $\rho_{R}$:

$\rho_{I}=1-\rho_{R}-\left[\frac{1}{(1-\rho_{0})^{k}}+kR_{0}(\rho_{R}-\rho_{0})\right]^{-\frac{1}{k}}$
(9)

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

For $k=0$, equation 7 solves as a simple SIR model which, for
$\rho_{R}(0)=\rho_{0}$ and $\rho_{I}(0)=\delta\rightarrow 0$, leads to a
$\rho_{I}$ _vs_ $\rho_{R}$ trajectory given by:

$\rho_{I}=1-\rho_{R}-(1-\rho_{0})e^{-R_{0}(\rho_{R}-\rho_{0})}$ (10)

In the long-term strategy with $k>0$, the contact reduction coefficient
$a(\rho_{I},\rho_{R})$ is a decreasing function of $\rho_{R}+\rho_{I}$, which
in turn can only increase over time. This means that $a(\rho_{I},\rho_{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(\rho_{I},\rho_{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 $\delta>0$ after the main
outbreak initiated. This comprises the fact that, during the time evolution
given by equations 4 to 6, $\rho_{I}(t)$ tends to but never actually reaches
zero. For the calculation of the $\rho_{I}$ _vs_ $\rho_{R}$ trajectories
though, we can consider $\delta\rightarrow 0$. After the system is released,
the trajectory follows Eq. 10 with $\rho_{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.

Figure 2: Analytical trajectories of the homogeneously mixed long-term (LT)
strategy with system release after the first outbreak (i.e., complete lift of
the contention measures). The trajectories are combinations of equations 9
(solid lines) and 10 (dashed lines) with $\rho_{0}$ set to the size of the
first outbreak. $R_{0}$ is set to $1.5$.

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 $\rho_{R}$ and $\rho_{I}$), it may leave the
system under its herd immunity threshold, and thus vulnerable to new outbreaks
Lu2021Mar . In this homogeneous model, the threshold for herd immunity is
given by $\rho_{herd}=1-R_{0}^{-1}$, meaning that for
$\rho_{R}\geq\rho_{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
$\rho_{R}^{*}=1-R_{0}^{-1}$, as noticed by Eksin and others
eksin2019systematic .

###  III.2 System resets

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(\rho_{I},\rho_{R})$ coefficient to be by default strictly decreasing is its
dependence on $\rho_{R}$, which holds a “long-term memory”. We can modify and
possibly reset such memory by subtracting a constant $\rho_{0}$ from the
prevalence that’s considered for the contact reduction coefficient
$a(\rho_{I},\rho_{R})$. If we do this right after the main outbreak and set
$\rho_{0}$ to the recovered fraction at that time, then $a(\rho_{I},\rho_{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 $\rho_{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:

$\frac{\text{d}\rho_{S}}{\text{d}\rho_{R}}=-R_{0}\cdot(\rho_{S}-\rho_{0})^{k}\rho_{S}$
(11)

which is still separable, but for $\rho_{0}>0$ solves into a less insightful
expression:

$R_{0}\cdot(\rho_{R}-\rho_{0})=P_{\rho_{0},k}(\rho_{S})$ (12)

where we define:

$\displaystyle P_{\rho_{0},k}(\rho_{S})$
$\displaystyle=\int_{\rho_{S}}^{1-\rho_{0}}\frac{\text{d}u}{u(u+\rho_{0})^{k}}=$
(13)
$\displaystyle=\left.-\frac{1}{k(u+\rho_{0})^{k}}\left(1+\frac{\rho_{0}}{u}\right)^{k}{}_{2}{F_{1}}\left(k,k;k+1,-\frac{\rho_{0}}{u}\right)\right|_{\rho_{S}}^{1-\rho_{0}}$

where ${}_{2}F_{1}$ is the Gaussian hypergeometric function. Equations 12 and
13 provide $\rho_{R}$ for any given $\rho_{S}$. For each of these values,
$\rho_{I}$ is calculated as $\rho_{I}=1-\rho_{S}-\rho_{R}$, which finally
allows the construction of a $\rho_{I}$ _vs_ $\rho_{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.

Figure 3: Analytical trajectories of the homogeneously mixed long-term (LT)
strategy with system resets. The first outbreak (solid lines) follows equation
9, while subsequent outbreaks (dashed lines) are traced using equations 12 and
13. The inset shows a sub-region of the main plot (indicated by a blue dotted
rectangle) where secondary outbreaks are more easily visible, showing how
there may be multiple and successively smaller peaks. $R_{0}$ is set to 1.5.

###  III.3 Short-term strategy

For the short-term (ST) strategy, the coefficient of contact reduction depends
only on $\rho_{I}$, being written as
$a(\rho_{I})=(1-\rho_{I})^{k}=(\rho_{R}+\rho_{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) dormand1980family  to integrate
the equations in time, and plot the $\rho_{I}$ _vs_ $\rho_{R}$ trajectories.

Figure 4: Numerical trajectories of the homogeneously mixed short-term (ST)
strategy, obtained from Runge-Kutta integration of equations 4 to 6 starting
with $\rho_{I}(0)=1\cdot 10^{-5}$ and $\rho_{R}(0)=0$. $R_{0}$ 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 $\rho_{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.1. Note also that even for very large values of $k$ the final
fraction of removed individuals is beyond the herd immunity threshold.

##  IV Global strategies and the heterogeneity of local features

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 subpopulations 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\pi d^{2}\approx 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)=\lfloor N_{pop}k_{i}/(2M)\rfloor$, 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.

###  IV.1 Long-term (LT) global strategy

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.1). 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 $\sim 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.

Figure 5: Comparison of the average outbreak size in each subpopulation
between a) global LT strategy and b) global constant response with threshold.
For the LT strategy, a response strength $k=2$ was used, while the constant
response uses $a_{0}=0.78$ and a threshold of $\varepsilon_{I}{}=10^{-3}$,
meaning that the contact reduction is triggered after $10^{4}$ overall
infections (out of $N_{pop}{}=10^{7}$ individuals) are registered. Both
simulations use $\tau=1\cdot 10^{-3}$ as the travel parameter. Figure 6: Local
outbreak size (given by the maximum of $R_{i}/N_{i}$ over time) in each
subpopulation as a function of its shortest path length to the seed, for
global LT strategy with response strength $k=2$ (circles) and for constant
response with factor $a_{0}=0.78$ (squares). The y-axis values represent
averages over $1000$ independent executions. The boxes show the inner
quartiles of each set. The mobility coefficient is set to $\tau=1\cdot
10^{-3}$.

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 Pollan2020Aug ; Riccardo2020Dec ;
Chinazzi2020Apr . In contrast, a uniform strategy would lead to more
homogeneous outcomes.

##  V The efficiency of local and global strategies

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.

###  V.1 Short-term (ST) strategy

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
$\rho_{I}$) yet longer outbreak duration. The outbreak size (i.e., maximum of
$\rho_{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
$\rho_{I}^{\text{max}}$ 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:

$J{}_{i}=\sum_{t=0}^{\infty}(1-a_{i}(t))$ (14)

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 $\rho_{I}^{\text{max}}$ 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.

Figure 7: Parametric curves for the effectiveness and cost of local and global
ST strategy, for values of $k$ ranging from $3$ to $90$. The first plot is for
the seeded subpopulation, whereas other plots are arithmetic averages over
sets of subpopulations according to their shortest path length to the seeded
subpopulation. The last plot is an arithmetic average over all subpopulations.
The arrows indicate the direction towards which $k$ increases. The mobility
factor is set to $\tau=1\cdot 10^{-3}$. $\rho_{I}^{\text{max}}$ is the maximum
over time of the local infection prevalence $\rho_{I}$, while $J$ is given by
equation 14.

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 $\tau=10^{-3}$, the same conclusions are obtained for
$\tau=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 $\tau$ values, the dynamics of each subpopulation is less
coupled by mobility, making local strategies more practical.

###  V.2 Long-term (LT) strategy

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 ($\rho_{R}^{\text{max}}$) or by the
peak size ($\rho_{I}^{\text{max}}$). For consistency with section V.1, we
chose to work with $\rho_{I}^{\text{max}}$ 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
$\overline{a}_{\text{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
$\rho_{R}^{\text{max}}$ x $\overline{a}_{\text{max}}$ parametric curves for
local and global strategies, and for three values of the mobility master
parameter $\tau=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
$\tau=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.

Figure 8: Parametric curves for the effectiveness and cost of local and global
LT strategy, for the mobility parameter (a) $\tau=10^{-2}$, (b) $\tau=10^{-3}$
and (c) $\tau=10^{-4}$. The $k$ values range from $0.25$ to $50$. The leftmost
panels represent the seeded subpopulation, the middle panels represent an
arithmetic average over its immediate neighbors, and the rightmost panels show
the averages over all subpopulations. The arrows indicate the direction
towards which $k$ increases. The cost $\overline{a}_{\text{max}}$ is the
maximum value of $1-a_{i}{}$ reached during the simulation.

##  VI Secondary outbreaks in the global long-term strategy

We adapt the metapopulation model with a long-term global strategy to consider
system resets, as explained in section III.2. 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 $\rho_{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 $\rho_{I}$
surpasses an activation threshold $\varepsilon_{I}^{\text{in}}=10^{-4}$
(activation threshold), and then deactivated (possibly after an outbreak) if
$\rho_{I}$ goes under another threshold $\varepsilon_{I}^{\text{out}}=0.8\cdot
10^{-4}<\varepsilon_{I}^{\in}$. 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-(\rho_{I}(t)+\rho_{R}(t)-\rho_{0}))^{k}$, with $\rho_{0}$ set
to the value of $\rho_{R}$ at the moment of the activation event. This
mechanism is slightly different from that used in section III.2 for the rate
equations model, but the distinction
$\varepsilon_{I}^{\text{out}}<\varepsilon_{I}^{\text{in}}$ is needed in
stochastic simulations so that small fluctuations around $\varepsilon_{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$.

Figure 9: Time series of the fraction of infectious individuals in each
subpopulation (colored lines), as well as in the whole population (black
line). The global value of $a(t)$ is also shown as a red dot-dashed line (with
the upper limit $a(t)=1$ as a red thin horizontal line). The parameters
(described in text) are set to: $k=5$, $\varepsilon_{I}^{\text{in}}=10^{-4}$,
$\varepsilon_{I}^{\text{out}}=0.8\cdot 10^{-4}$, $\tau=10^{-3}$.  Figure 10:
Map of the final outbreak sizes ($\rho_{R}^{\text{max}}$) of each
subpopulation, for a single execution with the memory reset mechanism (main
panel). The panels show the time series of $I_{i}/N_{i}$ for some sample
subpopulations, showing that the mechanism introduces diversity in the
epidemic trajectory over the population. The average number of detected
outbreaks for this run was $1.20$. Figure 11: Histograms of (a) the number of
local outbreaks (at subpopulation level) for all subpopulations, and (b) of
the average number of outbreaks in each execution (the average is taken over
subpopulations, but not over executions). The parameters of the simulations
are set to: $k=5$, $\varepsilon_{I}^{\text{in}}=1\cdot 10^{-4}$ and
$\varepsilon_{I}^{\text{out}}=0.8\cdot 10^{-4}$. The number of system reset
events was constrained to at most 10.

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 $\tau=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
$\rho_{I}^{\text{max}}{}$) of each subpopulation, along with the time series
of the fraction of infectious individuals 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
Costa2020Dec ; Sun2020Dec ; Dong2020May ; Starnini2020Jun .

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
$\phi_{\text{in}}=1\cdot 10^{-3}$ from bellow and, posteriorly, another
threshold $\phi_{\text{out}}=0.5\phi_{\text{in}}<\phi_{\text{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.

##  VII Conclusions

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
country-wide lockdowns even if only one region was severely affected
Pollan2020Aug ; Riccardo2020Dec ; Chinazzi2020Apr . 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 Costa2020Dec ; Sun2020Dec ;
Dong2020May ; Starnini2020Jun .

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.

##  VIII Acknowledgments

P.C.V. acknowledges the financial support of FAPESP through grants
2016/24555-0 and 2019/11183-5. Y.M. acknowledges partial support from the
Government of Aragón, Spain through grant E36-20R, by MINECO and FEDER funds
(grant FIS2017-87519-P). A.A. and Y.M. acknowledge support from Intesa
Sanpaolo Innovation Center. F.A.R. acknowledges CNPq (grant 309266/2019-0) and
FAPESP (grant 19/23293-0) for the financial support given for this research.

## References

  * [1] Sebastian Funk, Marcel Salathé, and Vincent A. A. Jansen.  Modelling the influence of human behaviour on the spread of infectious diseases: a review.  J. R. Soc. Interface, 7(50):1247–1256, Sep 2010.
  * [2] J. M. Read, W. J. Edmunds, S. Riley, J. Lessler, and D. A. T. Cummings.  Close encounters of the infectious kind: methods to measure social mixing behaviour.  Epidemiol. Infect., 140(12):2117–2130, Dec 2012.
  * [3] Jared K. Wilson-Aggarwal, Laura Ozella, Michele Tizzoni, Ciro Cattuto, George J. F. Swan, Tchonfienet Moundai, Matthew J. Silk, James A. Zingeser, and Robbie A. McDonald.  High-resolution contact networks of free-ranging domestic dogs canis familiaris and implications for transmission of infection.  PLOS Neglected Tropical Diseases, 13(7):1–19, 07 2019.
  * [4] Joël Mossong, Niel Hens, Mark Jit, Philippe Beutels, Kari Auranen, Rafael Mikolajczyk, Marco Massari, Stefania Salmaso, Gianpaolo Scalia Tomba, Jacco Wallinga, Janneke Heijne, Malgorzata Sadkowska-Todys, Magdalena Rosinska, and W. John Edmunds.  Social Contacts and Mixing Patterns Relevant to the Spread of Infectious Diseases.  PLoS Med., 5(3):e74, Mar 2008.
  * [5] Sergio Arregui, Alberto Aleta, Joaquín Sanz, and Yamir Moreno.  Projecting social contact matrices to different demographic structures.  PLoS Comput. Biol., 14(12):e1006638, Dec 2018.
  * [6] Marcel Salathé, Maria Kazandjieva, Jung Woo Lee, Philip Levis, Marcus W. Feldman, and James H. Jones.  A high-resolution human contact network for infectious disease transmission.  Proc. Natl. Acad. Sci. U.S.A., 107(51):22020–22025, Dec 2010.
  * [7] Ceyhun Eksin, Keith Paarporn, and Joshua S Weitz.  Systematic biases in disease forecasting–the role of behavior change.  Epidemics, 27:96–105, 2019.
  * [8] Paulo Cesar Ventura da Silva, Fátima Velásquez-Rojas, Colm Connaughton, Federico Vazquez, Yamir Moreno, and Francisco A. Rodrigues.  Epidemic spreading with awareness and different timescales in multiplex networks.  Phys. Rev. E, 100(3), Sep 2019.
  * [9] Nathalie Stroeymeyt, Anna V. Grasse, Alessandro Crespi, Danielle P. Mersch, Sylvia Cremer, and Laurent Keller.  Social network plasticity decreases disease transmission in a eusocial insect.  Science, 362(6417):941–945, Nov 2018.
  * [10] Gillian K. SteelFisher, Robert J. Blendon, Mark M. Bekheit, and Keri Lubell.  The Public’s Response to the 2009 H1N1 Influenza Pandemic.  N. Engl. J. Med., 362(22):e65, Jun 2010.
  * [11] Jewel Park.  Changes in Subway Ridership in Response to COVID-19 in Seoul, South Korea: Implications for Social Distancing.  Cureus, 12(4), Apr 2020.
  * [12] Ruth Waitzberg, Nadav Davidovitch, Gideon Leibner, Nadav Penn, and Shuli Brammli-Greenberg.  Israel’s response to the COVID-19 pandemic: tailoring measures for vulnerable cultural minority populations.  Int. J. Equity Health, 19(1):1–5, Dec 2020.
  * [13] Anton Gollwitzer, Cameron Martel, William J. Brady, Philip Pärnamets, Isaac G. Freedman, Eric D. Knowles, and Jay J. Van Bavel.  Partisan differences in physical distancing are linked to health outcomes during the COVID-19 pandemic.  Nat. Hum. Behav., 4(11):1186–1197, Nov 2020.
  * [14] Alun L. Lloyd and Robert M. May.  Spatial Heterogeneity in Epidemic Models.  J. Theor. Biol., 179(1):1–11, Mar 1996.
  * [15] Vittoria Colizza and Alessandro Vespignani.  Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations.  J. Theor. Biol., 251(3):450–467, Apr 2008.
  * [16] Alberto Aleta, Andreia N. S. Hisi, Sandro Meloni, Chiara Poletto, Vittoria Colizza, and Yamir Moreno.  Human mobility networks and persistence of rapidly mutating pathogens.  R. Soc. Open Sci., 4(3):160914, Mar 2017.
  * [17] Marco Ajelli, Bruno Gonçalves, Duygu Balcan, Vittoria Colizza, Hao Hu, José J. Ramasco, Stefano Merler, and Alessandro Vespignani.  Comparing large-scale computational approaches to epidemic modeling: Agent-based versus structured metapopulation models.  BMC Infect. Dis., 10(1):1–13, Dec 2010.
  * [18] Andrea Apolloni, Chiara Poletto, José J. Ramasco, Pablo Jensen, and Vittoria Colizza.  Metapopulation epidemic models with heterogeneous mixing and travel behaviour.  Theor. Biol. Med. Model., 11(1):1–26, Dec 2014.
  * [19] Paolo Bajardi, Chiara Poletto, Jose J. Ramasco, Michele Tizzoni, Vittoria Colizza, and Alessandro Vespignani.  Human Mobility Networks, Travel Restrictions, and the Global Spread of 2009 H1N1 Pandemic.  PLoS One, 6(1):e16591, Jan 2011.
  * [20] Matteo Chinazzi, Jessica T. Davis, Marco Ajelli, Corrado Gioannini, Maria Litvinova, Stefano Merler, Ana Pastore y. Piontti, Kunpeng Mu, Luca Rossi, Kaiyuan Sun, Cécile Viboud, Xinyue Xiong, Hongjie Yu, M. Elizabeth Halloran, Ira M. Longini, and Alessandro Vespignani.  The effect of travel restrictions on the spread of the 2019 novel coronavirus (COVID-19) outbreak.  Science, 368(6489):395–400, Apr 2020.
  * [21] Guilherme S. Costa, Wesley Cota, and Silvio C. Ferreira.  Outbreak diversity in epidemic waves propagating through distinct geographical scales.  Phys. Rev. Res., 2(4):043306, Dec 2020.
  * [22] Alberto Aleta, Qitong Hu, Jiachen Ye, Peng Ji, and Yamir Moreno.  A data-driven assessment of early travel restrictions related to the spreading of the novel COVID-19 within mainland China.  Chaos, Solitons Fractals, 139:110068, Oct 2020.
  * [23] Alberto Aleta and Yamir Moreno.  Evaluation of the potential incidence of COVID-19 and effectiveness of containment measures in Spain: a data-driven approach.  BMC Med., 18(1):1–12, Dec 2020.
  * [24] Matt J. Keeling, Ottar N. Bjørnstad, and Bryan T. Grenfell.  Metapopulation Dynamics of Infectious Diseases.  In Ecology, Genetics and Evolution of Metapopulations, pages 415–445. Academic Press, Cambridge, MA, USA, Jan 2004.
  * [25] Frank Ball, Tom Britton, Thomas House, Valerie Isham, Denis Mollison, Lorenzo Pellis, and Gianpaolo Scalia Tomba.  Seven challenges for metapopulation models of epidemics, including households models.  Epidemics, 10:63–67, Mar 2015.
  * [26] Lin Wang and Xiang Li.  Spatial epidemiology of networked metapopulation: an overview.  Chin. Sci. Bull., 59(28):3511–3522, Oct 2014.
  * [27] Nicola Perra, Duygu Balcan, Bruno Gonçalves, and Alessandro Vespignani.  Towards a Characterization of Behavior-Disease Models.  PLoS One, 6(8):e23084, Aug 2011.
  * [28] Daniela Calvetti, Alexander P. Hoover, Johnie Rose, and Erkki Somersalo.  Metapopulation Network Models for Understanding, Predicting, and Managing the Coronavirus Disease COVID-19.  Front. Phys., 8, Jun 2020.
  * [29] Piero Manfredi and Alberto D’Onofrio.  Modeling the Interplay Between Human Behavior and the Spread of Infectious Diseases.  Springer, New York, NY, New York, NY, USA, 2013.
  * [30] Dan Lu, Alberto Aleta, Marco Ajelli, Romualdo Pastor-Satorras, Alessandro Vespignani, and Yamir Moreno.  Data-driven estimate of SARS-CoV-2 herd immunity threshold in populations with individual contact pattern variations.  medRxiv, page 2021.03.19.21253974, Mar 2021.
  * [31] John R Dormand and Peter J Prince.  A family of embedded runge-kutta formulae.  Journal of computational and applied mathematics, 6(1):19–26, 1980\.
  * [32] Marina Pollán, Beatriz Pérez-Gómez, Roberto Pastor-Barriuso, Jesús Oteo, Miguel A. Hernán, Mayte Pérez-Olmeda, Jose L. Sanmartín, Aurora Fernández-García, Israel Cruz, Nerea Fernández de Larrea, Marta Molina, Francisco Rodríguez-Cabrera, Mariano Martín, Paloma Merino-Amador, Jose León Paniagua, Juan F. Muñoz-Montalvo, Faustino Blanco, Raquel Yotti, Faustino Blanco, Rodrigo Gutiérrez Fernández, Mariano Martín, Saturnino Mezcua Navarro, Marta Molina, Juan F. Muñoz-Montalvo, Matías Salinero Hernández, Jose L. Sanmartín, Manuel Cuenca-Estrella, Raquel Yotti, José León Paniagua, Nerea Fernández de Larrea, Pablo Fernández-Navarro, Roberto Pastor-Barriuso, Beatriz Pérez-Gómez, Marina Pollán, Ana Avellón, Giovanni Fedele, Aurora Fernández-García, Jesús Oteo Iglesias, María Teresa Pérez Olmeda, Israel Cruz, Maria Elena Fernandez Martinez, Francisco D. Rodríguez-Cabrera, Miguel A. Hernán, Susana Padrones Fernández, José Manuel Rumbao Aguirre, José M. Navarro Marí, Begoña Palop Borrás, Ana Belén Pérez Jiménez, Manuel Rodríguez-Iglesias, Ana María Calvo Gascón, María Luz Lou Alcaine, Ignacio Donate Suárez, Oscar Suárez Álvarez, Mercedes Rodríguez Pérez, Margarita Cases Sanchís, Carlos Javier Villafáfila Gomila, Lluis Carbo Saladrigas, Adoración Hurtado Fernández, Antonio Oliver, Elías Castro Feliciano, María Noemí González Quintana, José María Barrasa Fernández, María Araceli Hernández Betancor, Melisa Hernández Febles, Leopoldo Martín Martín, Luis-Mariano López López, Teresa Ugarte Miota, Inés De Benito Población, María Sagrario Celada Pérez, María Natalia Vallés Fernández, Tomás Maté Enríquez, Miguel Villa Arranz, Marta Domínguez-Gil González, Isabel Fernández-Natal, Gregoria Megías Lobón, Juan Luis Muñoz Bellido, Pilar Ciruela, Ariadna Mas i. Casals, Maria Doladé Botías, M. Angeles Marcos Maeso, Dúnia Pérez del Campo, Antonio Félix de Castro, Ramón Limón Ramírez, Maria Francisca Elías Retamosa, Manuela Rubio González, María Sinda Blanco Lobeiras, Alberto Fuentes Losada, Antonio Aguilera, German Bou, Yolanda Caro, Noemí Marauri, Luis Miguel Soria Blanco, Isabel del Cura González, Montserrat Hernández Pascual, Roberto Alonso Fernández, Paloma Merino-Amador, Natalia Cabrera Castro, Aurora Tomás Lizcano, Cristóbal Ramírez Almagro, Manuel Segovia Herná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.  Lancet, 396(10250):535–544, Aug 2020.
  * [33] Flavia Riccardo, Marco Ajelli, Xanthi D. Andrianou, Antonino Bella, Martina Del Manso, Massimo Fabiani, Stefania Bellino, Stefano Boros, Alberto Mateo Urdiales, Valentina Marziano, Maria Cristina Rota, Antonietta Filia, Fortunato D’Ancona, Andrea Siddu, Ornella Punzo, Filippo Trentini, Giorgio Guzzetta, Piero Poletti, Paola Stefanelli, Maria Rita Castrucci, Alessandra Ciervo, Corrado Di Benedetto, Marco Tallon, Andrea Piccioli, Silvio Brusaferro, Giovanni Rezza, 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, Italy, 28 January to 31 March 2020.  Eurosurveillance, 25(49):2000790, Dec 2020.
  * [34] Feinuo Sun, Stephen A. Matthews, Tse-Chuan Yang, and Ming-Hsiao Hu.  A spatial analysis of the COVID-19 period prevalence in U.S. counties through June 28, 2020: where geography matters?  Ann. Epidemiol., 52:54–59.e1, Dec 2020.
  * [35] Ensheng Dong, Hongru Du, and Lauren Gardner.  An interactive web-based dashboard to track COVID-19 in real time.  Lancet Infect. Dis., 20(5):533–534, May 2020.
  * [36] Michele Starnini, Alberto Aleta, Michele Tizzoni, and Yamir Moreno.  Impact of the accuracy of case-based surveillance data on the estimation of time-varying reproduction numbers.  medRxiv, page 2020.06.26.20140871, Jun 2020.