Resilience management during large-scale epidemic outbreaks.

Assessing and managing the impact of large-scale epidemics considering only the individual risk and severity of the disease is exceedingly difficult and could be extremely expensive. Economic consequences, infrastructure and service disruption, as well as the recovery speed, are just a few of the many dimensions along which to quantify the effect of an epidemic on society's fabric. Here, we extend the concept of resilience to characterize epidemics in structured populations, by defining the system-wide critical functionality that combines an individual's risk of getting the disease (disease attack rate) and the disruption to the system's functionality (human mobility deterioration). By studying both conceptual and data-driven models, we show that the integrated consideration of individual risks and societal disruptions under resilience assessment framework provides an insightful picture of how an epidemic might impact society. In particular, containment interventions intended for a straightforward reduction of the risk may have net negative impact on the system by slowing down the recovery of basic societal functions. The presented study operationalizes the resilience framework, providing a more nuanced and comprehensive approach for optimizing containment schemes and mitigation policies in the case of epidemic outbreaks.

Introduction

Data-driven models of infectious diseases[1-15] are increasingly used to provide real- or near-real-time situational awareness during disease outbreaks. Indeed, notwithstanding the limitations inherent to predictions in complex systems, mathematical and computational models have been used to forecast the size of epidemics[16-19], assess the risk of case importation across the world[10,14,20], and communicate the risk associated to uncurbed epidemics outbreaks[21-23]. Despite contrasting opinions on the use of modelling in epidemiology[24], in the last few years a large number of studies have employed them to evaluate disease containment and mitigation strategies as well as to inform contingency plans for pandemic preparedness[11,13,15,24,25]. Model-based epidemic scenarios in most cases focus on the “how many and for how long?” questions. Furthermore, mitigation and containment policies are currently evaluated in the modelling community by the reduction they produce on the attack rate (number of cases) in the population. These studies aim at identifying best epidemic management strategies but typically neglect the epidemic and mitigation impact on the societal functions overall.

The evaluation of vulnerabilities and consequences of epidemics is a highly dimensional complex problem that should consider societal issues such as infrastructures and services disruption, forgone output, inflated prices, crisis-induced fiscal deficits and poverty[26,27]. Therefore, it is important to broaden the model-based approach to epidemic analysis, expanding the purview by including measures able to assess the system resilience, i.e. response of the entire system to disturbances, their aftermath, the outcome of mitigation as well as the system’s recovery and retention of functionality[28-30]. Most important, operationalizing resilience[29-31] must include the temporal dimension; i.e. a system’s recovery and retention of functionality in the face of adverse events[30,32-35]. The assessment and management of system resilience to epidemics must, therefore, identify the critical functionalities of the system and evaluate the temporal profile of how they are maintained or recover in response to adverse events.

Even though the assessment and management of adverse events resilience of complex systems is the subject of active research[32,33,35,36], its integration in the computational analysis of epidemic threats is still largely unexplored[27,37,38].

Here, we introduce a resilience framework to the analysis of the global spreading of an infectious disease in structured populations. We simulate the spread of infectious diseases across connected populations, and monitor the system–level response to the epidemic by introducing a definition of engineering resilience that compounds both the disruption caused by the restricted travel and social distancing, and the incidence of the disease. We find that while intervention strategies, such as restricting travel and encouraging self-initiated social distancing, may reduce the risk to individuals of contracting the disease, they also progressively degrade population mobility and reduce the critical functionality thus making the system less resilient. Our numerical results show a transition point that signals an abrupt change of the overall resilience in response to these mitigation policies. Consequently, containment measures that reduce risk may drive the system into a region associated with long-lasting overall disruption and low resilience. Interestingly, this region is in proximity of the global invasion threshold of the system, and it is related to the slowing down of the epidemic progression. Our study highlights that multiple dimensions of a socio-technical system must be considered in epidemic management and sets forward a new framework of potential interest in analyzing contingency plans at the national and international levels.

Results

We provide a general framework for the analysis of the system-level resilience to epidemics by initially considering a metapopulation network (Fig. 1A). In this case we consider a system made of V distinct subpopulations. These form a network in which each subpopulation i is made of N individuals and is connected to a set k of other subpopulations. A complete description of the networked systems is given in the Methods section. The notation and the description of the parameters used in our simulations are reported in Table 1.

