Metapopulation models imply non-Poissonian statistics of interevent times

Interevent times in temporal contact data from humans and animals typically obey heavy-tailed distributions, and this property impacts contagion and other dynamical processes on networks. We theoretically show that distributions of interevent times heavier-tailed than exponential distributions are a consequence of the most basic metapopulation model used in epidemiology and ecology, in which individuals move from a patch to another according to the simple random walk. Our results hold true irrespectively of the network structure and also for more realistic mobility rules such as high-order random walks and the recurrent mobility patterns used for modeling human dynamics.

However, these models do not explain the genesis of heavy-tailed IET distributions in the presence of mobility of individuals. Human and animal individuals move around to meet others, transmitting information and disease, and then depart from each other. In fact, most of face-to-face contact or proximity data have been collected from mobile individuals, and such data show heavy-tailed IET distributions [25,[39][40][41][42][43][44][45][46][47]. Previous modeling studies showed that encountering events of random walkers moving on a two-dimensional plane with heterogeneous attractiveness values generate heavy-tailed IET distributions [48][49][50][51][52]. These studies crucially assumed that each walker is endowed with a randomly assigned at- * naokimas@buffalo.edu tractiveness value and that walkers tend to slow down when they approach attractive walkers nearby. In the present study, we are interested in analytically accounting for heavy-tailed IET distributions under mobility in a simpler manner, i.e., without introducing the concept of attractiveness or any heterogeneity among walkers.
In the present study, we provide accounts of heavytailed IET distributions when individuals move around, meet, and separate according to a standard metapopulation network model, in which the individuals perform simple random walks, and its more realistic variants. We analytically show that the IET distribution for a pair of individuals is a mixture of exponential distributions, which robustly holds true for different network structures and mobility rules. Crucially, mixtures of exponential distributions produce IET distributions that are substantially closer to heavy-tailed distributions observed in empirical data than exponential distributions do.

II. MODEL
We consider a metapopulation network with N subpopulations (i.e., nodes). The network may be directed and/or weighted. The individuals that populate the network are assumed to be simple random walkers in continuous time. Different walkers move independently of each other and may produce contact events (events for short) between them only when they are copresent in the same subpopulation.
Specifically, the time t for which each walker waits until it leaves the current subpopulation is independently arXiv:2106.10348v2 [physics.soc-ph] 7 Feb 2022 distributed according to an exponential distribution with rate D. The mobility and contact dynamics are illustrated in Fig. 1. When the walker leaves the ith subpopulation, it moves to a neighboring jth subpopulation with probability A ij / N =1 A i , where A ij is the weight of edge (i, j). Events between a pair of walkers in the same subpopulation occur according to a Poisson process with rate λ, which epidemic process models on metapopulation networks usually assume. In other words, the distribution of time t to the next event between a pair of walkers in the same subpopulation is given by φ e (t) = λe −λt . When the walkers are in different subpopulations, no event occurs between them.

A. General solution
Consider two walkers visiting the same subpopulation. We derive the probability density of an IET, τ , between an arbitrary pair of individuals, which are random walkers. When the two walkers are copresent in a subpopulation, events (including the case of 0 event) occur between them before either walker leaves the subpopulation. After that, a series of n copresences occurs before the next event occurs (see Fig.1).
We let p(τ, n) be the probability of generating an IET of length τ during the nth copresence after either walker leaves the current subpopulation. The probability density of the time at which either walker leaves the subpopulation where the two walkers are copresent is given by φ d (t) = 2De −2Dt . We obtain By Laplace transforming p(τ, 0), we obtain .
We denote the distribution of time t between consecutive copresences of the two walkers, which we refer to as the inter-copresence time (ICT), by φ c (t). Then, the joint probability density of τ and n = 1 is given by where S e (t) = e −λt and S d (t) = e −2Dt are the survival functions of φ e (t) and φ d (t), respectively, and S c (t) = ∞ t dt φ c (t ) and h c (t) = φ c (t)/S c (t) are the survival and hazard functions of φ c (t), respectively. By Laplace transforming p(τ, 1), we obtain Similarly, the Laplace transform of the joint probability density of τ and n is given bŷ Therefore, the IET distribution in the frequency domain is given bŷ Using Eq. (6), we obtain the following expression of the coefficient of variation (CV), defined as the standard deviation divided by the mean, of τ (see Appendix A for the derivation): whereφ c (0) andφ c (0) are the first and second derivatives ofφ c (s) evaluated at s = 0, respectively.

B. Exponential ansatz
The asymptotic form of the distribution of first-passage time of random walks in networks with short relaxation time has an exponential tail [66]. With φ c (t) = αe −αt , which we refer to as the exponential ansatz, Eq. (6) is reduced top The inverse Laplace transform of Eq. (8) is given by where ∆ ≡ (λ + 2D + α) 2 − 4αλ > 0. Therefore, the distribution of IET is a mixture of two exponential distributions (see Appendix B for a proof). By substituting φ c (t) = αe −αt in Eq. (7), we obtain FIG. 1. Schematic illustration of the dynamics of two walkers in a network with four subpopulations. Initially, at t = t0, the walkers are copresent in subpopulation 1. Then, events between them occur before t1. At t = t1, one walker moves from subpopulation 1 to subpopulation 4. At t = t2, the other walker moves from subpopulation 1 to subpopulation 2. At t = t3, the same walker moves to subpopulation 4 to be copresent with the other walker. Therefore, the first ICT is given by t3 − t1.
Although the walkers are copresent for t ∈ (t3, t4), no event occurs. At t = t4, one walker moves to subpopulation 3. At t = t5, the other walker also moves to subpopulation 3, producing a second ICT, which is equal to t5 − t4. Finally, an event occurs during their copresence in subpopulation 3, yielding an IET, denoted by τ * . In this example, this IET is followed by several much shorter IETs.

C. Intercopresence time distribution of the walkers
To mathematically derive the distribution of ICTs, φ c (t), we denote by s = (m, n) the state of the system in which one walker is in the mth subpopulation and the other walker is in the nth subpopulation. Because the walkers are indistinguishable, the system has N = N (N + 1)/2 states, i.e., (m, n), where 1 ≤ m ≤ n ≤ N .
We consider the first-passage time of a continuoustime Markov chain from a state (i, i), i.e., a copresence, to a state (j, j), i.e., the next copresence, where i, j ∈ {1, 2, . . . , N }. A similar problem of first-passage time to the next copresence of two random walkers appears in the theory of cooperation in evolutionary games on networks [67,68]. We refer to the states in which the two walkers are in the same subpopulation as absorbing states and those in which the walkers are in different subpopulations as transient states. Let A = {(m, n); m = n} and B = {(m, n); m < n} be the set of absorbing and transient states, respectively. We denote the probability that the dynamics starting in state (i, i) at time 0 is in state s = (m, n), where m ≤ n, at time t by p i, s (t), or equivalently, p i,(m,n) (t). We set At time 0, one of the walkers leaves subpopulation i such that the system leaves the absorbing state (i, i) and enters a transient state. Therefore, the distribution of the initial state of the Markov process is given by where k out A il is the out-degree of subpopulation i, and δ ij is the Kronecker delta.
The master equation for a transient state s = (m, n) is given by where L ≡ 2DI −Q, matrix Q is the matrix of transitions rates between transient states, and I is the N T × N T identity matrix (see Appendix C for details). Equation (13) leads to The master equation for the absorbing states is given by where R is the matrix of transition rates from a transient state to an absorbing state (see Appendix C). By substituting Eq. (14) into Eq. (15), we obtain The probability density function of the first-passage time from (i, i) to a state s, denoted by f i, s (t), satisfies where f A i, s (t), or equivalently f A ij (t), is the probability density of the first-passage time from (i, i) to (j, j), because p s, s (t − t ) = 1 for any t ≥ t . Therefore, the probability density of the first-passage time to the absorbing states is given by where is the probability density of the first-passage time from (i, i) to (j, j).
We define the N × N T matrix P T (t) such that its ith row is given by where P T 0 = P T (0). We are interested in the IET distribution in the equilibrium. In the equilibrium, the initial copresence state (i, i) should be the stationary probability of the discretetime random walk on A, where the transition probability T ij from (i, i) to (j, j) is given by Using Eqs. (20) and (21), we obtain the N × N matrix We obtain the stationary distribution q * for the absorbing states as which we can numerically solve for arbitrary networks of subpopulations.
If Q is diagonalizable, we can write where v j and w j are the right and left eigenvectors, respectively, associated with eigenvalue γ j of Q. By combining Eqs. (20), (23), and (24), we obtain the ICT distribution weighted by the stationary probability of the initial location of the two copresent walkers as where c j = q * P T 0 v j w j 1 N T , α j = 2D − γ j > 0, and 1 N and 1 N T are the column vectors of size N and N T , respectively, in which all the elements are 1 [see Appendix D for a proof of Eq. (25) and that α j > 0]. Because is a mixture of exponential distributions, including the case in which c j < 0.
By substituting the Laplace transform of Eq. (25) into Eq. (6) and Eq. (7), we obtain the IET distribution in the frequency domain and the CV aŝ and respectively. This solution depends on the network structure through φ c (t).
To compare CV ansatz with CV mixture we impose that the mean ICT is equal between the two. This condition combined with φ c (t) = αe −αt and Eq. (25) yields Then, CV mixture is larger than CV ansatz because 2 , which is satisfied by the Sedrakyan's inequality (also known as Titu's lemma) [69]. Therefore, the exponential ansatz, which has led to the mixture of two exponential distributions for IETs, gives a lower bound in terms of the dispersion of IETs. Additionally, in the case of the network of two subpopulations connected to each other, Eq. (25) is reduced to the exponential ansatz (see Appendix D). Therefore, although the time-domain solution of the IET distribution is not available in general, the IETs are guaranteed to be distributed according to a distribution with a heavier tail than a mixture of two exponential distributions.
In fact, if the system starts in state (i, i) and reaches (j, j) to produce an ICT, it restarts from (j, j) to produce the next ICT. Therefore, Eq. (26), which assumes the independence of different ICTs, is only approximate. Nevertheless, the following numerical simulations support that the difference between the approximate solution [i.e., Eqs. (26) and (27)] and the exact solution (see Appendix E) is negligible.

