key: cord-0788320-pv5k6kwm
authors: Meng, Lingqi; Masuda, Naoki
title: Epidemic dynamics on metapopulation networks with node2vec mobility
date: 2021-06-15
journal: Journal of theoretical biology
DOI: 10.1016/j.jtbi.2021.110960
sha: 042dbf4e27624770dbcfe02a30fdf8fb49ca010a
doc_id: 788320
cord_uid: pv5k6kwm

Metapopulation models have been a powerful tool for both theorizing and simulating epidemic dynamics. In a metapopulation model, one considers a network composed of subpopulations and their pairwise connections, and individuals are assumed to migrate from one subpopulation to another obeying a given mobility rule. While how different mobility rules affect epidemic dynamics in metapopulation models has been studied, there have been relatively few efforts on comparison of the effects of simple (i.e., unbiased) random walks and more complex mobility rules. Here we study a susceptible-infectious-susceptible (SIS) dynamics in a metapopulation model, in which individuals obey a parametric second-order random-walk mobility rule called the node2vec. We map the second-order mobility rule of the node2vec to a first-order random walk in a network whose each node is a directed edge connecting a pair of subpopulations and then derive the epidemic threshold. For various networks, we find that the epidemic threshold is large (therefore, epidemic spreading tends to be suppressed) when the individuals infrequently backtrack or infrequently visit the common neighbors of the currently visited and the last visited subpopulations than when the individuals obey the simple random walk. The amount of change in the epidemic threshold induced by the node2vec mobility is in general not as large as, but is sometimes comparable with, the one induced by the change in the diffusion rate for individuals.

One reason for the scientific community to study contact networks to great extents is their relevance to infectious diseases [1] [2] [3] [4] [5] [6] . Even in just the last two decades, we have witnessed several global health threats due to infectious diseases, such as SARS pandemic from 2002 to 2003, influenza H1N1 that broke out in 2009, West Africa Ebola outbreak from 2013 to 2016, and the pandemic of COVID-19 that is impacting the entire globe as of 2021.

In the case of sexually transmitted infections, partnerships are relatively stable over time. In contrast, mobility of individuals causes contact networks to vary over time and affects dynamics of the aforementioned infectious diseases and many others. Metapopulation models provide a mathematically tractable description of epidemic processes under mobility of individual agents in networks [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] . They take advantage of the merits of both the fully mixed property of individuals within each subpopulation (also called patch) and the specific structure of connectivity between the subpopulations. A metapopulation model assumes that a population of individuals is distributed over subpopulations, which correspond to the geographical locations, such as local gathering places, cities, or counties. On the microscopic scale, the individuals are assumed to be fully mixed within each subpopulation. An infectious individual infects each susceptible individual in the same subpopulation with the same rate/probability. This assumption is practical in the absence of detailed data on the structure of interactions among the individuals within each subpopulation. On the macroscopic scale, the individuals traverse edges in the network to travel from one subpopulation to another according to a mobility rule.