Figure 1

Schematic representation of the metapopulation model. The system is composed of a network of subpopulations or patches, connected by diffusion processes. Each patch contains a population of individuals who are characterized with respect to their stage of the disease (e.g. susceptible, exposed, susceptible with fear, infected, removed), and identified with a different color in the picture. Individuals can move from a subpopulation to another on the network of connections among subpopulations. At each time step individuals move with a commuting rate c from subpopulation i to subpopulation j. (B) Schematic illustration of the system’s critical functionality. The system if fully functional (CF(t) = 1) during ordinary conditions when all the subpopulations are healthy and the number of real commuters is equal to the number of virtual commuters, i.e. D(t) = 0 and C(t) = Z(t). After the outbreak takes place (T0) the system’s functionality decreases because of the disease propagation and the eventual travel reduction. Next the system starts to recover until the complete extinction of the epidemic (T) which corresponds to the time when no more infected individuals are in the system. The curves (a) and (b) represent the critical functionality of scenarios corresponding to high and low values of resilience.

Table 1

Notation and description of the parameters used in our simulations.

Notation	Description
V	Number of subpopulations in the metapopulation network
N	Number of individuals in the system
〈k〉	Average degree of the metapopulation network
D	Number of diseased populations
H	Fraction of healthy populations
A	Fraction of active travelers in the system
p	The parameter that regulates the system wide travel restrictions
r	System’s resilience
CF	System’s critical functionality
T c	Resilience control time
S	Susceptible individuals
S F	Susceptible individuals with fear
E	Exposed individuals
I	Infected individuals
R	Recovered individuals
R 0	Basic reproduction number
λ	The rate at which an ‘exposed’ person becomes ‘infected’
μ	The rate at which an ‘infected’ recovers and moves into the ‘recovered’ compartment
β	The parameter controlling how often a ‘susceptible’-‘infected’ contact results in a new ‘exposed’
β F	The parameter controlling how often a ‘susceptible’-‘infected’ contact results in ‘susceptible individual with fear’
α	The parameter controlling how often a ‘susceptible’-‘susceptible individuals with fear’ contact results in a new ‘susceptible individual with fear’
r b	The parameter that modulates the level of self-induced behavioral change that leads to the reduction of the transmission rate
μ F	The rate at which individual with fear moves back into the ‘susceptible’ compartment

Diffusion Processes

The edge connecting two subpopulations i and j indicates the presence of a flux of travelers i.e. diffusion, mobility of people. We assume that individuals in the subpopulation i will visit the subpopulations j with a per capita diffusion rate d on any given edge[39] (see the Methods section for further details). We define the total number of travelers Z between the subpopulations i and j at time t as Z(t) = dN(t), so that when the system is fully functional, the total number of travelers at time t from the node i is . Under these conditions, the total number of travelers in the metapopulation system at time t is simply Z(t) = ∑Z(t). In the following we assume that infected individuals do not travel between subpopulations, thus reducing the actual number of travelers.

Reaction Processes

