Critical behavior in interdependent spatial spreading processes with distinct characteristic time scales

The spread of an infectious disease is well approximated by metapopulation networks connected by human mobility flow and upon which an epidemiological model is defined. In order to account for travel restrictions or cancellation we introduce a model with a parameter that explicitly indicates the ratio between the time scales of the intervening processes. We study the critical properties of the epidemic process and its dependence on such a parameter. We find that the critical threshold separating the absorbing state from the active state depends on the scale parameter and exhibits a critical behavior itself: a metacritical point – a critical value in the curve of critical points – reflected in the behavior of the attack rate measured for a wide range of empirical metapopulation systems. Our results have potential policy implications, since they establish a non-trivial critical behavior between temporal scales of reaction (epidemic spread) and diffusion (human mobility) processes. The time scales of individual contagion and human mobility are both relevant to understand the evolution of epidemics. Here, the authors study a reaction-diffusion process for an epidemic metapopulation model, showing that modulating the relative time scale of these two processes produces different epidemic outcomes and can reproduce the effect of non-pharmaceutical interventions.

T he spatio-temporal spread of infectious diseases is inherently interdependent with human behavior [1][2][3][4] . On the one hand, individuals move between spatial patches, effectively creating human flows through which pathogens travels from one place to another. On the other hand, social dynamics can favor or hinder potential transmission channels: a typical example for the latter relies on risk perception and spread of awareness 3,[5][6][7] , determining the choice to take individual actions-such as social distancing and wearing masks-or mid-scale non-pharmaceutical interventions-such as school or workplace closure and household quarantine-and have the direct effect of reducing the probability that a pathogen is transmitted from one infectious individual to susceptible ones in their social neighborhood. Convincing evidence for the role of human behavior in the spreading of several infectious diseases has been reported in the last decades for Foot-and-mouth disease 8 , Influenza and its strains [9][10][11][12][13] , SARS 14,15 , Ebola 16 , Dengue 17,18 , Zika 19 and, very recently, for COVID-19 [20][21][22][23][24] .
Consequently, a huge effort has been devoted to understand the fundamental dynamical mechanisms where it is possible to take action through policy in order to drive a system subjected to an ongoing epidemics toward a controlled state, rather than an uncontrolled one. A standard approach is to coarse-grain the system under investigation in terms of spatial patches, named metapopulations, which account for groups of individuals localized at a given geographic scale. Metapopulations can be thought as nodes of a complex network of spatial patches, where links encode human flows from one place to another and are responsible for between-patch transmission 25 . In fact, a metapopulation can be one neighborhood within a city 26,27 , a city 28,29 , a subnational region or higher scales [30][31][32][33][34] . Models for human mobility can range from Markov chains encoding standard diffusion 35 to more sophisticated dynamics encoding higher-order memory 36,37 : they are often encoded into origin-destination (OD) matrices which are later encoded into the essential diffusive part in the dynamical equations that describe the spatio-temporal evolution of epidemics in terms of reaction-diffusion processes.
However, quite often those models assume that an epidemic unfolds slower than the average time people need to move around 38 . This assumption, known as separation of scale approximation, is verified in practice, since the incubation period and the average duration of many diseases is of the order of days, which is much longer than the typical time scale of human movements, ranging from about 30 min for urban trips to 8 h for interurban ones 39,40 . Nevertheless, there are specific scenarios which can alter the typical time scales of the intervening dynamical processes, e.g., nonpharmaceutical interventions employed for disease containment can dramatically change the time scales of social dynamics. Specifically, changes in human mobility can be due to extraordinary travel restrictions, cancellations or to the effects of other restrictions such as school and workplace closure or delayed access to point of interests, from restaurants and gyms to shopping malls and museums, as happened during the COVID-19 pandemic [41][42][43][44] .
In this work, we propose a model able to cope with this scenarios by explicitly accounting for the existence of two distinct time scales, each one characterizing the intervening reaction and diffusion processes. The time-scale separation is therefore encoded into a parameter ϵ, which is defined in terms of the ratio between the diffusion and reaction time scales. This parameter allows us to tune the behavior of a system between two extremal scenarios-such as reaction faster diffusion (ϵ > 1, Fig. 1a) and diffusion faster than reaction (ϵ < 1, Fig. 1b)-and all the intermediate ones. For simplicity, in the following we will refer to this model as the ϵ-model. To allow for an analytical treatment of the problem, we focus our attention on the emblematic Susceptible-Infected-Recovered (SIR) epidemic model for a metapopulation network, where particles flow at each time step according to transition rules encoded by empirical human flows. Using meanfield approximations and perturbative analysis, we find an analytic expression for the basic reproduction number R 0 (ϵ), which indicates if the epidemics will reach a finite fraction of the population (R 0 (ϵ) > 1) or to die out (R 0 (ϵ) < 1) in the stationary state, effectively separating the active state from the absorbing state, respectively. We then find that the co-evolution of reaction and diffusion on different time scales reveals the existence of a phase transition characterizing the behavior of the critical threshold R 0 (ϵ) as a function of the parameter ϵ: the value ϵ ⋆ separates the regime where R 0 (ϵ) is constant (ϵ < ϵ ⋆ ) from the regime where it grows monotonically (ϵ > ϵ ⋆ ). Note that for this behavior to be observed it is required that transmission rates are distributed among nodes in a nonuniform way. We identify a quantity R min 0 that parametrizes such nonuniform distributions so that the point ðR min 0 ; ϵÞ ¼ ð1; ϵ ? Þ is a metacritical point, which separates all epidemic outbreaks between certain or possible based on the separation of time scales, the infectiousness distribution and the topology of the network.