D. Numerical results
We simulated the model to validate our theory. We show the IET distribution produced by two walkers on the Barabási-Albert (BA) and Watt-Strogatz (WS) networks in Figs. 2(a) and 2(b), respectively. The figure shows that the IET distributions obtained from numerical simulations and ansatz have heavier tails than the exponential distribution whose mean IET is equal to that for the numerical simulations. The distribution obtained from the simulation decays more smoothly than the exponential ansatz.
To be more quantitative and general in terms of the parameter values and variety of networks, we compared the numerically and theoretically obtained CV of IET on Survival function of IETs, S(τ ). We show the results for the direct numerical simulation of the model (simulation), the exponential ansatz [ansatz; Eq. (9)], and the single exponential distribution (exponential) whose mean is equal to that for the direct simulations. In Eq. (10), we set  10) and (27)], which only depends on D/λ. Therefore, we set λ = 1 without loss of generality and varied D.
The results shown in Fig. 3 indicate that the CV is large when D/λ is small for all four model networks and two empirical networks. Moreover, the CV obtained from our theory [Eq. (27)  . Therefore, spatiality, or the large average path length between nodes, which is present in the last two networks but not in the other networks, may negatively affect the accuracy of the exponential ansatz.
Next, we investigated the effects of the network size on the IET. We set λ = 1 and D = 0.25, and computed the CV value for each of the six networks shown in Fig. 3 for different numbers of nodes, i.e., N = 50, 100, 200, 500, 1000, and 2000. We used the same parameter values as those used in Fig. 3 for each network model.
The CV values for the different network models and different N values are shown in Fig. 4. We find that the results are roughly independent of N , except for the geographical threshold graph. For the geographical threshold graph, the CV is large at small N . Nevertheless, the CV value for the geographical threshold graph depends little on N for larger networks.
bility [76,77]. Our theory and the result that IETs obey a heavier-tailed distribution than the exponential distribution holds true for second-and higher-order random walks (see Appendix F).

E. Work-home model
More realistic metapopulation models incorporate the recurrent nature of human mobility where individuals iterate between a home subpopulation and a work subpopulation, which are different for different individuals [78][79][80][81][82][83]. In this section, we analyze the IET distribution for such a mobility rule.
We emulate this situation by assigning to each walker a home subpopulation and allowing it to visit subpopulations adjacent to the home subpopulation. Then, each walker is confined to a star graph, of which the central subpopulation is the home subpopulation and the leaves are other locations such as work. Without loss of generality, we consider two random walkers on the respective star graphs to generate sequences of IETs.
We examine two cases. In the first case, the two walkers share the home subpopulation and all leaf subpopulations, as shown in Fig. 5(a). In other words, the two walkers are confined to the same star graph, and events between them may occur in any subpopulation. Therefore, our theory directly applies.
In the second case, we assume that the two walkers have different home subpopulations and share some but not necessarily all the leaf subpopulations, as shown in Fig. 5(b). This case mimics, for example, the situation in which the two individuals are coworkers living in different cities. Events may occur between the walkers only when they are copresent in any of the shared leaf subpopulations. To derive the states and the transition rates between the states, we distinguish the two walkers and the subpopulations. There are three types of subpopulations: the home subpopulations, which are unshared, the shared leaf subpopulations, which are work locations, for example, and the unshared leaf subpopulations, which represent other locations that one but not both walkers visits. In a network with N subpopulations, two subpopulations are home subpopulations, N s subpopulations are the shared leaf subpopulations, and walkers 1 and 2 have N 1 and N 2 unshared leaf subpopulations, respectively, such that N = 2 + N s + N 1 + N 2 .
We define the home subpopulation nodes by h 1 and h 2 , and the set of shared leaf subpopulations by W = {w 1 , w 2 , . . . , w Ns }. There are four types of transient states. In the first type of transient states, both walkers are in different shared leaf subpopulations. There are N s (N s − 1) such states. In the second type, walker 1 is in a shared leaf subpopulation and walker 2 is not There are N s (N 1 + 1) such states. In the third type, walker 2 is in a shared leaf subpopulation and walker 1 is not There are N s (N 2 + 1) such states. In the fourth type, both walkers are not in any of the shared leaf subpopulations. There are (N 1 + 1)(N 2 + 1) such states. Therefore, there are where χ W is the indicator function defined by In an absorbing state, the two walkers are copresent in one of the shared leaf subpopulations. Therefore, there are N s absorbing states. Each element of the N T × N w matrix of transition rates from a transient state to an absorbing state, R s s , where s = (m 1 , m 2 ) ∈ B and s = (m 1 , m 2 ) ∈ A, is given by The distribution of the initial state of the Markov process, is given by Then, with Eqs. (30), (32), and (33) as inputs to our theory, the calculation steps to derive the distribution and CV of IET remain the same.
The CV values obtained from the simulation and the theory for the case in which the two walkers share the home and the case in which they have different homes are shown in Figs. 5(c) and 5(d), respectively. In both cases, the results are similar to those shown in Fig. 3, i.e., the CV is large when D/λ is small, and the theory agrees with the simulation. The exponential ansatz solution yields slightly smaller values of CV than the numerical simulations, which we also observe in Fig. 3.

IV. DISCUSSION
We have shown that IET distributions are mixtures of exponential distributions for mobile individuals in metapopulation networks. The results hold true under mild conditions, i.e., for various structures of the metapopulation network, a work-home mobility rule, and higher-order random walks. Although a mixture of exponential distributions is technically not heavytailed, it often approximates heavy-tailed distributions reasonably well over scales [35,37,84,85]. Therefore, the present results provide a compelling explanation of heavy-tailed IET distributions in human and animal contact data. Additional mechanisms such as circadian or weekly rhythms [33] and dynamics of individuals' internal states (e.g., high-activity versus low-activity states) [39] on top of mobility and metapopulation networks may make IET distributions more smooth and more powerlaw-like.
We assumed that two walkers in the same subpopulation meet each other according to a Poisson process. In other words, the waiting time until they have the next event obeys an exponential distribution if neither walker leaves the current subpopulation. In fact, this assumption is unnecessary for our results to hold true. One can assume other types of distributions for the mentioned waiting time to obtain qualitatively the same results. This is because the result that the IET distribution is a mixture of exponential distributions owes to the fact that the ICT distribution is a mixture of exponential distributions.
The present work bridges statistics of IETs and metapopulation models, which have been extensively but in most cases separately investigated topics in network science and related research fields. For example, heavy-tailed IET distributions suppress epidemic spreading in a majority of scenarios [7][8][9][10][11][12][13][14]. The present results allow a new interpretation of these results. That is, different results of epidemic spreading associated with non-Poissonian IET statistics may owe in large part to whether the structure of the metapopulation network plays a significant role (which would yield heavier-tailed IET distributions) or the population is sufficiently wellmixed (which would yield an exponential IET distribution). The effects of IET statistics other than their heavy-tailed distributions [2,3,86] on epidemic and other dynamical processes may also be due to the underlying metapopulation network. Reexamining the role of IETs in contagion and other dynamical processes from the viewpoint of metapopulation networks warrants future work. The coefficient of variation (CV) is defined as where · is the expectation. From Eq. (6), we obtain the first and second moments of τ as and (A3) In this section, we show that Eq. (9) is a mixture of two exponential distributions. We first rewrite Eq. (9) as where and The two exponents in Eq. (B1) are positive because the weights C 1 and C 2 are positive. Equations (B2) and (B3) also imply that Therefore, Eq. (B1) is a mixture of two exponential distributions.

Appendix C: Transition rate matrix
An element of the transition rate matrix W from a state s = (m, n) to a state s = (m , n ) is given by The first term on the right-hand side of Eq. (C1) represents the transitions between transient states. We define Q as the N T × N T matrix of transition rates between transient states, i.e., W s s = Q s s , where s ∈ B and s ∈ B. Then, The second term on the right-hand side of Eq. (C1) represents the transitions from a transient state to an absorbing state. We define R as the N T × N matrix of transition rates from a transient state to an absorbing state, i.e., W s s = R s s , where s ∈ B and s ∈ A. Then, Note that for any s ∈ B. Equation (C4) represents the fact that the system leaves any transient state at rate 2D owing to the movement of each walker, which occurs at rate D.
The master equation for a transient state s = (m, n), which is equivalent to Eq. (13), is given by In Eq. (C5), we used the fact that the sum of transition rates from any transient state to other states is equal to 2D, as shown in Eq. (C4).

Appendix D: Distribution of inter-copresence times
To derive the ICT distribution, we first note that Eq. (C4) leads to Then, by combining Eqs. (20), (23), and (D1), we obtain By substituting Eq. (24) into Eq. (D2), we find In our model, all transient states are reachable from any transient state. Therefore, Q is irreducible and, by definition, non-negative. Hence, by the Perron-Frobenius theorem, the spectral radius of Q, denoted by ρ(Q), satisfies ρ(Q) ≤ max s s ∈B Q s s . Using Eq. (C4), we obtain ρ(Q) ≤ 2D, which implies that 2D − γ j ≥ 0, for all j.
In the particular case of two subpopulations connected to each other (N = 2), the number of transient states is N T = 1. Therefore, Eq. (D3) becomes which is equivalent to the exponential ansatz. In this case, there is one transient state and two absorbing states. Therefore, Q = 0, R = D D , L −1 = 1/2D, and 1 N T = 1. From Eq. (22), we obtain such that q * = 1/2 1/2 . Then, we obtain α 1 as Therefore, for the case of two subpopulations, the ICT distribution is given by Appendix E: Exact solution for the interevent time distribution In this section, we derive the exact solution for the IET distribution. The probability density with which the two walkers are copresent in subpopulation j after time t given that the last copresence terminated at time 0 in subpopulation i is given by f A ij (t), which is given via Eq. (20). In other words, f A ij (t) is the probability density of an ICT, t, from (i, i) to (j, j). Note that f A ij (t) = φ c (j, t|i) and the normalization is given by N j=1 ∞ 0 dt f A ij (t) = 1. We denote by p ij (τ, n) the joint distribution of an IET, τ , and the number of copresences, n, such that the event has occurred in subpopulation i and the next event occurs in subpopulation j. The normalization is given by N j=0 ∞ n=0 ∞ 0 dτ p ij (τ, n) = 1. We derive the Laplace transform of p ij (τ, n) by extending Eq. (5) as follows: We define the N × N matrixΦ(s) by [Φ(s)] ij = φ c (j, t|i). Using Eq. (20), we obtain We also define the N × N matrixΠ(s, n) by [Π(s, n)] ij =p ij (s, n). Using Eq. (E1), we obtain (E3) We denote byĝ(s) the stationary distribution in the frequency domain of IET weighted by the stationary probability of the initial location of the two copresent walkers. By combining Eqs. (E2) and (E3), we obtain where q * is given by Eq. (23). With Eq. (E4), the first and second moments of IET are given by and respectively, wherê By substituting Eqs. (E5) and (E6) into Eq. (A1), we obtain Appendix F: Higher-order random walks In this section, we show that our theory holds true when the individuals move according to higher-order random walks. We focus on second-order random walks and then explain how our theory generalizes to higher orders.
In a second-order random walk, the probability with which a walker visits the next node depends on its current and last visited nodes [76,77]. Therefore, the state of a single walker is defined by a pair of nodes (m − , m), where m is the currently visited node, and m − is the node that the walker visited just before arriving in m. It should be noted that m = m − and that we distinguish between (m − , m) and (m, m − ). In other words, the state of each walker is specified by a directed edge. For example, the sequence of a walker's positions from node i to j and then to k is given by (i, j) → (j, k). Therefore, in a second-order random walk, we can regard the movement of the walkers as a first-order random walk from directed edge (i, j) to directed edge (j, k) instead of between nodes. We assume in the following text that there are N subpopulations and M undirected edges. Then, a walker moves among 2M directed edges.
The state of the system of two second-order random walkers is defined by a pair of directed edges, one for each walker. For example, if one walker is currently in subpopulation m and the other walker in subpopulation n, we denote the state of the system by ((m − , m), (n − , n)). Because the walkers are indistinguishable the system has N = 2M (2M + 1)/2 = M (2M + 1) states in total.
We denote the next state of the system by ((m − , m ), (n − , n )) and suppose that the walker at subpopulation m has moved to subpopulation m new . Therefore, the new state ((m − , m ), (n − , n )) is equal to ((m, m new ), (n − , n)). If m new = n, then the next state (i.e., ((m − , m ), (n − , n ))) is an absorbing state. Otherwise, it is a transient state.
The number of absorbing states with which the two walkers meet at node i is given by k i + ki 2 = k i (k i +1)/2, where k i is the degree of node i. Note that ((j, i), (j, i)), where j is a neighbor of i, is also a valid absorbing state. Then, the system has i absorbing states. The number of transient states is given by i . Because we have identified all the transient states, absorbing states, and state transition rules, one is able to define the N T × N T matrix Q of transition rates between transient states and the N T × N A matrix R of transition rates from a transient state to an absorbing state. Then, the theory that follows is the same as that developed in Section III C and in Appendix E.
It is straightforward to extend the same procedure to the case of directed metapopulation networks and higher-order (i.e., third-order or higher) random walks. For a third-order random walk, for instance, the state of the system is described by a pair of triples, i.e., ((m −2 , m −1 , m), (n −2 , n −1 , n)), where m is the subpopulation that the first walker currently visits, m −1 is the node that the same walker visited just before m, and m −2 the subpopulation that the walker visited just before m −1 . The definitions are analogous for (n −2 , n −1 , n).
As an example of a second-order random walk, we simulated the non-backtracking random walk. By definition, a non-backtracking random walker on an undirected unweighted network that has moved from m − to m, moves to any of the neighbors except m − with the equal probability in the next move [87,88]. The CV of IET for two non-backtracking random walkers on various networks is shown in Fig. 6. The results are similar to those for the simple random walk shown in Fig. 3. In other words, the CV is substantially larger than 1 and large when D/λ is small. In general, the CV values for the non-backtracking random walk are somewhat smaller than those for the simple random walk.