In analyzing contagion processes we extend the compartmental scheme of the basic SEIR model[40,41] (see Methods and Supplementary Information (SI) for a detailed description). Indeed an important element in the mitigation of epidemics is self-initiated behavioral changes triggered in the population by awareness/fear of the disease[42,43]. These generally reduce the transmissibility and spreading. Examples of behavioral changes include social distancing behaviors such as avoidance of public places, working from home, decrease of leisure and business travel etc. In order to include behavioral changes in our model, we consider a separate behavioral class within the population[44], defining a special compartment of susceptible individuals, S, where F stands for “fearful”. In particular, individuals transition to this compartment depending on the prevalence of infected and other fearful individuals according to a rate β. This rate mimics the likelihood that individuals will adopt a different social behavior as a result of the increased awareness of the disease as perceived from the number of infected and fearful individuals present in the system. Clearly, spontaneous or more complex types of transitions (for example indirectly linked to the disease transmission due to mass media effects[44]) could be considered. However, they would require more parameters and introduce other non-trivial dynamics. We leave the study of other behavioral changes models for future works. It follows that in each subpopulation the total number of individuals is partitioned into the compartments S(t),S(t),E(t),I(t),R(t) denoting the number of susceptible, fearful, exposed, infected, and removed individuals at time t, respectively. The transition processes are defined by the following scheme: S + I → E + I, S + I → S + I, S + S → 2S, S + I → E + I, E → I and I → R with their respective reaction rates, β, β, αβ, rβ, λ and μ. Analogously, individuals in the S compartment may transition back in the susceptible compartment with a rate μ, S + S → S. The model reverts to the classic SEIR if β = 0 (the detailed presentation of the dynamic is reported in the SI). The basic reproductive number of an SEIR model is R0 = β/μ. This quantity determines the average number of infections generated by one infected individual in a fully susceptible population. In each subpopulation the disease transmission is able to generate a number of infected individuals larger than those who recover only if R0 > 1, yielding the classic result for the epidemic threshold[45]; if the spreading rate is not large enough to allow a reproductive number larger than one (i.e., β > μ), the epidemic outbreak will affect only a negligible portion of the population and will quickly die out (the model details are reported in the Methods section).

System’s resilience

Here, we introduce a quantitative measure that captures and implements the definition of resilience in epidemic modelling, similarly to what proposed in Ganin et al.[32,34]. Among the many possible elements defining the resilience of a system, we consider the system-wide critical functionality as a function of the individual’s risk of getting the disease and the disruption to the system’s functionality generated by the human mobility deterioration. For the sake of simplicity, in our model we assume that infected individuals do not travel. The extension to models in which a fraction of infected individuals are traveling is straightforward[4] with the only effect of decreasing the timescale for the disease spreading, but not altering the overall dynamic of the system. Furthermore, as discussed below, the system might be subject to other travel limitations. As a result, during the epidemic we have an overall decrease in the mobility flows with respect to a disease-free scenario. It follows that the number of travelers between subpopulations i and j at time t is , where c is the adjusted diffusion rate, , and the total number of commuters in the metapopulation system at time t is given by . Note that in general, c < d. This can be naturally related to a deterioration of the system-level critical functionality as it corresponds to economic and financial losses as well as logistic and infrastructural service disruption. In order to evaluate the system’s loss of critical functionality related to the travel restrictions, we define the fraction of active travelers at time t as A(t) = C(t)/Z(t).Analogously, we characterize the system’s risk related to the disease propagation as the fraction of healthy subpopulations H(t) = 1 − D(t)/V, where D(t) is the number of diseased subpopulations at time t and Vis the total number of subpopulations in the system. The number of diseased subpopulations accounts for the amount of risk posed to individuals in the system, which we assume to be proportional to the overall attack rate and expresses the vulnerability of the networked system[36,46,47]. Here, as the model assumes statistically equivalent subpopulations, the attack rate is proportional to the number of subpopulations affected by the epidemic. At time t, we define the critical functionality, CF(t), (Fig. 1B), as the product of the fraction of active travelers A(t) and the fraction of healthy populations H(t), i.e. CF(t) = H(t) ⋅ A(t). Per our earlier definition of resilience[32] r, we evaluate it as the integral over time of the critical functionality, normalized over the control time T so that r ∈ [0,1]:The control time T corresponds to the maximum extinction time T for different values of epidemic reproductive number R0 (see the Supporting Information for further detail). Resilience, therefore, also includes the time dimension, in particular, the time to return to full functionality, as defined by the system’s critical elements. In reference[32] we provided an operational definition of resilience starting from the concepts advanced by the National Academy of Sciences in USA. In this paper, we apply such general framework to the case of disease spreading. Furthermore, we extend it to reaction-diffusion processes on metapopulations. In the following, we will quantitatively characterize different containment/mitigation interventions via a critical functionality analysis. Desirable (optimal) strategies correspond to high (maximum) value of r. It is worth remarking that, for the sake of simplicity, we use here a definition of critical functionality that weights equally the two components A(t) and H(t). Thus, our findings are constrained by such choice. The two contributions could be weighted differently, i.e. CF(t) = H(t) ⋅ A(t). However, our aim is to highlight the importance of going beyond “model-based” approach to epidemic analysis and move towards system resilience assessments. In this spirit, we opted for the simplest definition of critical functionality able to capture the two most used metrics in model-based approaches: epidemic risk and mobility. We used the multiplication of the two quantities because it makes the critical functionality more sensitive to small changes of the values. Furthermore, by multiplying two ratios we don’t need to add a normalization factor (the critical functionality is defined in the interval [0,1]). In more realistic context, and depending on the precise cost-benefit analysis, the various terms may be weighted differently and more complex functional form for the critical functionality can be defined. Among other things, these type of analysis could consider: (i) the details of the disease spreading in the population such as mortality, infectiousness, recovery time, and possible residual immunity (ii) the preparedness, measured in terms of availability of vaccines, anti-virals, hospital beds, or intensive care units, (iii) the socio-economical costs induced by a major outbreak and by interventions such as travel bans, school closures etc. (iv) politics and public perception of risk.