Results and discussion
Modeling of epidemics through multiscale reaction-diffusion processes. Let us consider the scenario of M geographical patches where individuals move from patch i to patch j (i, j = 1, 2, …, M) with probability rate P ij , encoding the entries of the corresponding mobility matrix P. Let us also consider a standard SIR model for the spreading of an infectious disease, subjected to the conditions , and R i (t) are the number of susceptible, infected and recovered agents in the i-th patch at time t, while N i (t) is the overall population of the i-th patch at time t, with ∑ i N i = N the total population, assumed to be constant over time. The variation of the number of infected individuals in each patch due to the epidemic spreading is given by where β i encodes the transmission rate in patch i for a communicable disease with empirical recovery rate given by μ. The possibility of having patch-dependent transmission rates is plausible for a metapopulation network where each spatial patch is characterized by a population with socioeconomic or clinical differences, such as average age or income. During the COVID-19 pandemic, for instance, such differences were reflected in the variance of infection fatality rates 45 . If the population is allowed to move among patches according to the mobility matrix, the variation due to diffusion is given by for infected individuals. Since we are working under the assumption that the epidemiological state of an agent does not affect their diffusive behavior, the equations for susceptible and recovered individuals have the exact same form. The time scales Δt epi and Δt mob should be interpreted as the minimum time necessary to have a certain small variation in the number of infected individuals, because of the epidemic process or the mobility one, respectively. Usually, those two distinct dynamics happen at the same time but, in principle, the corresponding time scales of epidemics spreading and mobility dynamics are not necessarily the same. This richness of dynamic scenarios can be encoded into a tunable parameter ϵ that slows down or speeds up diffusion with respect to epidemic spreading. We define it as the ratio ϵ = Δt mob /Δt epi between the two time scales, so that the differential equations for the model become dI i ðtÞ dt ¼ β i S i ðtÞ I i ðtÞ N i ðtÞ À μI i ðtÞ see Methods for a more explicit derivation of this equation. In the remaining of this paper, we will refer to this model as the ϵmodel, since the presence of the ϵ parameter, characterizing the interplay between diffusive and reactive time scales, is its most distinctive feature. At this point, it is worth remarking that we can redefine both P and ϵ without loss of generality such that ∑ M j¼1 P ij ¼ 1 and all the information on the time scale of the mobility lies in 1/ϵ. The reason for doing this will become clear in the next section.
Note that having a parameter that multiplies the mobility term in a metapopulation model is not uncommon 29,35 . In these cases such a parameter indicates the fraction of people that move between patches in a given time step. This is exactly equivalent to our interpretation of ϵ as the ratio between the time scales of the two processes of reaction and diffusion.
Calculation of the R 0 . In a simple SIR model without spatial structure, where every agent interacts with every other agent and there is only one transmission rate β for all spatial patches, it can be shown that the basic reproduction number is simply R 0 = β/μ. However, in our framework we assume that there is a β i for each patch i and we expect to have some dependence also on the mobility matrix and on ϵ. To this aim, we use the next-generation matrix formalism proposed by van den Driessche et al. 46 to calculate the critical threshold R 0 (ϵ).
First, it is convenient to rewrite Eq. (3) in terms of the concentration of infectious and susceptible individuals, indicated respectively as ρ i (t) = I i (t)/N i (t) and σ i (t) = S i (t)/N i (t). For simplicity, it is also useful to write the differential equation in matrix form as where L > ¼ I À P > and B ¼ diag ðβ 1 ; ; β n Þ. Let us consider the infection-free steady state, i.e., the case right above the critical point, meaning that there is a negligible fraction of infected individuals (σ i = 1, ∀i). The resulting system of equations reduces to _ ρðtÞ ¼ JρðtÞ; where J ¼ B À μI À ϵ À1 L > is the Jacobian matrix, whose analysis allows one to gain insights about the linear stability of the original system dynamics. By decomposing the matrix into J ¼ 1 ϵ ðF À VÞ, it is possible to demonstrate that the highest eigenvalue of the so-called next-generation matrix K = FV −1 is in fact the basic reproduction number R 0 46,47 . In our case F ¼ ϵB However, solving the eigenvalue problem for this M × M matrix is, in general, a very complicated task and tells you very little about the general dependences of K on the parameters of the problem. To gain some insights, we look for approximated solutions in the extreme cases in which mobility is either much slower or much faster than the epidemic, in order to characterize the behavior of the system from an analytical perspective. represented with dots whose colors indicate their state with respect to the disease. The big circles represent different spatial patches and the probability of going from patch j to i is represented with the black arrows and encoded in the matrix element P ij . The agents can either move around in the network, and in that case they are represented along the arrows outside the spatial patches, or stay where they are and undergo one of two epidemiological "reactionsˮ: infection or recovery. The infections are represented with the transformations with red arrows, while for the recoveries we used blue arrows. Finally according to the time scale separation parameter ϵ we can distinguish two regimes of this model: in (a) ϵ > 1, so reaction events are more frequent that diffusive ones while in (b) ϵ < 1, so the opposite is true.
Diffusion slower than reaction. In this first case we have that ϵ ≫ 1, so we can use perturbation theory (PT) using ϵ −1 as our parameter and writing the next-generation as The eigenvalues and eigenvectors of the unperturbed part are respectively β i /μ and the canonical basis {e i }. The resulting spectral radius, quantifying R 0 , is given by where the index max indicates the patch with the highest transmission rate (see Methods). Equation (8) highlights that when there is no mobility at all (ϵ → ∞) one obtains R 0 ¼ β max =μ: i.e., the condition for the epidemic to have a finite-size outbreak is for that to happen at least in one spatial patch, as expected. It should be noted that in this case we have a strong dependence of the systems on the initial distribution of infected individuals, since the epidemic cannot grow unless the patches with β i > μ have at least one infected individual at time t = 0. Furthermore, notice how the mobility network P only appears in the second order term and therefore its effects become more and more negligible as ϵ increases.
Diffusion faster than reaction. When the mobility is much faster than the epidemic (meaning that ϵ ≪ 1) we cannot use the same approach, since for ϵ → 0 the next-generation matrix becomes illdefined as it is defined in terms of the inverse of the Laplacian matrix L, which is singular. Nevertheless, one can solve this problem by using a corollary to Perron-Frobenious theorem which states that if one has a non-negative irreducible matrix, such as K, then its highest eigenvalue corresponds to the spectral radius ρ(K) and the following inequalities hold: It is possible to show (see Methods) that in the limit ϵ → 0 the next-generation matrix can be found exactly, leading to A closer inspection of the above results reveals that in this scenario R 0 is the weighted average of the local basic reproduction numbers-defined by β i /μ in the case of isolated spatial patcheswhere the i-th weight is given by the relative strength s i of node i, i.e., the sum of all the flows involving that node, with respect to the overall human flow. It is straightforward to verify that the above relation reduce to the expected R 0 = β/μ when β i = β for all spatial patches. It is also worth noticing that Eq. (10) tells us that R min 0 does not depend on degree correlations or on topology in general and, as a result, as ϵ vanishes, every network with the same strength sequence leads exactly to the same result. That also means that, in the case of unweighted metapopulation networks, the configuration model constitutes a good approximation with respect to R 0 when ϵ → 0. The notation R min 0 encodes the fact that, since R 0 (ϵ) is a monotonically increasing function, the basic reproduction number for ϵ → 0 corresponds to the lower bound (see Fig. 2).
Diffusion comparable with reaction. It is plausible to wonder about what happens in the intermediate regime ϵ ≃ 1. Figure 2 shows that R 0 (ϵ) can be written as the piecewise-definite function where R min 0 is given by Eq. (10) and a, b and ϵ ⋆ are three parameters, one of which is fixed by the continuity condition of f(ϵ) in ϵ = ϵ ⋆ . We can then safely assess that R 0 (ϵ) displays itself a phase transition, and as such it is possible to distinguish between the ordered and disordered phase just by looking at the initial condition in the number of infectious individuals in each patch. For instance, for ϵ < ϵ ⋆ the diffusion is so fast compared to the reaction that no matter what the initial conditions are in the longtime limit the number of infected individuals will be uniformly spread among all the patches; therefore this phase is ordered or, equivalently, invariant under the choice of initial conditions. Conversely, when ϵ > ϵ ⋆ the long-time limit will strongly depend Barabási-Albert metapopulation network, models for homogeneous and scale-free connectivity of spatial patches, respectively. The solid gray curve indicates the values of R 0 found numerically by solving the eigenvalue problem for the next-generation matrix K, whereas the red line corresponds to the value analytically derived in Eq. (8) by using perturbation theory. The approximation is excellent for the homogeneous connectivity, whereas in the case of the heterogenous one it works only for ϵ > 1. The vertical dashed violet line indicates the position of the metacritical point ϵ ⋆ obtained analytically, that separates the regime in which R 0 stays constant and equal to its minimum value R min 0 (see Eq. (10)) from the regime in which it start growing. Finally the horizontal dotted black lines correspond to the lower bound on R 0 , which in this case is R min 0 ¼ 1, and the upper bound β max =μ ¼ 2, where μ = 2 is the recovery rate and β max ¼ 4 is the maximum transmission rate among those in the network.
on initial conditions, therefore the aforementioned invariance is broken and the phase is disordered. This behavior may be seen as analogous to the so called "global invasion thresholdˮ 48 , which indicates the condition under which the epidemic goes from being localized only in a few geographical patches to invade the whole network.
To estimate the location of the critical point in the parameter space we can either obtain it from fitting Eq. (11)-as we have done for drawing the dotted lines in Fig. 2-or look for the interception between the sigmoid in Eq. (8) and the horizontal line R 0 ¼ R min 0 . In particular for R min 0 ¼ 1 and its corresponding ϵ ⋆ we get a metacritical point: i.e., a value of the critical threshold that divides absorbing states from active ones (see Fig. 3). In this case we have ϵ ? ¼ 1=ðβ max À μÞ, obtained by intersecting the line R 0 = 1 with the first term in Eq. (8) while assuming that there are no self-loops. We obtain a really good estimate for synthetic mobility networks characterized by homogeneous connectivity and our approximation still works as an upper bound for the ones characterized by heterogeneous connectivity. This is true because, as shown in Fig. 2, in heterogeneous networks the second order term in PT starts being relevant, and since this term can only be positive the net effect is to lowering the metacritical point. Therefore, we can conclude that the maximum possible value for R 0 should be the one obtained only with the first order in PT, that is very close to that of homogeneous networks (see again Fig. 2). Finally, note that if Eq. (11) is correct (as it is for homogeneous networks) then the critical exponent ν given by R 0 $ ðϵ À ϵ ? Þ ν for ϵ → ϵ ⋆ , is equal to 1.
Metacritical point observed in real human flows. We now apply our theoretical framework to analyze empirical metapopulation networks, where nodes are either counties of a U.S. state or other types of administrative regions in other countries, and where links are observed human mobility flows.
We use a series of datasets describing the anonymized mobility flows in different areas of the world, aggregated at administrative levels. All sources originated by the usage of mobile phones: this kind of passively collected data has been successfully validated for epidemics modeling 49 and became standard in recent years. The flows for France (FRA), Spain (ESP) and Portugal (PRT) 49 , Ivory Coast (CIT) 7 and Turkey (TUR) 33 are reconstructed from the Call Detail Records, the billing dataset held by mobile phone companies. The flows for USA (Florida and Colorado) are estimated on the base of data kindly provided by Cuebiq, a location intelligence and measurement platform that collects privacy enhanced GPS trajectories from mobile apps users who have opted-in to provide access to their aggregated location data anonymously. Also those of Italy (ITA) are derived from Cuebiq data, and are obtained from a published dataset 50 . These flows have been used for assessing change in mobility during the COVID-19 pandemics in 2020. 51,52 , but for our purposes we analyze here only a static snapshot of the mobility in the area. See Fig. 4(a), (b) and (c) for an illustration of the mobility flows in Italy, Florida, and Colorado. For these datasets, the population of each node is obtained from census data.
The Philippines. In this case, the population at each site is assumed to be proportional to the total out-bound flow, upscaled so that the sum of the node populations matches the country population.
In order to verify the validity of our theoretical framework when dealing with real-world mobility flows, we assign to each node a transmission rate β i which is proportional to its population. We choose to do this because it is reasonable to think that more populated areas are linked to a larger number of interaction and therefore to a higher risk of transmission. However we cannot expect this mechanism to hold for arbitrarily large populations, since the number of contact each person can have has a limit that does not depend on the size of the population itself (we got the idea behind this from 54,55 ). For this reason we set the maximum value to β max ¼ 2μ. In addition we set R min 0 ¼ 0:5 to be sure that we are in the regime where a transition happens for some value of ϵ (see Fig. 3). Each simulation was run until the variation in the number of recovered individuals after 20 days was 0.01% of the total population (note that for China and India this fraction was lowered to 0.0001% due to the large number of people living inside each node). For each of simulation, initial conditions were set to ten infected individuals per node. Finally the heterogeneity of every weighted network was calculated using the Gini coefficient G of the strength distribution of the network.
The results are shown in Fig. 5, highlighting the behavior of the attack rate at different stages of the epidemic, for various realword mobility networks. Here with epidemic stages we simply mean the status reached by the epidemic after different times, measured in multiples of the time necessary to reach the maximum number of infected (or infectious peak). We chose to do this so that, even if the absolute duration of each simulation had varied, we would have been at the same point (or "stageˮ) of the evolution of the epidemics. The results display, qualitatively, a similar behavior: for low values of ϵ the attack rate is approximately zero, meaning that there is no epidemic outbreak at all, while for larger ϵ this number grows rapidly, confirming the existence of a phase transition in the curve of critical thresholds R 0 (ϵ). Notice that such a transition point ϵ ⋆ always remains We recall that ϵ is the parameter that indicates the ratio between the time scale of the mobility and that of the epidemic, while R min 0 is the weighted average of the local basic reproduction numbers-defined in the case of isolated spatial patches-using the patches' strengths as weights (Eq. (10)). The dashed lines correspond to approximately below the upper bound 1=ðβ max À μÞ, which is equal to 10 in this case. Finally, when ϵ is very large there is a drop in the attack rate, which happens because inter-node travels are so limited that the epidemic can only spread inside the nodes where β i > μ, quickly reaching the stationary state.
Note that here, when taking the long-time limit or stationary state, we are actually referring to the stage in which the epidemic has slow down so much that variations are extremely small. Curves in Fig. 5 weakly depend on such a condition, even though their qualitative shape would remain the same: an abrupt growth for low values of ϵ followed by a sharp drop for high values of ϵ.
Finally, from a quantitative perspective, the different behavior of the attack rate across countries depends mainly on the population distribution at their sites and on the degree distribution of the corresponding metapopulation networks. In particular, countries with a more homogeneous degree distribution (such as Zambia or Philippines) tend to have a transition point ϵ ⋆ closer to the upper bound 1=ðβ max À μÞ, while for heterogeneous networks (such as Portugal or Turkey) this value is lower, as confirmed by Fig. 2. Then, as ϵ increases, the population distribution becomes more and more important, since in this limit people cannot move so much and they are bound to stay in their original patch. In this case, the outbreak will spread only in those sites where the local R 0 is larger than one. Since we have defined transition rates to be proportional to the local population, the nodes with the highest β i are also the most populated ones. Therefore, if we have a network with a very heavy-tailed population distributions (such as Turkey) the attack rate in the limit ϵ → ∞ will be large while, conversely, in the case of a network with a more homogeneous population distribution (such as Nigeria) this limit is low.

Conclusion
In this work, we introduced a model for epidemic spreading on metapopulation network with arbitrary topology. It is structured as a reaction-diffusion process where the reaction part is that of a standard SIR model with site-dependent transition rates, while the diffusive part corresponds to a random walk whose dynamics is scaled by a parameter which allows us to tune the speed of mobility with respect to the epidemic spread. This framework is more convenient than classical single-scale reaction-diffusion models when conditions such as non-pharmaceutical interventions and behavioral changes become non-standard, altering the characteristic temporal scales of human mobility.
We provided analytic predictions for the basic reproduction number R 0 in two extreme cases, namely the "diffusion much faster than reactionˮ regime, thanks to a perturbative approach, and the "diffusion much slower than reactionˮ regime, where we make use of the Perron-Frobenius theorem. In the latter case we discovered that the global R 0 can be well approximated by the weighted average of the local R 0 with the patches' degrees as their weights.
In the intermediate regime, the most interesting one, we discovered the existence of a metacritical point on synthetic networks and on real human mobility flows. Such a point separates two distinct regimes observed in the behavior of the critical threshold, identifying the region of the parameter space where the conditions for having a non-vanishing attack rate change abruptly, providing useful information on how the different interventions to control or mitigate the epidemic in each site work in combination with travel restrictions.
To the best of our knowledge this type critical behavior in the curve of critical points has been observed in the presence of intertwined dynamical processes such as awareness and epidemics 6 -where the spreading of awareness about a disease is able to inhibit the epidemic spreading-and in competing epidemics spreading 56 -where the spreading of a disease is able to facilitate or inhibit infection by another disease-on multilayer networks 57 , but it is the first time that it is observed for multiscale processes in classical networked systems.

Methods
Derivation of the ϵ-model. Here we explain how to go from Eq. (1) and (2) to Eq. (3). First of all we recall that Δt epi and Δt mob are the time intervals that we have to wait in order to get a variation in the number of infected individuals equal to some arbitrarily fixed value. This value however should be small, because the transition rates in both Eq. (1) and (2) depend on time, so waiting too long would mean to change the value of the transition rates. For instance we could decide that Δt epi and Δt mob are those that you need to have ΔI epi i and ΔI mob i equal to 1, respectively. Now we ask ourselves: at each time step how much of the the variation in infected individuals is due to epidemics and how much to mobility? In order to answer that we compute the joint variation in infected individuals after an arbitrary time Δt From this equation one can see that the contribution to ΔI epi i due to the epidemics is always a factor 1/ϵ larger (or smaller) than the contribution ΔI mob i due to the mobility, which justify the presence of the 1/ϵ factor in Eq. (3). Finally if we incorporate the 1/Δt epi factor in the definitions of the rates and then send Δt to zero we end up with the differential equation of the ϵ-model.
Perturbation theory. First of all we recall that the general form of the nextgeneration matrix for the ϵ-model is K ¼ ϵB ϵμI þ L > À Á À1 and that we can write it us a sum K ¼ K 0 þ ϵ À1 K 0 where the terms are defined in Eq. (7). We also recall that the aim of PT is that of writing the eigenvalues as a power series in the parameter ϵ −1 . In our case it will be sufficient to consider only the expansion up to second order, that is The unperturbed eigenvalues are simply λ ð0Þ i ¼ β i =μ, which are the "localˮ basic reproduction numbers, meaning the values that R 0 would assume inside each patch if they were all isolated, while the unperturbed eigenvectors are those of the With epidemic stages we simply mean the stages of development of the epidemic measured using as time unit the time necessary for an epidemic to reach the maximum number of infected (or infectious peak). In all cases, the critical point ϵ ⋆ is clearly visible, since it corresponds to the value of ϵ such that the attack rate is larger than zero. It is worth noticing that such a critical point is approximately always lower than our analytically found upper bound 1=ðβ max À μÞ represented with a vertical dashed black line (with β max ¼ 0:2 and μ = 0.1 in this case). The G in each sub-figure indicates the Gini coefficient of the strength distribution of the corresponding network and it is used to quantify the heterogeneity of each network. Finally, we notice a drop in the attack rate for large values of ϵ: an effect due to the slowing down of the mobility.
canonical basis {e i }. On the other hand the first order contribution is Putting together the first two terms of the power expansion (Eq. (13)) we get where the max index points to the node with the highest transmission rate β max . For the second order in PT the key factor to compute is Putting this in the formula for the second order contribution and adding it to Eq. (15) we obtain Eq. (7).
The ϵ → 0 limit. We now approach the problem of computing the next-generation matrix K = FV −1 for an arbitrary mobility network. We proceed step by step by computing one piece at the time, starting with the elements of V ¼ ϵμI þ L > . This matrix in general is not symmetric, meaning that the right and left eigenvalues are different where φ (n) and ψ (n) are respectively the right and left eigenvectors of both V and L > , while the eigenvalues ν (n) of V differ from the eigenvalues λ (n) of L > by a term ϵμ. Even though we do not know neither the eigenvectors nor the eigenvalues it is useful to write the matrix in the following form Finally if the matrix V is invertible, and we know it is thanks to the term ϵμI, we can write the first building block of the next-generation matrix, that is The second step consists in multiplying what we have with the matrix ϵB and then take the limit for ϵ → 0. This will lead to where first we wrote the eigenvalues of V in function of those of L > and then we noticed that the only term in the sum over n who was not multiplied by the the vanishing ϵ was the one associated to the null eigenvalue of the Laplacian matrix λ (1) = 0, that corresponds to ν (1) = ϵμ. Fortunately for us we know that the eigenvectors corresponding to the stationary states are where s i is the strength of the i-th node. Finally we can apply the corollary to Perron-Frobenius theorem mentioned in the main article, but instead of doing it on K we do it on K ⊤ which has the same spectrum and the property of having all the rows equal to each other, so that upper and lower bound coincide and we get that the spectral radius of K is exactly that written in Eq. (10).
Mobility networks data. Our analysis of real human flows are based on (i) mobility flow networks and (ii) local population that were obtained from a broad range of different data sources. We have in general two types of flows: estimated and empirical. In the estimated category fall the flows for Senegal (SEN), Nigeria (NGA), Zambia (ZMB), China (CHN), India (IND), Mexico (MEX), Brazil (BRA) and the Philippines (PHL). In this case the data are therefore not representing real flows but the estimated internal human migration flows in a series of countries where malaria is endemic. The reconstruction is the product of the complex datafusion of census-based migration microdata, which are assimilated through gravity modeling. For further detail please refer to the Data Descriptor 53 . Being the local population resident in the nodes not part of the data shared, we used as a proxy the total out-bound flow, which we upscaled so that the sum of all nodes populations was equal to the total population of the country. All other datasets represent instead empirical flows and are associated with node population estimated from census data. The FRA, ESP and PRT 49 , CIT 7 and TUR 33 are reconstructed from Call Details Records, which provide information on the position (closest antenna) and time for each phone call, texting, or internet access of mobile users, which are recorded by the mobile companies for billing purposes. Please refer to the cited papers for the data availability of the respective dataset. The flows for USA (Florida and Colorado) are estimated on the basis of data kindly provided by Cuebiq, a location intelligence and measurement platform that collects privacy enhanced GPS trajectories from mobile apps users who have opted-in to provide access to their aggregated location data anonymously. Cuebiq provided us with information, on a weekly basis, about: (i) the ratio of users who have been seen moving between county i and county j to the number of users seen in county i ratioZ ij = (users flow) ij /(county users) i (ii) the ratio of devices that we see in county i to all the devices seen across all counties U i = (county users) i /(total users). We therefore upscaled to the values of flows one would have if the movements observed were to represent the entire population with the formula where Pop is the total population across the country considered here, the United States of America. For this analysis, we selected only the flows associated to the first week of 2020, so before the diffusion of the COVID-19 epidemics in the USA. These data are available from the authors upon request. Finally, also the flows of Italy (ITA) are derived from Cuebiq data, which in this case are again publicly available and we invite to refer to the Data Descriptor 50 for more details on the methodology. The italian flows shared capture the fraction of the total userbase moving between two italian provinces. We upscaled to the total population by simply multiplying by the population data obtained from census. In this case, we analyzed the flows associated with the 18th January 2020, before the beginning of the COVID-19 pandemic. The fundamental descriptors of these networks can be found in the in Table 1. The bold font indicates the type of information contained in the corresponding column. Here the word "Giniˮ refers to the Gini coefficient of the strength distribution of the corresponding network, while the symbol # stands for the word "numberˮ.
Simulation methodology. Here we put all the information on the Figs. 2 and 3 that we could not fit in the main text. For instance in Fig. 2 both the networks had n = 500 nodes. The Erdős-Rényi network was undirected and sampled using a wiring probability p = 0.04; the Barabási-Albert network was also undirected and sampled adding at each time step a single edge. In both networks the weights were random integer numbers between 1 and 50, drawn from a uniform distribution. Finally in both the networks the transmission rates {β i } were chosen in such a way that the R min 0 2 ½0:95; 1:05 while β max ¼ 2μ (where μ was set equal to 2). In order to do that we proceeded as follow: first we sampled {β i } from a normal distribution centered in μ and with a standard deviation of 0.2. Then we set the maximum value equal to 2μ. Finally we added a cycle in which we computed the R min 0 of the current distribution and if it did not belong to the mentioned interval every β i except the maximum one were rescaled by the current value of R min 0 . The process was repeated until R min 0 was in the desired interval.
The same sampling technique for the {β i } distribution was used in Fig. 3, but in this case we had μ = 2 and β max ¼ 6μ while R min 0 was a variable. The reason why we chose to sample the transmission rates in this way is that every curve R 0 (ϵ) is uniquely identified by its maximum and minimum value, so if one's aim is to observe how the critical point depends on R min 0 we have to keep everything else constant, even β max . If we had allowed β max to vary too we would have ended up comparing completely different curves and any information about the metacritical point would have been lost.
Finally for the simulations in Figs. 3 and 5 we used the R package called "deSolveˮ with a time step of minð1; ϵÞ=10.  49 , those of Ivory Coast (CIT) from 7 and those of Turkey (TUR) from 33 . The data regarding Italy (ITA) come from Cuebiq data (see the Data Descriptor 50 for more details on the methodology). Finally the data on USA (Florida and Colorado) were also provided by Cuebiq and are available upon request.

Data availability
Code availability