Mobility patterns impact epidemic spreading, perhaps as considerably as the network structure. The simplest mobility rule is to assume that each individual moves from the currently visited subpopulation to one of its adjacent subpopulations with the same probability or with probability proportional to the weight of the edge connecting the two subpopulations. Most of early researches of metapopulation models assumed this mobility rule, which one often refers to as the simple random walk in the case of unweighted networks (i.e., no weight on edges). Then, one often lumps together the subpopulations having the same degree (i.e., number of edges that the subpopulation has) to facilitate analytical calculations [14, 15, 18, 19, 25] . However, mobility patterns for both human and animal individuals relevant to epidemic spreading are considered to be more complex than is described by the simple random walk, which has prompted studies of different mobility rules. For example, the propensity to move may depend on the degree of the subpopulation that the individual currently visits and the number of the individuals in the subpopulation [15] . Furthermore, empirical mobility patterns may be better approximated by non-simple random walks [18, [26] [27] [28] [29] [30] [31] or recurrent mobility patterns [18, 22, 23, 27, [32] [33] [34] . Other extensions of metapopulation network models include multilayer ones, in which individuals having different mobility patterns are assigned to different network layers [23, 35, 36] .

The effect of the diffusion rate of the individuals on the extent of epidemic spread is known for metapopulation models. Perhaps counterintuitively, a larger diffusion rate yields a larger epidemic threshold, which implies that epidemic spreading is less likely to occur when individuals diffuse at a faster rate [17, 22, 23] . An intuitive explanation of this result is that infectious individuals arising in subpopulations containing relatively many individuals more easily spill over to other subpopulations with fewer individuals if the diffusion rate is higher. Effects of higher-order mobility rules on epidemic spreading in metapopulation models have also been examined [31] . In contrast, how mobility rules compare with each other in terms of epidemic spreading has been underexplored. In the present study, we study the effect of the mobility rule modeled by a second-order random walk called the node2vec on epidemic spreading in metapopulation networks. The node2vec is a second-order Markovian random walk first proposed for improving the performance of data mining tasks such as classification of nodes and link prediction [37] . With the node2vec, one can tune the propensity that the individuals backtrack, which roughly corresponds to the frequency of recurrent mobility, and the weight of the local versus global search of the network. These two properties each correspond to the two parameters of the node2vec random walk. We note that the node2vec includes the simple random walk and non-backtracking random walk as special cases. There have been some theoretical results on the node2vec random walks. Qiu et al. showed that one can factorize a matrix related to the stationary distribution and transition probability of the node2vec as the length of random walks tends to infinity [38] . Furthermore, we previously showed that a node2vec random walker diffuses faster than a simple random walker in general when it infrequently travels backwards or visits the common nodes of the currently visited and the last visited nodes [39] . In the present study, we derive the epidemic threshold for the susceptible-infectious-susceptible (SIS) model when individuals obey the node2vec mobility rule. Then, we analyze the epidemic threshold for synthetic and empirical networks of subpopulations. We find that for various networks with a heterogeneous degree distribution, the epidemic threshold is large when individuals tend to avoid backtracking and visiting subpopulations that are an immediate neighbor of the presently visited and last visited subpopulations.

We examine the continuous-time SIS model on the metapopulation network when individuals move from one subpopulation to another obeying the node2vec random walks. We consider a connected, undirected, and unweighted network with N subpopulations, M edges, and N ρ individuals. Therefore, there are ρ individuals per subpopulation on average. The generalization of the following formulation to directed and weighted networks is straightforward. A subpopulation is a container of individuals and represents a habitat, a college dormitory, an Node v represents the currently visited subpopulation. Node u represents the subpopulation that the individual visited just before moving to v by traversing directed edge ei (i.e., by traversing undirected edge (u, v) in the direction from u to v). The transition probability from v to one of the four neighbors of v in the next move is given by a/(a + b + 2), b/(a + b + 2), or 1/(a + b + 2) in this example. (b) Composition of individuals in a subpopulation. urban area, or a province, for example. Each individual is in either the susceptible or infectious state at any given time. Each infectious individual independently infects each susceptible individual in the same subpopulation at rate β, which we refer to as the infection rate. This implementation of the infection process for the metapopulation model is known as the type-I reaction [5, 14] . Each infectious individual independently recovers to transit to the susceptible state at rate µ, which we refer to as the recovery rate.

The individuals also independently perform a continuous-time node2vec random walk on the network of subpopulations as follows. For infinitesimal time ∆t, an individual leaves the currently visited subpopulation, denoted by v, with probability D S ∆t or D I ∆t depending on whether the individual is susceptible or infectious, respectively. When an individual leaves v, it moves to one of the neighboring subpopulations of v. The individual moves back to the subpopulation where the individual was located just before arriving in v, which we denote by u, with the probability proportional to a. The individual moves to a subpopulation that is adjacent to both v and u with the probability proportional to b. The individual moves to a different subpopulation with the probability proportional to 1. For example, in Fig. 1(a) , the individual in its next move backtracks to u with probability a/(a + b + 2), moves to w 1 , which is adjacent to v and u, with probability b/(a + b + 2), or moves to w 2 or w 3 with probability 1/(a + b + 2) each. Note that the node2vec mobility reduces to the simple random walk if a = b = 1.

We apply the method used in Ref. [14] to simulate the stochastic SIS dynamics. We set ρ = 50. In each simulation, one individual that is selected uniformly at random is initially infectious, and the other N ρ − 1 individuals are initially susceptible. We use the rejection method with time steps of length ∆t = 10 −4 . Specifically, we repeat the following three steps. First, each infectious individual recovers to become susceptible with probability µ∆t. Second, each susceptible individual in subpopulation i becomes infectious with probability 1 − (1 − β∆t) I i , where I i is the number of infectious individuals in subpopulation i. Third, each susceptible and infectious individual independently leaves its current subpopulation with probability D S ∆t and D I ∆t, respectively. If this event happens, then the individual migrates to a neighboring subpopulation according to the node2vec mobility rule.

In each simulation, for each pair (a, b) and infection rate β, we run the simulation for 3 × 10 6 steps such that the total simulated time is 300, and then calculate the fraction of infectious individuals in the equilibrium. We calculate the average and standard deviation of the fraction of infectious individuals on the basis of 100 simulations for each parameter set.

In this section, we explain the synthetic networks used in our analysis. For each synthetic network model, except the extended ring network, we set the number of nodes, N , to 100 and repeat generating the network until we obtain a connected network.

In the Erdős-Rényi (ER) model, we place M edges selected uniformly at random between N (N − 1)/2 pairs of nodes. We set M = 300 such that the mean degree, k , is 6.

In the Barabási-Albert (BA) model, we sequentially add new nodes each with m = 3 edges that connect to existing nodes according to the linear preferential attachment rule [40] . We started the growth process from N 0 = 3 isolated nodes. We connected the first added node to each of the N 0 isolated nodes. The degree distribution p(k) approximately obeys p(k) ∝ k −3 , where ∝ represents "proportional to", in the limit of N → ∞. With N 0 = 3 and m = 3, there are 291 edges, which implies k = 5.82 ≈ 6.

The random clustered graph is specified as follows [41, 42] . For each node i, parametert i represents the number of triangles in which node i participates, ands i represents the number of edges other than those belonging to the triangles. One draws N random vectors (t i ,s i ), where i = 1, . . . , N , from a distribution F(t,s). To form a network, N i=1t i must be a multiple of 3 such that there are N i=1t i /3 triangles in the network, and N i=1s i must be an even number such that no half-edge connected to a node is left unconnected to another node. The generated network may have self-loops or multiple edges. We remove the self-loops and duplicated multiple edges from the generated network. We let F(t,s) be a doubly Poisson degree distribution [42] , i.e.,

where the parameter µ and ν are the average number of the triangles per node and the average number of edges which do not form a triangle per node, respectively. We set µ = 2.5 and ν = 1. With self-loops and multiple edges, the expected degree k = 2µ + ν = 6. After removing self-loops and duplicated multiple edges from a generated network, there are 309 edges in the final network, which implies k = 6.18. The power-law cluster graph generates a network with degree distribution p(k) ∝ k −3 with tunable clustering, i.e., density of triangles [43] . This model shares the initialization and preferential attachment rule with the BA model. However, when a new node v with m edges joins the existing network, we only carry out the preferential attachment step for the first new edge to connect v to an existing node u. For each of the other m − 1 edges from v, we carry out a triad formation step with probability p or the preferential attachment step with probability 1 − p. In a triad formation step, we add an edge from v to a uniformly randomly chosen neighbor of u. If u is one of the N 0 nodes that initially exist, all the neighbors of u may be already adjacent to v when one attempts to add an edge between v and a neighbor of u using the triad formation step. Then, no more triangles involving both v and u can be formed. In this case, we carry out the preferential attachment even if we have selected the trial formation step with probability p. We set m = 3, N 0 = 3, and p = 0.5.

The geographical threshold graph model places N nodes independently and uniformly at random in a bounded subset of the d-dimensional Euclidean space [44] . Then, we assign to each node v a weight w v . One joins each pair of nodes, v and v , by an edge if and only if

We do not assume the periodic boundary condition. We let d = 2 and the bounded subset be a unit square. We set h(r) = r −2 , θ = 80, and drew the node weights from the exponential distribution with rate parameter λ = 1 independently for the different nodes. The final network has 303 edges such that k = 6.06.

The Lancichinetti-Fortunato-Radicchi (LFR) model generates networks with community structure [45] . The degree is designed to obey a power-law distribution with power-law exponent γ, and the size of the community obeys a power-law distribution with power-law exponent κ. The model also requires the maximal degree k max and mean degree k as input. The mixing parameter µ ∈ (0, 1) specifies the fraction of edges that connect different communities. A small value of µ leads to strong community structure. We set γ = 3, κ = 2, k = 6, k max = 100, and µ = 0.1.

The extended ring network, an example of which is shown in Fig. 2 , is a network with a homogeneous degree distribution. Each node is connected to two neighbors on each side of the ring. 

We use the following two empirical networks. The US airport network has N = 332 airports and 2126 edges [46] , and it has a heterogeneous degree distribution. We ignore the weight of edges in the data set. Each subpopulation represents an airport. Two airports are connected if there exists a direct commercial flight between them.

In the Manizales network, each subpopulation represents a location that individuals visit during working days [47] . Two subpopulations are connected if at least one individual travels between the two locations. The authors of Ref. [47] categorized the individuals in the data set by their socio-economic statuses, from the lowest-income to the wealthiest, into six classes, such that the data set forms a six-layer network. We use the sixth (i.e., wealthiest) layer, which has 45 subpopulations and 194 edges. We ignore the weight of edges.

To calculate the epidemic threshold of the SIS model on metapopulation networks under the node2vec mobility, we start by considering a pair of directed edges in the opposite directions corresponding to each undirected edge. Let E = {e 1 , . . . , e 2M } be the set of directed edges. We denote each directed edge by e i = (e i (0), e i (1)) ∈ E, where e i (0) and e i (1) are the source and target node of the directed edge, respectively. One can formulate the node2vec mobility rule, which is a second-order mobility rule in the network of subpopulations, as a first-order Markov chain on a network of the 2M directed edges [29, 30, 39] . The transition-probability matrix T of the first-order Markov chain in the network of the directed edges is given by

The normalization is given by

The total infection rate in subpopulation e i (1) at time t is equal to β

where ρ S,i (t) is the number of susceptible individuals that are located in subpopulation e i (1) at time t and were located in subpopulation e i (0) just before arriving at e i (1). A similar definition applies to ρ I,i (t) (see Fig. 1(b) ). We omit t in the following text. The master equations that describe the dynamics of the number of susceptible and infectious individuals in the different subpopulations are given by

and

where i = 1, . . . , 2M . The first term on the right-hand side of Eq. (4) represents the infection events within subpopulation e i (1). The sum k;e k (1)=e i (1) ρ I,k is equal to the number of infectious individuals at subpopulation e i (1) regardless of the subpopulation where they were located just before moving to e i (1). The second term represents the recovery. The third term is the rate of diffusion out of e i (1) for the susceptible individuals. The last term is the rate of diffusion into e i (1). Note that ρ S,k T ki is the number of susceptible individuals that are located in subpopulation e k (1) and move to subpopulation e i (1). Note that T ki = 0 if e k (1) = e i (0), which excludes the possibility that individuals move to a subpopulation that is not adjacent to the currently visited subpopulation. A similar interpretation applies to Eq. (5).

To assess the accuracy of Eqs. (4) and (5) to describe stochastic SIS dynamics, we run stochastic numerical simulations of the SIS model with the node2vec mobility on a network generated by the Barabási-Albert (BA) model [40] with 100 nodes. We assume that one individual that is selected uniformly at random is initially infectious and that all the other individuals are initially susceptible. The simulation results are shown by the circles and error bars in Fig. 3 . The solid curves in Fig. 3 indicate the fraction of infectious individuals in the equilibrium that we have obtained by integrating Eqs. (4) and (5) . We confirm that the results obtained from stochastic numerical simulations are sufficiently close to those obtained from Eqs. (4) and (5) . Therefore, we will exclusively use Eqs. (4) and (5) in the following analyses.

To derive the epidemic threshold, we rewrite Eqs. (4) and (5) as follows: [48, 49] applied to Eqs. (4) and (5) with step size 0.01. We terminate the Euler method when t exceeds 300 and the difference between the current time step and the last time step in terms of the fraction of infectious individuals is smaller than 10 −9 for the first time.

where ρ S = (ρ S,1 , . . . , ρ S,2M ) , ρ I = (ρ I,1 , . . . , ρ I,2M ) , represents the transposition, • represents the Hadamard product, and (M ij ) is a 2M × 2M matrix given by

Let p * i (with i = 1, . . . , 2M ) be the stationary probability for directed edge e i , i.e.,

In other words, p * i is the stationary probability that an individual visits subpopulation e i (1) and the same individual visited subpopulation e i (0) just before e i (1). The disease-free equilibrium is given by

We apply a standard linearization technique to analyze the stability of the nonlinear dynamics represented by Eqs. (6) and (7) at the equilibrium (ρ * S ; ρ * I ). The stability is determined by the eigenvalues of the Jacobian matrix [49, 50] . We examine the Jacobian matrix of the nonlinear dynamics given by Eqs. (6) and (7) at the disease-free equilibrium. It is given by

where L is the random-walk normalized Laplacian matrix [51] , i.e.,

I is the 2M × 2M identity matrix,

and

In Eqs. (13) and (14), diag() represents the diagonal matrix whose diagonal elements are given by the arguments. Matrix J is isospectral toJ

Because L is a Laplacian matrix, all the eigenvalues of L have non-negative real parts. Because we have assumed that the original network of subpopulations is connected, the directed network induced by matrix T is strongly connected, which implies that matrix T is irreducible. By the Perron-Frobenius theorem, matrix L has 0 as a simple eigenvalue, i.e., the multiplicity of eigenvalue 0 is equal to 1. Because this zero eigenvalue corresponds to the probability flow of the node2vec random walk in the stationarity, it does not affect the stability of the diseasefree equilibrium or the determination of the epidemic threshold. Therefore, the epidemic threshold β c is given by

where spec(J 22 ) denotes the spectrum (i.e., the set of all the eigenvalues) of matrix J 22 , (spec(J 22 )) is the set of real parts of the complex numbers belonging to spec(J 22 ), and max denotes the maximum value.

In this section, we numerically examine the epidemic threshold given by Eq. (16) as a function of the weight of backtracking, a, the weight of visiting the common neighbor of the currently visited and last visited subpopulations, b, and the diffusion rate for infectious individuals, D I , for the networks introduced in sections 2.3 and 2.4. Note that the epidemic threshold is independent of the diffusion rate for susceptible individuals, D S . Because simultaneously multiplying β, µ, D S , and D I by a positive constant c is equivalent to replacing t by ct and not changing β, µ, D S , or D I , we set µ = 1 without loss of generality. Equation (14) indicates that there are three parameters a, b and D I that determine the epidemic threshold. Infection rate β and average number of individuals per subpopulation ρ only occur as their product. Therefore, we set ρ = 1 without loss of generality.

We use the bisection method to find the epidemic threshold given by Eq. (16) . More precisely, given parameters a, b, and D I , we set β lower = 0.01 and β upper = 1.5, which guarantee that max( (spec(J 22 ))) < 0 and max( (spec(J 22 ))) > 0, respectively, for all the following simulations. The eigenvalues are continuous functions of the entries of the matrix, which guarantees that spec(J 22 ) is a continuous map in terms of β [52, 53] . Because the composition of continuous maps is still continuous, the function max( (spec(J 22 ))) is continuous in terms of β. Therefore, the bisection method converges to a root. With tolerance = 10 −4 , we ran the bisection method with the initial bracketing interval [β lower , β upper ] to obtain the epidemic threshold β c with a truncation error less than .

In Fig. 4(a) , we show the epidemic threshold for the Erdős-Rényi (ER) random graph with average degree k = 6, D I = 1 and various values of a and b. We observe that the epidemic threshold increases as a or b decreases, which suggests that epidemic spreading is suppressed when the individuals travel without frequent backtracking or frequent visiting to common neighbors of the presently visited subpopulation and the last visited subpopulation. We observe that the dynamic range of the epidemic threshold β c when we fix a and vary b ∈ [0, 5] is larger than when we fix b and vary a ∈ [0, 5]. Therefore, the weight of visiting the common neighbors more strongly influences the epidemic threshold than the weight of backtracking. The epidemic threshold when we set b = 1 and vary a and D I and when we set and a = 1 and vary b and D I is shown in Fig. 4(b) and Fig. 4(c) , respectively. These two figures indicate that a larger diffusion rate D I of the infectious individuals suppresses epidemic spreading, which is consistent with the previous results for the first-order mobility rule (i.e., simple random walk) on the network of subpopulations [17, 22, 23] . Figures 4(b) and 4(c) also indicate that changes in D I ∈ [0, 10] have a larger impact on the epidemic threshold than changes in a ∈ [0, 5] and b ∈ [0, 5], whereas these parameters have different units.

We show the epidemic threshold for the BA model [40] in Figs. 4(d)-(f), random clustered graph [41, 42] in Figs. 4(g)-(i), power-law cluster graph [43] in Figs. 4(j)-(l), geographical threshold graph model [44] in Figs. 4(m)-(o), and Lancichinetti-Fortunato-Radicchi (LFR) model [45] in Figs. 4(p)-(r), for the same parameter regions as those used for the ER random graph. Figures 4(d) , (g), (j), (m), and (p), in which we vary a and b with D I = 1, suggest that the epidemic threshold largely decreases as a or b increases for all model networks, which is similar to the results for the ER random graph shown in Fig. 4(a) . With an exception of the geographical threshold graph model at small a and b values (see Fig. 4(m) ), the epidemic threshold is largest near a = b = 0. We also find for these networks that the impact of b on the epidemic threshold is stronger than that of a, which is consistent with Fig. 4(a). Figures 4(b) , (c), (e), (f), (h), (i), (k), (l), (n), (o), (q), and (r) indicate that a large diffusion rate suppresses epidemic spreading. Figures 4(b) , (e), (h), (k), (n), and (q) indicate that the effect of a on the epidemic threshold is much smaller than that of D I . These results are consistent with those for the ER random graph. In contrast, the effect of b is smaller than D I in some cases (Figs. 4(i) and (l)), which is similar to the case of the ER random graph (Fig. 4(c) ), but is comparable to that of D I for some networks (Fig. 4(f), (o), (r) ).

The observation above may fail for networks with a homogeneous degree distribution. To examine this possibility, we consider an extended ring network shown in Fig. 2 , in which each node is connected to two neighbors on each side of the ring. For this network, we have derived the following theorem.

Theorem 1: The epidemic threshold for the extended ring network is equal to µ/ρ, independently of the value of a, b, and D I .

Note that the statement of the Theorem does not hold true in general if the network is regular (i.e., a network in which all the nodes have the same degree) but is different from the extended ring network. The reason for this difference is probably that the extended ring network is not only regular but also vertex-transitive [39, 54, 55] ; casually speaking, all the nodes are equivalent to each other. We provide the proof of Theorem 1 in Appendix. This result contrasts to the numerical results shown in Figs. 4 and 5, in which the epidemic threshold systematically depends on the a, b, and D I values. The reason for this difference is probably that the networks used in Figs. 4 and 5 are relatively heterogeneous in terms of the node's degree, whereas all the nodes have the same degree in the extended ring network.

We find similar results for two empirical networks. The dependence of the epidemic threshold on a, b, and D I for these two networks is shown in Fig. 5 . The results are similar to those for the model networks shown in Fig. 4 .

We investigated the SIS dynamics on metapopulation networks with the node2vec mobility. The node2vec is a second-order random walk, which one can transform into a first-order Markov chain by considering the set of all the directed edges and transitions among them [29, 30, 39] . We built the master equations for the SIS dynamics involving the transition probability matrix of the node2vec random walk, T . We derived the epidemic threshold for any connected undirected networks. We observed that, for various synthetic and empirical networks with heterogeneous degree distributions, the epidemic threshold was large (and therefore, epidemic spreading is suppressed) when individuals explore the network without frequently backtracking (i.e., small a) or without frequently visiting common neighbors of the currently visited and the last visited subpopulations (i.e., small b).

The convergence of the node2vec random walk towards the stationary density is fast when a or b is small [39] . In the present study, we have shown that the epidemic threshold tends to be large for small a or b values. Therefore, a fast convergence in the node2vec random walk is associated with a large epidemic threshold. When individuals obey the simple random walk in metapopulation models, the epidemic threshold for the SIS dynamics is known to increase as the diffusion rate increases [17, 22, 23] , which is consistent with our present numerical results. Because the rate of the convergence is generally proportional to the diffusion rate, a fast convergence in the simple random walk is associated with a large epidemic threshold. Therefore, our present result that the epidemic threshold tends to be large for small a, small b, or large D I values collectively supports the idea that the epidemic threshold is large when the random walk rapidly converges towards the stationary density.

When individuals move according to Markovian random walks with recurrent mobility, the epidemic threshold does not necessarily increase monotonically as a function of the diffusion rate, and the relationship between the epidemic threshold and the diffusion rate depends on which subpopulations the individuals use as their home subpopulation [22, 23, 33] . We did not find such a nonmonotonicity because our models do not describe recurrent mobility with which different individuals use different subpopulations as home. It may be interesting to combine the recurrent mobility modeling and the node2vec random walk. Then, we may observe nonmonotonic dependence of the epidemic threshold on the a or b values as well as on the diffusion rate. This warrants future work.

Let us mention other possible directions for future research. First, we may be able to exploit our observation that the node2vec mobility with small a and b suppresses the spread of infections to inform intervention methods. A family of methods to intervene into epidemic dynamics on networks of subpopulations is to lessen the infection rate within subpopulations having large degrees [56] [57] [58] . Under the simple random walk, containment of the epidemic has also been examined in combination with adaptive dynamics, with which individuals cancel their travel in response to their infection status [28, 59] or adjust their movement based on the information about the safety level measured by the number of susceptible individuals in various subpopulations [59] [60] [61] [62] . Other containment methods include the shutdown of some edges connecting subpopulations [12, 63] , usage of antibiotics and antiviral drugs [63, 64] , and finding the most influential spreaders [65, 66] . It is worth considering these containment methods assuming the node2vec mobility because the node2vec random walk may change the efficiency of these intervention methods. Furthermore, forcing individuals to use small a or b values itself may also be used for intervention, which may be combined with the aforementioned intervention strategies. Second, in some multilayer networks, the epidemic threshold is not a smooth function of the diffusion rate when the interlayer coupling is weak [23] . Examining the possibility of similar non-smooth behavior of the epidemic threshold in terms of the a and b under the node2vec mobility may be interesting. Finally, it is straightforward to extend our modeling framework to analyze the effects of the node2vec mobility on other dynamical processes on metapopulation model networks, such as the susceptible-infectious-recovered (SIR) model [25, 57] , evolutionary games [67, 68] , and prey-predator dynamics [69, 70] . matrices B i , where i = 1, 2, . . . , n, we define the kn × kn block circulant matrix bcirc(B 1 , B 2 , . . . , B n ) by bcirc (B 1 , B 2 

We order the directed edges in E in the extended ring network as illustrated in Fig. 6 . Then, the transition probability matrix, T , is block circulant and given by 

where

and

In this case, matrix

is also block circulant, where

and

We also obtain

where

and N is the normalization constant. By combining Eqs. (14) , (18) , (23) , and (29), we obtain Next, let

denote the N th roots of 1, where i is the imaginary unit and j = 0, 1, 2, ..., N − 1. Using ρ j , we define 4 × 4 matrices

Theorem 3 in Ref. [71] guarantees that

Therefore, the epidemic threshold β c is given by β c = max{β | max( (spec(J 22 ))) = 0} = max{β | max( (spec(H j ))) = 0, j = 0, . . . , N − 1}.

The max in the second line of Eq. (34) represents the maximum among the four eigenvalues of H j and among the N values of j. The construction up to here is similar to that in our previous work for computing the spectral gap of node2vec random walks on networks [39] .

To derive β c , we further proceed as follows. Using Eq. (31) and (32), we obtain

where j = 0, . . . , N − 1, and s = 2 + 2(a + 2b + 1)

In particular, we obtain 

The main objective of the remainder of the proof is to show that the eigenvalue of J 22 with the largest real part is completely determined by the matrix H 0 . The sum of each column of H 0 is equal to βρ − µ. Therefore, we distinguish three cases and claim the following statements:

(S1) If β < µ/ρ, then each Gershgorin circle of H 0 created by each of its columns is contained in the left half of the complex plane excluding the imaginary axis. Therefore, the real part of each eigenvalue of H 0 is negative, i.e., (spec(H 0 )) < 0.

(S2) If β = µ/ρ, then the row vector 1 = (1, 1, 1, 1) is a left eigenvector of H 0 associated with eigenvalue 0, i.e., 1H 0 = 0.

(S3) If β > µ/ρ, then βρ − µ is a positive eigenvalue of H 0 with the associated eigenvector 1.

To verify statement (S2), we evaluate each coordinate of row vector 1H 0 :

where i = 1, 2, 3, 4. This finishes the proof of statement (S2). One can verify statement (S3) in the same manner.

To prove (S1), we assume that β < µ/ρ. We consider the Gershgorin circle created by the first column of H 0 as an example. The coordinate of the center of this Gershgorin circle is (C 1 , 0), where

In the first inequality in Eq. (38), we used β < µ/ρ. In the second inequality, we used a > 0, b > 0, and s > 2; note that Eq. (36) implies s > 2. The radius of the corresponding Greshgorin circle, denoted by R 1 , is given by

Note that C 1 + R 1 is the sum of the first column of H 0 , which is equal to βρ − µ < 0. Therefore, the Gershgorin circle is contained in the left half of the complex plane without touching the imaginary axis. Similarity, it is straightforward to show that all the Gershgorin circles of H 0 created by its columns are contained in the left half of the complex plane without touching the imaginary axis. This completes the proof of (S1). Now we assess the Gershgorin circles generated by H j with j = 1, . . . , N − 1. When β < µ/ρ, the diagonal elements of H j satisfy

for any i = 1, . . . , 4 and j = 1, . . . , N − 1, because multiplication by ρ j is a rotation about the origin in the complex plane. For example, one obtains

Equation (40) implies that the center of the Gershgorin circle for H j created by its ith column is located to the left of the center of the corresponding Gershgorin circle for H 0 in the complex plane. It also holds true that 4 =1; =i

for any i = 1, . . . , 4 and j = 1, . . . , N − 1. 

Equation (42) implies that the radius of the Gershgorin circle for H j created by its ith column is smaller than the radius of the corresponding Gershgorin circle for H 0 . Because each Gershgorin circle for H j , where j = 1, . . . , N − 1, has its center to the left of that of the corresponding circle for H 0 and its radius is smaller than that of the corresponding circle for H 0 , all the Gershgorin circles for H j (with j = 0, 1, . . . , N − 1) are contained in the left half of the complex plane without touching the imaginary axis. Recall that spec(J 22 ) = N −1 j=0 spec(H j ). Therefore, all the eigenvalues of J 22 have negative real part when β < µ/ρ. This situation is schematically shown in Fig. 7 . 41) and (43) when β < µ/ρ. The larger Gershgorin circle with radius R = 4 =1; =i |(H0) 1 | is centered at ((H0)11, 0). The smaller Gershgorin circle with radius r = 4 =1; =i |(Hj) 1 | is centered at ( ((Hj)11), ( D I a+b+2 ρ 2 j )), where and represent the real and imaginary part of the complex number, respectively.

The case β = µ/ρ is similar to the case β < µ/ρ. The only difference is that the four Gershgorin circles generated by H 0 are tangent to the imaginary axis. Therefore, we obtain max( (spec(H j ))) < 0 for j = 1, . . . , N − 1, and max( (spec(H 0 ))) ≤ 0. Furthermore, by statement (S2), H 0 has an eigenvalue 0. Therefore, when β = µ/ρ, all the eigenvalues of J 22 have negative real part except the zero eigenvalue. Statement (S3) implies that J 22 has a positive eigenvalue βρ − µ when β > µ/ρ. These results verify that β c = µ/ρ is the epidemic threshold.

The empirical network data sets are open resources and available at [46, 47] . The Python codes and synthetic network data sets used in the present study are available on Github [72] .

Networks and epidemic models

Networks and the epidemiology of infectious disease

Predicting and controlling infectious disease epidemics using temporal networks

Eight challenges for network epidemic models

Epidemic processes in complex networks

An immunization model for a heterogeneous population

Spatial heterogeneity and the design of immunization programs

Spatial heterogeneity in epidemic models

(Meta) population dynamics of infectious diseases

Cities and villages: infection hierarchies in a measles metapopulation

Forecast and control of epidemics in a globalized world

The role of the airline transportation network in the prediction and predictability of global epidemics

Reaction-diffusion processes and metapopulation models in heterogeneous networks

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

Multiscale mobility networks and the spatial spreading of infectious diseases

Effects of diffusion rates on epidemic spreads in metapopulation networks

Phase transitions in contagion processes mediated by recurrent mobility patterns

Modelling dynamical processes in complex socio-technical systems

A metric of influential spreading during contagion dynamics through the air transportation network

The scaling of human contacts and epidemic processes in metapopulation networks

Critical regimes driven by recurrent mobility patterns of reaction-diffusion processes in networks

Spreading processes in multiplex metapopulations containing different mobility networks

Vectorborne epidemics driven by human mobility

Invasion threshold in heterogeneous metapopulation networks

Natural human mobility patterns and spatial spread of infectious diseases

Invasion threshold in structured populations with recurrent mobility patterns

Human mobility and time spent at destination: impact on spatial epidemic spreading

Memory in network flows and its effects on spreading dynamics and community detection

Causality-driven slow-down and speed-up of diffusion in non-markovian temporal networks

Assessing reliable human mobility patterns from higher order memory in mobile communications

Epidemic spreading in localized environments with recurrent mobility patterns

Impact of temporal scales and recurrent mobility patterns on the unfolding of epidemics

Epidemic spreading in heterogeneous networks with recurrent mobility patterns

Reaction-diffusion processes and metapopulation models on duplex networks

Epidemic spread on interconnected metapopulation networks

Proc. 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining

Network embedding as matrix factorization: Unifying deepwalk, line, pte, and node2vec

Analysis of node2vec random walks on networks

Emergence of scaling in random networks

Percolation and epidemics in random clustered networks

Random graphs with clustering

Growing scale-free networks with tunable clustering

Geographical threshold graphs with small-world and scale-free properties

Benchmark graphs for testing community detection algorithms

Pajek datasets

Rich do not rise early: Spatio-temporal patterns in the mobility networks of different socio-economic classes

An Introduction to Numerical Analysis

Nonlinear Dynamics and Chaos with Student Solutions Manual: With Applications to

Differential and Integral Inequalities

Random walks and diffusion on networks

Solutions of Equations in Euclidean and Banach Spaces

An alternate proof of the continuity of the roots of a polynomial

Algebraic Graph Theory

Algebraic Graph Theory

Random and targeted interventions for epidemic control in metapopulation models

Intervention threshold for epidemic control in susceptible-infected-recovered metapopulation models

Modelling the effect of heterogeneous vaccination on metapopulation epidemic dynamics

Modeling human mobility responses to the large-scale spreading of infectious diseases

Safety-information-driven human mobility patterns with metapopulation epidemic dynamics

The price of anarchy in mobility-driven contagion dynamics

Interplay between epidemic spread and information propagation on metapopulation networks

Strategies for mitigating an influenza pandemic

Modeling the worldwide spread of pandemic influenza: baseline case and containment interventions

Identification of influential spreaders in complex networks

Identification of influential spreaders in complex networks using hybridrank algorithm

Metapopulation model for rock-paper-scissors game: mutation affects paradoxical impacts

Evolutionary vaccination game approach in metapopulation migration model with information spreading on different graphs

Demographic and environmental stochasticity in predator-prey metapopulation dynamics

Intermediate fragmentation per se provides stable predator-prey metapopulation dynamics

Eigenvectors of block circulant and alternating circulant matrices

Data sets and Python codes for epidemic dynamics on metapopulation networks under node2vec mobility

The outline of the following proof is to show that, as β increases through β = µ/ρ, there is an eigenvalue of J 22 whose real part changes from negative, to zero, and then to positive.We first present an explicit expression for matrix J 22 when the network is an extended ring. Given k × k