Effects of system-wide travel restrictions

Epidemic containment measures, based on limiting or constraining human mobility, are often considered in the contingency planning and always re-emerge when there are new infectious disease threats[1]. The target of these control measures are travels to/from the areas affected by an epidemic outbreak and the corresponding decrease of infected individuals reaching areas not yet affected by the epidemic. At the same time, travel restrictions have a large impact on the economy and affect the delivery of services, including medical supplies and the deployment of specialized personnel to manage the epidemic. For this reason, travel restrictions must be carefully scrutinized to trade off the costs and benefits. We introduce the parameter p∈[10−5,1] that allows us to simulate policy-induced system-wide travel restrictions. In our settings, such measures are active until the disease is circulating in the system, i.e. there is at least one infected individual across all subpopulations. In the case of no travel restrictions and/or after the disease dies out, we have p = 1. In the case of travel restrictions (p < 1), we rescale travel flow so that mobility is a fraction of that in the unaffected system; i.e. c = p ⋅ d. To better understand the effect of such mitigation strategy, let us consider the classic SEIR model by setting β = 0. In the presence of travel restrictions and depending on the level of mixing, each subpopulation may or may not transmit the infection or contagion process to another subpopulation it is in contact with. In other words, the mobility parameter p influences the probability that exposed individuals will export the contagion process to other regions of the metapopulation network. Further, it introduces a transition between a regime in which the contagion process may invade a macroscopic fraction of the network and a regime in which it is limited to a few subpopulations. The transition is mathematically characterized by the global invasion threshold R*[45]. This is the analogue of the basic reproductive number at the subpopulations level and defines the average number of infected subpopulations generated by one infected subpopulation in a fully susceptible metapopulation system. In general, R* is a function of the basic epidemic parameters, including R0, and the mobility parameter p. The invasion threshold occurs at the critical value p for which R* = 1. In some cases, p can be evaluated analytically (see the Methods section). In general, it can be estimated numerically by measuring the number of infected subpopulations as a function of the parameter p.

Risk, as measured in terms of attack rate, is, therefore, monotonically decreasing due to increasingly restricted travel, and falls to virtually zero for values of p below the invasion threshold. Thus, from a risk perspective, the best strategy during a disease outbreak is to reduce the mobility. However, an inspection of the profile of resilience provides a different picture. In Fig. 2 we report the value of r obtained by sampling the phase space of the model p−R0 for different values of the travel diffusion parameter and the epidemic reproductive number in heterogeneous metapopulation systems (a comparison between homogenous and heterogeneous networks is reported in SI). Each point of the phase space is studied by performing 100 stochastic realizations. The 3D dimensional plot in the p,R0,r space reported in Fig. 2A indicates that the overall resilience profile is characterized by a sharp drop as we approach the invasion threshold, i.e. p → p.

Figure 2

