Vector-borne epidemics driven by human mobility

PHYSICAL REVIEW RESEARCH 2, 013312 (2020) Editors’ Suggestion Vector-borne epidemics driven by human mobility David Soriano-Paños ,1,2Juddy Heliana Arias-Castro,3Adriana Reyna-Lara ,1,2Hector J. Martínez ,3 Sandro Meloni ,4and Jesús Gómez-Gardeñes 1,2,5 1GOTHAM laboratory, Institute for Biocomputation and Physics of Complex Systems (BIFI), University of Zaragoza, 50018 Zaragoza, Spain 2Departamento de Física de la Materia Condensada, Universidad de Zaragoza, 50009 Zaragoza, Spain 3Department of Mathematics, Universidad del Valle, 760032 Santiago de Cali, Colombia 4IFISC, Instituto de Física Interdisciplinar y Sistemas Complejos (CSIC-UIB), 07122 Palma de Mallorca, Spain 5Center for Computational Social Science (CCSS), Kobe University, 657-8501 Kobe, Japan (Received 16 December 2019; accepted 24 February 2020; published 13 March 2020) Vector-borne epidemics are progressively becoming a global burden, especially those related to flaviviruses, and the effects of different factors such as climate change or the increase of human mobility can sensibly increase the population at risk worldwide. Such outbreaks are the result of the combination of different factors including crossed contagions between humans and vectors, their demographic distribution and human mobility among others. The current availability of information about all those ingredients demands their incorporation into current mathematical models for vector-borne disease transmission. Here, relying on a Markovian formulation of the metapopulation dynamics, we propose a framework that explicitly includes human-vector interactions, mobility, and demography. The analysis of the framework allows us not only to derive an expression of the epidemic threshold capturing the conditions for the onset of the epidemics but also to highlight some unseen features of vector-borne epidemics, such as abrupt changes in the unfolding patterns of the disease for small variations of the degree of mobility. Finally, driven by these insights, we obtain a prevalence indicator to rank populations according to their risk of being affected by a vector-borne disease. We illustrate the utility of this indicator by reproducing the spatial distribution Dengue cases reported in the city of Santiago de Cali (Colombia) from 2015 to 2016. DOI: 10.1103/PhysRevResearch.2.013312 I. INTRODUCTION The explosive dissemination of Zika virus across the Americas has been one of the major concerns of public health organizations across the world in recent years [1]. Zika’s global threat is, unfortunately, the last example of the extremely rapid dissemination of mosquito-borne flaviviruses over the past two decades. From Dengue to Zika, through West Nile and Chikungunya viruses, more than one billion people are infected and more than one million people die from vector-borne diseases (VBD) every year [2]. According to the World Health Organization (WHO), VBD are responsible for one-sixth of the illness worldwide and more than half of the human population live in risk areas for these diseases [3]. The threat of emergent VBD in tropical and equatorial regions progressively span across more temperate areas as a byproduct of climate change. As temperature rises, the areas that are conducive to mosquitoes expand, meaning more opportunities for VBD to spread [4–9]. Another important driver of many recent VBD outbreaks is the explosion of human Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. mobility we have witnessed over the last decades, including urban, regional and long-range displacements [10–14]. In spite of the social and economic benefits, the expansion of human mobility networks produced, as a consequence, the speed-up of epidemic waves and the emergence of correlated outbreaks in faraway regions. When facing these VBD outbreaks, public health systems have to respond quickly and efficiently to mitigate the spread of the disease. Although vaccines are a long-term solution, up to now, the ones available for some flaviviruses have limited impact [15,16]. Thus the most common way for preventing local Aedes-VBD outbreaks is tackling the vectors by the use of pesticides, larvivorous fishes, or Wolbachia bacteria [17]. The use of geolocalized control strategies, however, may not always be as effective when facing the threat of a globalscale pandemic [18]. On the contrary, the fast transcontinental movement of VBD demands a coordinated action of all the involved actors for the efficient use of local control means. This implies taking into account that those populations at risk are not isolated and, as for human-human transmission diseases [19], human mobility plays a key role in the spread of VBD across different populations [20]. Therefore, incorporating human mobility into disease transmission models has become a must when proposing mathematical frameworks aimed at capturing the contagions patterns observed in actual epidemic scenarios. So, the creation of new policies could rely on those refined models to generate better strategies involving human mobility to contain epidemics spreading. 2643-1564/2020/2(1)/013312(12) 013312-1 Published by the American Physical Society DAVID SORIANO-PAÑOS et al. PHYSICAL REVIEW RESEARCH 2, 013312 (2020) Metapopulation models, originally proposed in the field of ecology [21–23], enable the mixing of mobility and contagion dynamics into a single formulation. These models can be described as networks in which nodes account for geographic locations (such as neighborhoods in cities, cities within countries, etc.), i.e., subpopulations where large collectivities live and interact. In addition, the links of the network represent (and quantify) trips made by individuals between different subpopulations. From the first studies making use of metapopulations for the study of infectious diseases transmission [24–29], the field has advanced both in its theoretical grounds [30–35] and its use for large-scale agent-based simulations [36–38]. The latter approach incorporates many realistic aspects of human interactions with the goal of being useful for making epidemic forecast and the design of efficient prevention policies. On the other hand, the theoretical part has been spurred by the increasing spatiotemporal resolution of current data gathering techniques and many efforts have been devoted to bridge the gap between theory and realistic models during the last years [39]. First attempts in this direction involved displacement kernels [40,41] to model a local range of movement around an area that, in the last years, led the way to the inclusion of more sophisticated mobility patterns, such as the commuting nature of human mobility [42–46], the high order memory of human displacements [47], or the coexistence of different transportation behaviors [48]. These sophisticated theories allow to capture the temporal and geographical spread of diseases while providing insights about the mechanisms driving the observed patterns. In the case of VBD, the use of metapopulation models has been recently fostered due to recent outbreaks such as Zika and Chikungunya. This way, epidemic models for VBD transmission have abandoned mean-field and well mixed hypothesis to consider patchy environments subjected to human flows. On the theoretical side, and pretty much as for metapopulations models of human-human transmission diseases, the frameworks rely on important assumptions that allow its analytical study. One of these assumptions is to consider the random diffusion of humans across patches [49–51] or displacement kernels [52] instead of realistic mobility patterns. On the other hand, when actual ingredients of human mobility, such as its recurrent nature, are taken into account both random diffusion [53,54] and displacement kernels models [55] fail to provide insights about the role that real mobility networks play on the transmission of VBD. Thus metapopulation theories are still far from incorporating the many aspects influencing the onset of VBD outbreaks and lack the predictive power provided by data-driven agent-based simulations [56,57]. The main goal of this work is to provide a benchmark that allows the study of large-scale vector-borne epidemics in a unified way and, more importantly, enabling the test of coordinated control strategies at the light of available data about vector incidence together with human demography and mobility datasets. To this aim, we first elaborate a metapopulation model for the transmission of VBD that allows us to derive the conditions under which epidemics take place. The analytical expression of the epidemic threshold is revealed by a matrix encoding the probability of crossed infections between humans and vectors of different subpopulations. Importantly, the spectral analysis of this matrix reveals the risk associated with each patch, pointing out those subpopulations triggering the epidemic onset. We confirm these results by testing synthetic metapopulations and a real case, the city of Cali (Colombia). To round off, taking advantage of the theoretical insights and analytical expressions provided by the formalism, we propose a metric capturing the risk associated with each patch. We implement this metric in the city of Cali obtaining a very good agreement between the estimated risk and the actual distribution of Dengue incidence across districts, highlighting the important role of recurrent human mobility patterns for explaining the spatial dissemination of VBD. II. METAPOPULATION MODEL FOR VBD TRANSMISSION In the following, we focus on the description of a vectorborne contagion dynamics in a complex metapopulation. To this aim, we consider a set of Npopulations or patches in which contagion processes occur governed by the RossMacdonald (RM) model. The RM model [58–60] assumes that both vectors and individuals can be either susceptible of contracting the disease or infectious entities capable of fostering the disease. In this sense, this model was one of the first attempts to address the impact of human-vector crossed contagions on the spread of diseases. However, let us remark that, the original RM model has some limitations such as assuming that VBD do not confer long-lasting immunity to infected individuals after overcoming the disease or including a closed population of both humans and vectors. Nevertheless, despite these assumptions, it is able to capture the essential ingredients involved in the transmission of many VBD that do not confer immunity after infection, e.g., Malaria, or those diseases conferring immunity but having several co-circulating serotypes within the same population, e.g., Dengue. Additionally, here, we mainly focus on VBD provoked by flaviviruses, in which humans act as reservoir of the disease and mosquitoes are the vectors intermediating its transmission among individuals. We choose the original formulation of the RM model since it constitutes the minimal model of crossed infections enabling the analytical study of the interplay between human mobility and the spread of VBD. However, different disease-specific contagion patterns, such as humanhuman contagions in Zika [61,62], different epidemiological features such as passive immunity [63], or the effects of voluntary vaccination [64,65] can still be easily accommodated by introducing small refinements to our original framework. The relevant variables of the RM dynamics are (i) the fraction of infected humans at time t,ρH(t), and (ii) the fraction of vectors infected at time t,ρM(t). The evolution of these two variables is given as a product of the elementary processes described in Fig. 1(a). Namely, susceptible humans become infected with probability λMH after being bitten by an infected vector whereas healthy vectors become infectious with probability λMH when interacting with an infected human. In addition, we assume that each vector makes a number of βcontacts with (healthy or infected) humans. This way, no human-human or vector-vector direct infections are allowed. Finally, infected humans become susceptible with probability 013312-2 VECTOR-BORNE EPIDEMICS DRIVEN BY HUMAN … PHYSICAL REVIEW RESEARCH 2, 013312 (2020) FIG. 1. Ross-Macdonald model and metapopulation approach. (a) Schematic representation of the processes described in the RM model. Relevant parameters are (i) the probability that an infected vector transmits the disease to a healthy individual, λMH; (ii) the probability that an infected human transmits the disease to a healthy vector, λHM; (iii) the feeding rate of vectors β; (iv) the probability that an infected human recovers, μH; and (v) the mortality rate of vectors, μH. In (b) and (c), we show schematically the two basic stages of each Monte Carlo step in our metapopulation approach. As shown in (b), individuals are associated with one of the 3 nodes of a network. Namely, starting from the subpopulation in the top (1) and following clockwise we have populations composed of 8, 10, and 4 humans with 5, 3, and 9 vectors respectively. Once the movement has been done [see (c)] individuals mix and, consequently, the instant populations of humans at each node change to 6, 8, and 8. The RM dynamics then takes place among the individuals and vectors coexisting at that moment in the same subpopulation. Finally, after these interactions, the individuals go back to their respective associated nodes and the configuration is again the one of (b). Note that we consider that vectors do not move from their corresponding node. μH, while (healthy or infected) vectors die with probability μM, being replaced by newborn healthy ones. Although the RM dynamics captures the elementary contagion processes taking place inside each population, the dynamical evolution of each patch depends strongly on the others, since they are not isolated. On the contrary, many individuals with residence in one subpopulation may visit others during, for instance, their daily commuts to other geographical locations. On the other hand, the mobility of vectors is rather limited, hence they are assumed not to move from their original population. This assumption is valid for many VBD such as Dengue, Zika, or Chikungunya since their carriers, Aedes mosquitoes, typically fly an average of 400 meters [66]. Thus it is the mobility of infectious individuals (who pass the diseases to healthy vectors living in distant subpopulations) what triggers the propagation of local disease outbreaks across the whole system. To characterize the mobility of individuals, we denote each area as a patch in a metapopulation. Each patch ihas a population of niindividuals and mivectors and, importantly, they may be different from one population to the other, as they are derived from the demographic partition of the population and the observed vector prevalence in each patch. Each individual is associated to one subpopulation, say i, considered as her residence. This way, the population niof the patch iis the number of individuals whose residence is node i. Simultaneously, nodes are connected in pairs forming a complex weighted and directed network encoded in an adjacency matrix R, whose entries Rij account for the probability that a trip departing from patch ihas as destination population j [see Fig. 1(b)]. Matrix Rcan be computed from the observed number of trips between each pair of nodes (i,j), Wij,as Rij =Wij N l=1Wil .(1) Thus matrix Rencodes the information provided by mobility datasets. The former two dynamical processes at work—RM dynamics and human movements—interplay in each time step of the metapopulation dynamics as follows. We start with a small fraction of infected humans and/or vectors. The initial quantity, being small, can be homogeneously distributed across the populations or localized in one or few nodes in case of being interested in determining those patches boosting epidemic spreading in its early stage. Once the initial infectious seed has been placed, the following microscopic processes are considered. (1) At each time step, t, healthy agents decide to move from their residence with probability por remaining in it with probability (1 −p). Moreover, as symptoms associated to some VBD are severe, we also include the possibility of rescaling the infected agent mobility to αpwith α∈[0,1]. (2) If an agent leaves her residence, say i, she goes to a different subpopulation chosen among those connected to i. The choice is dictated by matrix Rin Eq. (1), being Rij the probability of moving from ito subpopulation j. (3) Once all the individuals have been placed in the patches [see Fig. 1(c)], humans and vectors that are currently at the same patch interact as dictated by the RM model. This way, both humans and vector update their dynamical states [susceptible or infected as shown in Fig. 1(a)] as a result of the contagion and recovery processes. 013312-3 DAVID SORIANO-PAÑOS et al. PHYSICAL REVIEW RESEARCH 2, 013312 (2020) (4) Once the epidemic state of the agents have been updated, each individual moves back to her residence and the process starts again for time t+1. III. MODEL EQUATIONS Once defined the basic steps of the mechanistic simulations, we now tackle the mathematical formulation of the processes described above. To this aim, for each patch i (i=1,...,N) we have two variables: the probabilities that humans with residence in i,ρH i(t), and vectors associated to i,ρM i(t), are infectious at time t, respectively. These 2Nvariables evolve according to the following Markovian equations: ρH i(t+1) =ρH i(t)(1 −μH)+1−ρH i(t)IH i(t),(2) ρM i(t+1) =ρM i(t)(1 −μM)+1−ρM i(t)IM i(t),(3) where IH i(t) and IM i(t) account for the probability that a healthy human with residence in subpopulation iand a healthy vector associated to iare infected at time trespectively. The former infection probability reads IH i(t)=(1 −p)PH i(t)+p N  j=1 RijPH j(t),(4) where PH i(t) is the probability that an agent placed in population iat time tis infected. This probability can be written as PH i(t)=1−1−λMHρM i1 neff i(ρH(t),α,p)βmi .(5) Finally, neff i(ρh(t),α,p), which is the number of humans placed in (but not necessarily residing in) population ican be expressed as: neff i(ρH(t),α,p)=1−p1−(1 −α)ρH i(t)ni +p N  j=1 Rji1−(1 −α)ρH j(t)nj.(6) In the same fashion, the expression for IM i(t)inEq.(3) reads IM i(t)=1−1−λHM ieff i(t) neff iβ (7) where ieff i(t) is the number of infected humans placed in population iat time t ieff i(t)=(1 −αp)niρH i(t)+αp N  j=1 RjinjρH j(t).(8) The above equations describe the time evolution for the VBD incidence,  ρH(t)={ρH i(t)}and  ρM(t)={ρM i(t)},ina collection of connected patches with arbitrary demographic,  n, and vector,  m, distribution. Thus, by iterating Eqs. (2) and (3) starting from a given initial condition  ρH(0) and  ρM(0), we can monitor the spatiotemporal propagation of VBD and evaluate the steady epidemic prevalence on each geographical area. Let us first study the steady epidemic prevalence across patches  ρHas a function of the contagion probabilities λMH,λ HM, and the degree of human mobility pin the population. To reduce the number of parameters and without loss of generality, let us define λHM =λMH =λ.Westart analyzing the case in which human mobility is governed by an unweighted undirected Barabási-Albert network (BA) [67] of N=50 patches and average degree k=4, all of them homogeneously populated by ni=1000 agents. Concerning vector distribution, we consider that the ratio between vectors and humans populations inside a patch i, denoted in the following as γi, is randomly drawn from a uniform distribution within the range γi∈[0.3,1.7]. The top panels in Fig. 2 reveals the great agreement between the predictions obtained by iterating Eqs. (2)and (3) and the results obtained from mechanistic simulations for both the cases in which infected agents mobility is totally restrained, i.e., α=0 [Fig. 2(a)], and when there is no influence of the disease on agent mobility, α=1 [Fig. 2(b)]. Finally, we check if the formalism is able to capture the spatiotemporal unfolding of VBD. For this purpose, we start by setting a seed localized in a single patch of the BA network, and then monitor the temporal evolution of the fraction of the population affected by the disease inside each area. The bottom panels of Fig. 2show the evolution obtained by iterating Eqs. (2) and (3) [Fig. 2(c)] and by tracking the individual state of each agent in Monte Carlo simulations [Fig. 2(d)]. Again, it becomes clear that the formalism also reproduces the different propagation pathways of VBD in the synthetic metapopulation, despite the noise induced by the stochastic nature of the mechanistic simulations. IV. ESTIMATION OF THE EPIDEMIC THRESHOLD The validation of the Markovian formalism offers the possibility of saving computational costs by iterating 2 ×N equations instead of performing lengthy agent-based simulations. However, the advance behind Eqs. (2)–(8) is that they also allow for deriving metrics and analytical results about the dynamical behavior of VBD in large metapopulations. A relevant quantity that can be analyzed is the epidemic threshold, i.e., those conditions that turns the epidemic state into a stable solution of Eqs. (2)–(8). To address it, let us assume that the disease has reached the stationary state so ρH i(t+ 1) =ρH i(t)=ρH i,ρM i(t+1) =ρM i(t)=ρM i∀i. In this case, Eqs. (2) and (3) turn into μHρH i=1−ρH iIH i,(9) μMρM i=1−ρM iIM i.(10) As we are interested in computing the epidemic threshold, we now assume that the impact of the disease is finite but negligible for all geographical areas, which in mathematical terms implies that ρH i=H i1,ρM i=M i1∀i.Thisassumption allows us to linearize Eqs. (9) and (10) yielding: H i= N  j=1 λMHβ μHpRij mj ˜neff j +(1 −p)δij mi ˜neff i   Mij M j,(11) M i= N  j=1 λHMβ μMαpRji nj ˜neff i +(1 −αp)δij ni ˜neff i    Mij H j,(12) 013312-4 VECTOR-BORNE EPIDEMICS DRIVEN BY HUMAN … PHYSICAL REVIEW RESEARCH 2, 013312 (2020) 0 50 100 150 200 250 Time (days) 0 50 100 150 200 250 Time (days) 1 10 20 30 40 50 Patch ρH λλ 0.4 0.3 0.2 0.1 0 50 30 40 20 10 1250 200 15010050 0250 200 15010050 0 (a) (b) (c) (d) Time Time 0 .05 0 .1 .15 0 .2 .25 1 0.8 0.6 0.4 0.2 0 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.25 0.20 0.15 0.10 0.05 0 FIG. 2. Validation of the Markovian formalism. (Top) Epidemic size ρHas a function of the contagion rate between vectors and humans λHM =λMH =λand the human mobility p(color code). Dots correspond to results of averaging 50 stochastic realizations whereas solid lines represent the theoretical predictions obtained by iterating the equations of the formalism. Model parameters have been fixed to (β,μH,μ M)= (1,0.3,0.3). The values assumed for the restriction of the mobility of infected agents are (a) α=0and(b)α=1. (Bottom) Temporal evolution of the fraction of infected agents (color code) inside each patch obtained by iterating Eqs. (2)–(8) (c) and by tracking the individual state of each agent while performing a single Monte Carlo simulation (d). Infectivity is fixed to λ=0.3 whereas mobility parameter pis set to p=0.3. The rest of parameters are the same as in (a). where ˜neff ihas been defined as ˜neff i=neff i(0,α,p). Note that the form of the elements of these two matrices depends on both the mobility properties (p,R) and the demographic distribution of both agents and vectors ( n, m). For the sake of clarity, let us write the former system of equations in a more compact way:  H  M=0βλMH μHM βλHM μM M0 H  M.(13) Equation (13) makes evident the bipartite nature of the processes involved in the spread of VBD with matrices Mand  Mcapturing vector-to-human and human-to-vector infections, respectively. Thus, in order to quantify indirect infections between humans mediated by vectors and vice versa, we should iterate Eq. (13) obtaining  H  M=β2λMHλHM μMμHM M0 0 MM H  M.(14) From Eq. (14) it becomes clear that nontrivial solutions for  Hcorrespond to the eigenvectors of matrix M M. Specifically, given a metapopulation defined by  n, m,R, and p,the stationary solutions with infinitesimal incidence correspond to eigenvectors of M Mwhose eigenvalues can be written as i=μMμH β2λMHλHM .(15) Under these conditions, the maximum eigenvalue max(M M) encodes the combination of the RM parameters that corresponds to the epidemic threshold, namely, β2λMHλHM μMμHmax =1.(16) The former equation reveals the minimum infectivities, either λHM or λMH, that trigger the epidemic outbreak. To derive a simple critical infectivity one can set λMH =δλHM, so that λMH c=μHμM δβ2max(M M).(17) To test the validity of Eq. (17) we have carried out extensive numerical simulations in synthetic metapopulations considering λMH =λHM =λ. The top panels in Fig. 3show the epidemic diagrams ρH(p,λ) by computing the fractions 013312-5 DAVID SORIANO-PAÑOS et al. PHYSICAL REVIEW RESEARCH 2, 013312 (2020) FIG. 3. Epidemic threshold and evolution of the leading eigenvector of matrix M M. (a)–(c) in the top show three epidemic diagrams ρH(p,λ) in a synthetic metapopulation of N=50 patches. Each panel corresponds to a different rescaling value (α) for the mobility of infected humans, namely: (a) α=0, (b) α=0.5and(c)α=1. In addition, we have set λHM =λMH =λwhile the rest of the RM parameters are: μH=0.3, μM=0.3, and β=1.0. The color code show the incidence ρHas obtained from agent-based simulations while the solid curve represent the prediction for the epidemic threshold λc, calculated from Eq. (17). (d)–(f) show the evolution, as a function of p,oftheN components of the eigenvector of matrix M Mcorresponding to maximum eigenvalue max(M M). of humans infected ρHas a function of λand p. From these diagrams it becomes clear that, for each value of p, there exists a critical value λcso that for λ>λ cthe epidemic phase appears. The border of this region [solid curves in Figs. 3(a)–3(c)] is the function λc(p) calculated with Eq. (17), showing an excellent agreement with the results obtained from numerical simulations. V. ABRUPT TRANSITIONS OF LEADING PACTHES Apart from the agreement between Eq. (17) and the numerical simulations, the evolution of the epidemic threshold λc(p), reported in the three upper panels points out a nontrivial dependence with the degree of human mobility. Contrary to what naively expected, human mobility can be detrimental to epidemics, as clearly illustrated in the panels for α=0.5 and 1.0. This counterintuitive effect of mobility was already found for SIR and SIS diseases in networked metapopulations [46] as a result of the redistribution of the effective populations across patches due to mobility. In the case of VBD, this process corresponds to an homogeneization of the effective ratios between vectors and humans so that a high risk patch with large γi=mi/nitends to decrease its effective value due to the increase of neff icaused by human mobility. A more striking phenomenon reported in the epidemic diagrams of Fig. 3is revealed by the sharp variations in the slope of the curves λc(p). These abrupt changes are the product of collisions between the two maximum eigenvalues of matrix M Mas pvaries. This way, the two maximum eigenvalues interchange their order at some critical mobility value pc. These collisions do not have an strong impact in the epidemic threshold since the function λc(p) is continuous. However, they are the fingerprint of a sudden change in the form of the eigenvector corresponding to the maximum eigenvalue,  vmax, of matrix M M. This abrupt transition is of utmost importance since the components of  vmax encode the most important patches driving the unfolding of the epidemics. Let us recall that matrix M Mincorporates the demographic information,  n, the mobility patterns, Rand the vector distribution,  m, having as unique parameter the degree of mobility p. Thus, for each value of pthe spectral analysis of M M gives us the epidemic threshold λcand the distribution of patches triggering the epidemic onset in the components of  vmax. The evolution of the components of  vmax as a function of pis shown in the bottom panels, (d–f), of Fig. 3.From these plots, it becomes clear that the discontinuities of the slope of λc(p) correspond to abrupt changes in the form of  vmax. Namely, in the three cases, patch number 26 is the one causing the epidemic onset for p=0 and p1. This is obvious since patch 26 is the one with largest ratio γi=mi/niin the synthetic metapopulation. However, as p increases, the leading patch changes, being replaced by patch 18 in the case of α=0 (d) and α=0.5 (e), while for α=1 (f) the leading patch is replaced by a collection of them. Remarkably, the case α=0 shows a second abrupt transition at p c≃0.92. From a practical point of view, these abrupt changes point out that containment strategies targeting a certain neighborhood can sharply change from efficient to useless due to subtle changes in human mobility. In Appendix A, we confirm numerically these abrupt changes in the suitability of targeted immunization policies. 013312-6 VECTOR-BORNE EPIDEMICS DRIVEN BY HUMAN … PHYSICAL REVIEW RESEARCH 2, 013312 (2020) VI. EPIDEMIC RISK OF GEOGRAPHICAL AREAS Spurred by the ability of the Markovian formalism to capture the dynamics of VBD and the insights provided by the spectral properties of matrix M Mfor identifying the areas triggering the onset of epidemics, we move one step further and evaluate the epidemic risk associated with each patch. To this aim, we propose a theory-driven prevalence indicator which serves as a proxy to determine the most exposed areas to the spread of VBD. For thispurpose, let us analyze the elements of the matrices Mand  M, defined in Eqs. (11) and (12). From these equations we realize that the elements Mij ( Mij) contain all the possible microscopic contagion processes from vectors (humans) associated with patch jto humans (vectors) associated with i. Therefore it is possible to estimate the effective number of human-human interactions mediated by vectors that an individual from subpopulation ireceives from those with residence in j. This quantity, denoted as Cij, can be obtained as Cij = N  k=1 Mik  Mkj .(18) After introducing the expressions of Mjand ˜ Mkj from Eqs. (11) and (12) and after some algebra, the explicit expression for the indirect human-human contacts reads Cij =δij(1 −αp)(1 −p)mini ˜neff i2 +αp(1 −p)Rji minj ˜neff i2 +(1 −αp)pRij mjnj ˜neff j2 +αp2 k RikRjk mknj ˜neff k2.(19) This expression takes into account all the possible infection pathways connecting humans in patch jto an individual with residence in node i. In particular, those infections may take place in four possible ways [each one encoded by the terms of Eq. (19)]: (i) when an infected individual from patch j remains at his/her residence, transmits the disease to one vector in jwhich then transmits the disease to another resident in j; (ii) when an infected individual from patch jvisits patch i and infects a vector that will later pass the disease to a resident of patch i; (iii) when a resident of patch jinfects a vector there and a healthy human traveling to from ito jgets infected; and (iv) when both the infected individual from jand the healthy individual from itravel to a contiguous third patch kwhere the infection takes place mediated by a vector. Finally, in order to make predictions about the impact of the VBD on a geographical area, i, we must account for all the possible infections from each patch of the metapopulation and to weight the resulting number by the population of i. This way, the epidemic risk indicator for each patch i,inthe following denoted as ERi, can be defined as ERi=ni N  j=1 Cij =ni N  j=1 (M M)ij .(20) The evaluation of ERi, as defined above, can be done directly without the need of neither numerical simulations nor making the integration of the Markovian equations. In fact, once data about demography, vector distribution and mobility patterns are available, one can estimate the epidemic risk of each subpopulation. To validate the Epidemic Risk measure, we now move to a real metapopulation, the city of Santiago de Cali (Colombia). With a population of more than 2 millions of inhabitants, it offers the possibility of comparing our predictions in a scenario for which severe epidemic outbreaks of VBD are recurrently found. In particular, due to its location and climate, Cali is a Dengue endemic city in which records of the historical incidence of this disease are available for comparison. To this aim, we collected demographic and mobility datasets [68] whereas vectors abundance across districts was obtained from entomological reports made yearly by the local authorities [69] (further information in Appendix B). With this information at hand, and by using Eq. (20), we assign the epidemic risk of each of the 22 districts in which the city of Cali is divided. These epidemic risk values are compared to the observed Dengue incidence across the 22 patches during the years 2015 and 2016 [70,71]. In Fig. 4, we show this comparison by normalizing the values of both epidemic risk and Dengue incidence by their maximum observed value (in both cases that of district 13). In particular, we find a coefficient of determination of R2=0.81, indicating that the proposed prevalence indicator is able to capture the spatial distribution of Dengue cases across the city. On more general grounds, this agreement points out that given the demography, the commuting patterns and the spatial distribution of vectors across a given population, one can use Eq. (20) to identify areas where containment measures should be promoted to reduce the impact of possible outbreaks. VII. CONCLUSIONS The control of infectious diseases represents one of the major societal challenges. Understanding the complex interdependency between human activity and contagion processes is key to explain the onset and development of large-scale epidemics. Here, focusing on VBD, we have integrated information from urban daily commutes and the geographic distribution of humans and vectors to estimate the epidemic risk associated with different connected regions. In particular, we have provided a metapopulation formalism to assess the role that the former ingredients play on the propagation of VBD. We have proved that this formalism constitutes a very reliable and time-saving platform, since its Markovian equations enable to reproduce very accurately not only the global incidence of VBD but also the spatiotemporal spreading patterns observed in Monte Carlo simulations. Based on this agreement, we have derived an analytical expression of the epidemic threshold that captures the critical conditions which leads to the onset of epidemics. Apart from the detrimental effect that mobility may have on the spread of diseases, the study of the epidemic threshold has revealed interesting phenomena such as the existence of abrupt changes in the way epidemics unfold. In particular, we have shown that the subset of patches leading the epidemic onset can 013312-7 DAVID SORIANO-PAÑOS et al. PHYSICAL REVIEW RESEARCH 2, 013312 (2020) FIG. 4. Real Dengue incidence vs estimated epidemic risk in the city of Cali (Colombia). (Top) Normalized epidemic risk (ER/ERmax) vs normalized Dengue incidence (I/Imax) for each of the 22 districts of Cali. Color encodes the Epidemic Risk, from the lower (yellow) to the highest (blue).The correlation between the two variables yields a coefficient of determination of R2=0.81. (Bottom) Spatial distributions of the normalized Dengue incidence in the city of Cali (left) and the normalized epidemic risk (right) according to Eq. (20). The parameters concerning agent mobility have been set to (p,α)=(0.36,0.75). suddenly change as human mobility varies. This phenomenon highlights the need of incorporating real human mobility patterns into the design of containment policies targeting specific geographical areas, for efficient policies can turn useless due to a small variation of human mobility habits. Finally, relying on the matrix containing the information about the effective number of human-human contagions, we have derived an epidemic risk indicator that allows us to classify the patches according to their exposure to VBD. By computing this epidemic indicator, we have reproduced with great accuracy the geographical distribution of Dengue incidence in the city of Santiago de Cali (Colombia), where Dengue in an endemic disease, thus being able to identify the most vulnerable areas where prevention measures should be promoted. In a nutshell, our results point out that the spread of VBD is the result of a delicate interplay between commuting flows, human census and vector distribution. This interplay is captured both in the analytical expression of the epidemic threshold and in the epidemic risk indicator. As a result, we have shown that small variations of the former ingredients, such as the degree of mobility, can lead to abrupt changes in the way epidemics unfold. Our framework, although containing several simplifying assumptions to allow the analytical treatment, has shown useful to integrate human and contagion dynamics and it can be readily implemented to identify those regions where immunization policiesshould be reinforced and to forecast the consequences of control strategies focused on mobility restrictions. As future work, the formalism here presented paves the way to the incorporation of additional features that, together with human mobility, are key for the dissemination of vectorborne diseases. Among these features we find (a) the inclusion of seasonality modulating both, the population and the bitting rate of vectors [72] or the effects of climate change on vectors census [6–9], both ingredients can be incorporated as an external driver for each patch through the use of the multiplex formalism [73–75]; (b) the time-varying nature of human contact patterns [76–79] and the adaptability of mobility behaviors in epidemic scenarios [80,81]; and (c) the possibility of incorporating more refined compartmental models with the addition of multiple states [82,83] to address the impact that direct human-human transmission in some VBD such as Zika has on the epidemic threshold [62], or the interplay of several VBD with competitive or cooperative interactions [84–86]. These ingredients can be readily incorporated into the Markovian model, improving its capacity to evaluate the risk of specific VBD in particular regions. ACKNOWLEDGMENTS We are especially grateful to A. Arenas and O. Vasilieva for useful comments and discussions. This work was partially supported by Universidad del Valle (grant CI-165), Gobierno de Aragon/Fondo Social Europeo (Grant No. E36-17R), Ministerio de Economia, Industria y Competitividad (MINECO) and Fondo Europeo de Desarrollo Regional (FEDER) (Grants No. FIS2015-71582-C2 and No. FIS2017-87519-P) and by the Spanish Ministry of Science, Innovation and Universities, the National Agency for Research Funding AEI and FEDER (EU) under the grant PACSS (RTI2018-093732-B-C22) and the Maria de Maeztu program for Units of Excellence in R&D (MDM-2017-0711). APPENDIX A: EFFICIENCY OF TARGETED IMMUNIZATION POLICIES The study of the critical properties of VBD in the main text has revealed the existence of some mobility values, denoted as pc, for which the components of the leading eigenvector change abruptly. Translated into epidemiological words, this phenomenon pinpoints the change in the most affected patches, which is of great relevance since targeted policies in specific areas can pass from useful to useless as human mobility varies. To prove it, we now study the effects of applying prevention measures in specific locations selected according to the largest components of the leading eigenvector of the critical matrix M˜ M. We consider the synthetic network used in the main text and set α=0.5, for which the change in the leading patch happens at p=pc=0.38 (see Fig. 3). To study the 013312-8 VECTOR-BORNE EPIDEMICS DRIVEN BY HUMAN … PHYSICAL REVIEW RESEARCH 2, 013312 (2020) FIG. 5. Temporal evolution, according to the Markovian equations, of the number of infected agents by a VBD which spreads over a BA network. The black line denotes the original curve where no policy has been implemented whereas dashed lines correspond to the case in which patches 26 (red) and 18 (blue) are immunized. The epidemic parameters have been set to (α, λHM,λ MH,μ H,μ M)=(0.5,0.25,0.25,0.3,0.3). The values for human mobility are (a) p=0.1(p<pc)and(b)p=0.5(p>pc). effects of containment measures, we assume that immunized agents are no longer susceptible of contracting the disease, so λMH =0 for them. Finally, we target two different subpopulations for the immunization policies: patch 26 that is the leading patch for p<pc, and patch 18 that sustains the epidemic outbreak when p>pc. Figure 5confirms that the effectiveness of targeted policies against VBD is strongly influenced by aspects concerning human mobility. In fact, immunizing agents from patch 18 leads to the extinction of the disease for p=0.5>pc, whereas it is almost ineffective for mobility values below this threshold (p=0.1<pc). This result, along with the others presented in the main text, highlights the importance of designing containment policies taking into account the complex interplay between human mobility patterns, census data and vector abundance. APPENDIX B: DESCRIPTION OF CALI DATASET One of the most important contributions of this work is the formulation of a new framework which can easily incorporate mobility data of real cities to address real epidemic scenarios. In the main text, we tackle the spread of VBD in the city of Cali (Colombia), whose geographical and meteorological features make it an endemic region for several VBD such as Dengue, Chikungunya or Zika. In particular, here we focus on the spread of Dengue. To assess the effect of mobility on the spread of Dengue in Cali, it is necessary to reconstruct the mobility network of its inhabitants from data. For this purpose, we divide the city into 22 districts, which correspond to the official administrative divisions called comunas. Regarding demography, the population distribution across comunas has been extracted from census data that the municipality facilitates [68]. Mobility flows connecting comunas are extracted from urban commuting surveys [87]. As a result, more than 105 trajectories were recorded, which suppose a representative sample of Cali’s commuting flows. Once all the data have been gathered, an origin-destination matrix, encoded in our formalism by matrix R, is computed as Rij =Wij N l=1Wil ,(B1) where the numerator corresponds to the number of trips between patches iand jwhile the denominator counts all the reported trips departing from patch i. The result is a weighted directed network encoding the probability that an agent visits other neighborhoods different from its residence. Apart from the mobility network, thedistribution of vectors across the city also plays a crucial role in the outcome of the disease. The number of vectors inside a geographical region is strongly linked to environmental features such as altitude, temperature and humidity but also to human-dependent factors like health and economic conditions. To model vector distribution across comunas, we use as a proxy the socalled recipient index. This quantity encodes the probability of finding vector pupae in different recipients which have been previously distributed across the city. A high value of the index means a higher probability of finding vectors. For this reason, we assume that the ratio between the number of vectors and humans inside each patch in our model is directly proportional to its recipient index, which is extracted from entomological data of the year 2015 [69]. [1] Zika situation report, World Health Organization Tech. Rep. (2016). [2] A global brief on vector-borne diseases, World Health Organization (2014). 013312-9