Resilience and final fraction of diseased populations in the heterogeneous metapopulation system with traffic dependent diffusion rates. (A) 3D surface representing resilience in a homogeneous metapopulation system as a function of local threshold R0and the diffusion rate p: the minimum value of resilience separates two regions associated to values very close to the optimal case. (B) Cross-sections (blue) of the 3D plot for R0 = 3.5 and its comparison with the final fraction of diseased populations (red): while the reduction of the diffusion rate p brings to a constant the fraction of diseased populations it also causes an initial decrease of resilience to a minimum value after which it starts increasing and the system returns to its optimal conditions. (C) The map of the final fraction of diseased populations D∞/V is shown as a function of the local epidemic threshold R0 and the travel diffusion p. We show that the minimum values of resilience (blue points) correspond to the theoretical value of the final fraction of diseased subpopulations D∞/V at the end of the global epidemic (black line).

Figure 2B shows that, while the risk decreases, the reduction of the diffusion rate p causes a reduction of r until the global invasion threshold, after which the resilience value rapidly increases. This effect is mainly due to the critical slowing of the epidemic spreading near the invasion threshold. Indeed, close to the threshold, the epidemic is still in a supercritical state, but it takes increasingly longer time to invade the system as the threshold is approached. This can be simply related to the divergence of the invasion doubling time T, which is defined as the time until the number of infected subpopulations doubles, relative to that at some other time. The doubling time is related to the subpopulation reproductive number as T ~ (R* − 1)−1, leading to a divergence of the doubling time as the invasion threshold is approached for R* → 1. Although the absolute risk is very low, the system remains in a state of deteriorated functionality (restrictions in travels) for longer and longer times[48]. The decrease of functionality is not offset by a corresponding decrease of risk, and the minimum in resilience is attained exactly at the global invasion threshold. The comparison between the theoretical values of the invasion threshold and the computed minimum values of resilience is reported in Fig. 2C.

Effects of self-initiated behavioural changes

In order to isolate the effects of behavioural changes, in this section the travel parameter is kept constant with p = 1. Individuals in the S compartment adopt travel avoidance so that β plays a similar role to the travel restriction as reported in Fig. 3. Furthermore, inside each subpopulation, individuals in the S compartment reduce their contacts, thus decreasing the likelihood to become infected. Overall, the presence of self-initiated behavioral changes in a population results in a reduction of the final epidemic size. In this setting, we have explored a phase space of parameters constituted of R0∈[1.01,3] and β∈[0,20] (see the Methods section for the other model parameters). In Fig. 3A we quantify resilience for different values of the fear parameter β in heterogeneous metapopulation systems. The 3D dimensional plot in the β,R0,r space shows a clear similarity with the travel restrictions scenario. Figure 3B shows that, while increasing β leads to a decrease in risk, it also induces a reduction of resilience. It is possible to observe that, even in this case, the minimum values of r are related to the invasion threshold. In Fig. 3C the phase diagram of the fraction of diseased populations at the end of the simulations D∞/V is reported in the β,R0 space. This picture shows that there is a critical value of the fear transmissibility parameter β, after which the fraction of diseased populations D∞/V starts to decrease (i.e. D∞/V < 1). The minimum value of resilience, in this case, corresponds to the value of the fear transmissibility, after which a reduction of the fraction of diseased populations is observed. Although the approach to this critical boundary corresponds to a reduction of the infection risk, similarly to the case of travel restrictions, the measured resilience of the system decreases and attains its minimum value right at the transition point.

Figure 3

Resilience and diseased populations in a heterogeneous metapopulation system with individual self-dependent travel reduction. (A) 3D surface representing resilience in a heterogeneous metapopulation system as a function of local threshold R0 and the fear parameter β: two areas of high values of resilience are separated with a narrow region of very low ones. (B) Comparison between resilience (blue) reported as cross-sections of the 3D plot for R0 = 1.3 and the final fraction of diseased populations D∞/V (red): while the increase of the fear transmissibility parameter β brings to a constant the fraction of the diseased populations it also causes an initial decrease of resilience to a minimum value after which the system bounces back to optimal conditions. (C) Even in this case the minimum values of resilience (blue points) correspond to the transition region from high to low final diseased populations. The colormap of the logarithmic of the healthy populations (log(1 − D∞/V)) is shown as a function of the local epidemic threshold R0 and the fear parameter β.

Effects of system-wide travel restrictions in data-driven simulations

As a further confirmation of the validity of the theoretical construct above described, we tested our results in a data-driven modelling setting. We considered the Global Epidemic and Mobility model (GLEAM)[3,49] which integrates high resolution demographic and mobility data by using a high-definition, geographically structured metapopulation approach. The model’s technical details and the algorithms underpinning the computational implementation have been extensively reported in the literature. GLEAM is a spatial, stochastic and individual-based epidemic model that divides the world population into geographic census areas, defined around transportation hubs and connected by mobility fluxes. The population of each census area is obtained by integrating data from the high-resolution population database of the ‘Gridded Population of the World’ project of the Socioeconomic Data and Application Center at Columbia University (SEDAC). The mobility among subpopulations is comprised of global air travel and the small-scale movement between adjacent subpopulations; i.e., the daily commuting patterns of individuals. Commuting and short-range mobility considers data from 80,000 administrative regions in 5 different continents. Here, we considered the Continental United States and simulated an SEIR contagion process, in which infected individuals do not travel. The number of infected subpopulations at the end of an outbreak and resilience as a function of the global mobility restrictions that result are shown in Fig. 4. The initial conditions of the epidemic were set with 5 infected individuals in the city of New York, assuming β = 0.48, λ = 0.66 and μ = 0.45. Mobility restrictions are implemented by reducing all the mobility flows connecting diseased subpopulations by a factor p, thus considering the heterogeneities of the subpopulations due to their different local mobility patterns (see SI). The control time T used in the calculation of r corresponds to the epidemic extinction time for the different values of the diffusion rate.

Figure 4

Resilience and epidemic size in the data-driven scenario. (A) The plot shows the difference between resilience (blue) and the final fraction of diseased populations (red) for different values of the diffusion rate p. Here, we can identify three critical regions of the system. (i) diffusion rate p = 0.1 above the critical invasion threshold. Even if the system is characterized by sub-optimal resilience, the disease spreads all over the system. (ii) the reduction of the diffusion parameter p results in a significant decrease of the number of diseased populations but also in a dramatic decrease of resilience; (iii) below the critical invasion threshold resilience goes back to high values as fraction of diseased populations approaches zero. (B) Epidemic size (red) and resilience (blue) for the different values of the diffusion parameter p corresponding to the three aforementioned regions. Python 2.7 (https://www.python.org/) and the Basemap library (https://pypi.python.org/pypi/basemap/1.0.7) were used to create these maps.

As with the theory-driven model here we observe that a reduction of the travel diffusion p brings a constant reduction of diseased populations, but also reduces resilience until a critical value p = 1.2 ⋅ 10−4. In Fig. 4B we illustrate the geographical spreading of the contagion process and the reduction of traveling of each subpopulation tracked by the model in the Continental USA for values of p corresponding to three different regions of the diagram of Fig. 4A. The figure clearly illustrates three regimes: i) for low travel reduction, a very severe epidemic hits all the subpopulations, but the short duration allows the system to go back to normal in a short time (high values of resilience); ii) for travel reduction close to the global invasion threshold, a small number of subpopulations are hit but the system critical functionality is compromised for a very long time, thus, resulting in a low values of resilience; iii) travel reduction above the critical threshold allows the system to contain the epidemic at the origin with low risk and high values of resilience. It is worth remarking that in the data-driven model, the minimum value of resilience is reached for travel restrictions that correspond to a reduction of mobility of three to four orders of magnitude. This is because in modern transportation networks the global invasion threshold, as already pointed out in other studies[10,14,20,39], is reached only for very severe travel restrictions that are virtually impossible to achieve. In other words, in realistic settings the feasible increase of travel restrictions appears always to decrease resilience. This calls for a careful scrutiny of the trade-off between individual’s risk and resilience, as the region where both are achieved is virtually not accessible.

Discussion

The realistic threat quantification is difficult and evaluation of vulnerabilities and consequences of new disease epidemics is certainly a challenge. We analyzed the impact of an infectious disease epidemic in structured populations by considering a definition of system resilience that takes into consideration not only the number of infected individuals but also society’s need for maintaining certain critical functions in space and time[37]. In particular, we assume that the limitations and disruptions of human mobility deteriorate the system’s functionality. We observe that containment measures, that limit individuals’ mobility, are advantageous in reducing risk but may deteriorate the system’s functionality for a very long time and thus correspond to low resilience. Although we have considered only two of the many dimensions encompassing the functionality of socio-technical systems[28,30], we show that study of resilience allows stakeholders to measure the impact of epidemic threats and differentiate between different management alternatives. It is straightforward to envision more realistic definition of the critical functionality. The components of critical functionality could be weighted according to objective/subjective evaluation of their relevance to stakeholders. Finally, cost-benefit analyses and ethical considerations should be included in the analysis of the societal impacts of disease that could lead to long lasting effects or even death of the affected individuals. This study highlights the importance of resilience-focused analysis for selecting intervention strategies. The natural tendency to be conservative in managing potentially inflated risks associated with new and emerging epidemics can result in unnecessary burdensome and possibly ineffective actions like quarantines as well as travel bans[50]. The emerging field of resilience assessment and management[29] and its implementation[32,34,35] could thus evaluate cross-domain alternatives to identify a policy design that enhances the system’s ability to (i) plan for such adverse events, (ii) absorb stress, (iii) recover, and (iv) predict and prepare for future stressors through necessary adaptation. To this end, the framework we presented can be of potential use for optimizing the policy response to a disease outbreak by balancing risk reduction with the disruption to critical functions that is associated with public health interventions.

Methods

Disease propagation and self-initiated behavioral changes

The metapopulation system is described by a scale-free network (SF) with a power-law degree distribution P(k) ∼ k(−, which is generated by the configuration model[51] with the minimum degree m = 2, γ = 2.1. (For the travel restriction scenario, in the SI, we report a comparison of the results between the heterogeneous networked system described above and a metapopulation system formed by a random network with Poisson degree distribution, which is generated by the Erdos–Rényi (ER) model[52]). The networks have V = 5000 nodes and average degree 〈k〉 ∼ 6, while the total number of individuals is N = V 2 = 25. 106 which are distributed among the subpopulations nodes proportional to their degree distribution. At the beginning, 10 populations are selected at random and 50 individuals are set as exposed. All other individuals across the system are initially susceptible. We study a compartmental scheme that extends the basic SEIR[40] model by considering separate behavioral classes within populations (see SI for the detailed description of the model). For this reason, we assume that individuals can spontaneously change their behavior because of the fear of the disease entering in a specific compartment S of susceptible individuals. In the case of travel restrictions, we set the transition rate from exposed to infected λ = 0.67 days−1 and recovery rate μ = 0.33 days−1. In the case of the behavioral model, we set the disease parameters λ = 0.3 days−1 and μ = 0.1 days−1 while we consider an infection probability reduction r = 0.15, the self-reinforcement parameter α = 0.1 and the relaxation parameter μF = 0.5. All the presented results are averaged over 100 simulations.

Mobility process

We adopt a simplified mechanistic approach that uses a Markovian assumption for modeling migration among subpopulations; at each time step, the movement of individuals is given according to a matrix d that expresses the probability that an individual in the subpopulation i is traveling to the subpopulation j. We assume that the diffusion rate on any given edge from a subpopulation node of degree k to a subpopulation node of degree k is proportional to k[39] and it is given by d = ω0 (kk)/T, where T = ∑ w = ∑ ω0(kk) represents the total flow in i, while θ and the exponent ω0 are system specific (e.g., and θ = 0.5 and in the world-wide air transportation network[53]). In this scenario, we consider θ = 0.5 and ω0 = 10−3.

Global invasion threshold

For the SEIR model it is possible to explicitly calculate the average number of infected subpopulations for each infected subpopulation in a fully susceptible metapopulation system as [3] where represents the average number of individuals in a subpopulation. The condition R* = 1 defines the invasion threshold for the system. Only for R* > 1 can the epidemic spread to a large number of subpopulations. The invasion threshold readily provides an explicit condition for the critical mobility p, below which an epidemic cannot invade the metapopulation system, yielding the equation [39].