The influence of a transport process on the epidemic threshold

By generating transient encounters between individuals beyond their immediate social environment, transport can have a profound impact on the spreading of an epidemic. In this work, we consider epidemic dynamics in the presence of the transport process that gives rise to a multiplex network model. In addition to a static layer, the (multiplex) epidemic network consists of a second dynamic layer in which any two individuals are connected for the time they occupy the same site during a random walk they perform on a separate transport network. We develop a mean-field description of the stochastic network model and study the influence the transport process has on the epidemic threshold. We show that any transport process generally lowers the epidemic threshold because of the additional connections it generates. In contrast, considering also random walks of fractional order that in some sense are a more realistic model of human mobility, we find that these non-local transport dynamics raise the epidemic threshold in comparison to a classical local random walk. We also test our model on a realistic transport network (the Munich U-Bahn network), and carefully compare mean-field solutions with stochastic trajectories in a range of scenarios.


Introduction
Transport processes such as the movement through a network of airlines on a global (Hufnagel et al. 2004;Colizza et al. 2007) or a network of different public transport modes on a local scale (Balcan and Vespignani 2011; Ruan et al. 2015) play a crucial role in the spread of an epidemic, and in recent years, there has been an increasing amount of work on that topic . Events that transiently bring together people that would not otherwise meet and interact in the community impose a genuine risk and can imply a surge of infections driving the epidemic (McCloskey et al. 2020;Parnell et al. 2020;Gilat and Cole 2020).
Traditional modelling of epidemics goes back to the seminal work of Kermack et al. (1927), who considered a population divided into different compartments, the susceptible ("S"), the infected ("I"), and the recovered ("R") (or removed -in the sense that they do not contribute the epidemic spread anymore). The dynamics between the different compartments are governed by a set of ordinary differential equations that can be used to accurately describe the dynamics of many real epidemics in the last century (Kermack et al. 1927;Anderson and May 1992;Chang 2017;Dehning et al. 2020;Prodanov 2021). The classical SIR-model has been subsequently extended to account for more features of epidemic spreading (Batista et al. 2021), ranging from the addition of single compartments, e.g. for individuals that have been exposed to the contagion and are carriers but not yet infectious (Li et al. 2001;Aleta and Moreno 2020) or for individuals that are in quarantine (Horstmeyer et al. 2020), to intricate models that include different disease progressions (Romano et al. 2020;Abrams et al. 2021), intervention and containment strategies (Dashtbali and Mirzaie 2021), and age (Zhao 2020). Using data-driven methods, these models can be used to make predictions for the course of an epidemic and to assess different mitigation strategies (Dehning et al. 2020;Parino et al. 2021;Reiner 2021). A fundamental assumption of these deterministic, compartmental models is population-wide random mixing, even though in a population every individual only has a finite number of contacts (Keeling and Eames 2005). Thus, an alternative and more realistic way to model an epidemic is as a stochastic process on a network of social relationships (Pastor-Satorras et al. 2015;Kiss et al. 2017). Akin to the compartmental models, the nodes of the network, generally representing individuals in a population, can be in one of several discrete states, in the classical case, susceptible ("S"), infected ("I"), or recovered ("R"). The disease is transmitted along the links of the network, e.g. representing social relationships, in a probabilistic way. While these network models are considerably more complex than the traditional, compartmental models, depending on the topology of the underlying network, their mean-field limits can often still be reduced to differential Eq. (Kiss et al. 2017), which are more complex as they have to take the network topology into account. In particular, it has been shown that network structure can give rise to new phenomena, such as a vanishing epidemic threshold (Pastor-Satorras and Vespignani 2001).
The effect of transport processes on epidemic spreading has so far mostly been considered in meta-population models (Hufnagel et al. 2004;Colizza et al. 2007;Balcan and Vespignani 2011;Colizza et al. 2007;Brockmann and Helbing 2013;Ruan et al. 2015;Linka et al. 2020;Calvetti et al. 2020;Chang et al. 2021). In these models, one considers a set of communities in each of which the epidemic spreading is modelled by one of the classical compartmental models together with a network that connects these communities and governs the interactions between the epidemic in the different communities via a flux of the contagion. Importantly, it has been demonstrated that data-driven models of this kind can inform government policies and mitigation strategies (Chang et al. 2021). In another recent study, a framework for epidemic spreading on a network model of time-varying encounter networks was developed and used to deduce control policies for public transport (Mo et al. 2021).
In studying complex systems and dynamical systems on networks in general as well as epidemics and contact processes in particular, multilayer network structures have recently attracted an increasing amount of interest (De Domenico et al. 2013;Kivelä et al. 2014;Boccaletti et al. 2014;De Domenico et al. 2016;Bianconi 2018). As a special kind of these complex structures, multiplex networks combine several classical networks on top of the same set of nodes in a layered structure together with interlayer links between corresponding nodes in the individual layers (Boccaletti et al. 2014). In the context of epidemic dynamics, multiplex structures have e.g. been used to study the spread across multiple layers (Saumell-Mendiola et al. 2012;Ferraz de Arruda et al. 2017), the interplay between awareness of an epidemic and the epidemic itself (Granell et al. 2013), or the competition of two contagions that spread in individual layers (Darabi Sahneh and Scoglio 2014;Sanz et al. 2014). Overall, multilayer network structures do generally give rise to a much richer phenomenology than their classical counterparts (De Domenico et al. 2016).
Transport dynamics are a central aspect of this work, especially those corresponding to human mobility (Barbosa et al. 2018) and it has been recognised that the latter follows simple reproducible patterns (González et al. 2008). One key observation in many studies has been that the distribution of jump-sizes tends to be heavy-tailed (Brockmann et al. 2006;Jiang et al. 2009), which is specifically reminiscent of Lévy flights that have been a popular model for human mobility (Barbosa et al. 2018;Zaburdaev et al. 2015). At the same time, it has also been shown that popular continuous-time random walk and Lévy-flight models do not satisfactorily describe certain features of human mobility such as the tendency to preferentially return to previously and recently visited sites as well as ultraslow diffusive exploration (Song et al. 2010; Barbosa et al. 2015). These features correspond to non-Markovian dynamics and microscopic models that account for those have been shown to reproduce mobility patterns observed in data more accurately.
In this work, we introduce a model of epidemic dynamics on a network in the presence of a transport process that will give rise to a multiplex structure. In view of many epidemics for which recovery only leads to temporary immunity, we will consider SIRS-epidemic dynamics, develop a mean-field description and study it from a dynamical systems' point of view, considering its equilibria and how the transport process affects the epidemic threshold.

Model and results
In the following, we will define the model and subsequently derive a mean-field description of first and second order, i.e. at the level of individuals and pairs, respectively. We will then go on and show that the mean-field solutions provide very good approximations for direct stochastic simulation of the network dynamics. We will characterise the long-term behaviour of the dynamics towards equilibrium in terms of an epidemic threshold and study the effect transport dynamics have on this threshold.

The model
We consider an epidemic in a population, whose individuals also take part in a transport process, moving through a transport network. In addition to the links every individual has to other individuals, e.g., through social relationships in the population, the transport process will generate additional transient links along which the contagion can spread from one individual to the other.
The model the consists of two networks. On one side we have the epidemic network, a multiplex network structure of two layers, N , E c. , E t. with N the common set of nodes of the multiplex structure and E c. , E t. ⊆ N × N and on the other side the transport network (X , A) with A ⊆ X × X . The two layers of the epidemic network, (N , E c. ) and N , E t. , are referred to as the community and transport layer, respectively (Fig. 1). In terms of the topology, the community layer of the epidemic network is assumed to be a static, undirected, unweighted, simple network, while the topology of the transport layer is not fixed and subject to the transport process on the trans- Fig. 1 A model combining epidemic dynamics on a network with a transport process. The model consists of two networks, the epidemic network N , E c. , E t. , a two-layer multiplex network structure, on the left and a transport network (X , A) on the right. The two layers of the epidemic network are the community layer N , E c. and the transport layer N , E t. . While the topology of the former is assumed to be static, the one of the latter is dynamically changing. The epidemic network as a whole supports the SIRS-epidemic dynamics, where the contagion can spread from one individual to another via either links in the community or the transport layer. At the same time, the individuals of the epidemic network are assumed to move through the transport network, performing a random walk so that whenever two of them occupy the same site they will generate by a link in the transport layer of the epidemic network that persists for as long as they both occupy the same site port network (see below). The transport network is assumed to be a static, undirected, connected, weighted network, with (weighted) adjacency matrix A. Corresponding to the adjacency matrix, we define a transition matrix P is the (weighted) degreematrix. Regarding the notation, by A(x, x ) we mean the entry of A corresponding at nodes x and x ∈ X .
While the epidemic network, the multiplex structure, will support the epidemic dynamics, the transport network will support a transport process. More specifically, we will consider SIRS-epidemic dynamics across the multiplex epidemic network while at the same time letting the individuals perform a Poissonian node-centric continuoustime random walk on the transport network (Masuda et al. 2017). This process then defines the topology of the transport layer of the epidemic network. Whenever any two individuals occupy the same site in the transport network, they generate a link between them in the transport layer that persists until one of them leaves the site. As a result of that, the transport layer is a collection of complete networks on subsets of the population corresponding to the individuals at the different sites of the transport network.
To make all this more precise, with every individual n ∈ N , we associate a state X n t , H n The latter process describing SIRS-dynamics can be seen as a generalisation of the more frequently studied SIS-and SIR-dynamics. More specifically, we recover those in the limits σ → ∞ and σ → 0, respectively. We will make this precise once we have derived a mean-field description. Intuitively, however, 1 σ is the average time an individual spends in the state R once it has transitioned into this state. Thus, in the case of SIS-epidemic dynamics where upon recovery an individual becomes immediately susceptible again, this time can be thought of as arbitrarily small. In contrast, in the case of SIR-epidemic dynamics where an individual gains permanent immunity, it can be thought of as arbitrarily large. These two extremes correspond to the limits σ → ∞ and σ → 0, respectively. Moreover, it is worthwhile noting that we introduced different rates for infections via the community and the transport layer. This is to account for the fact that encounters in the community are generally different to the ones in public transport in terms of the risk of transmission of a contagion that depends on the external circumstances.
In addition to the epidemic parameters β c. , β t. , γ , and σ for the infection rates in the community and transport layer, the recovery rate and the immunity loss rate, respectively, we have introduced parameters μ for the rate at which nodes move through the transport network to control their mobility. Also, we remark that since all of the transitions are subject to homogeneous Poisson-point-processes, the overall process is Markovian.

Derivation of the model's mean-field description
If we disregard the transport process and thus assume the transport layer of the epidemic network to be static, the mean-field description of the epidemic dynamics can be deduced via standard techniques after projecting the epidemic multiplex network down onto a single layer (Kiss et al. 2017). Now, since all transitions are governed by homogeneous Poisson-point-processes, at any point in time at most one transition can occur. Therefore, we can take the mean-field description that we obtained by considering the transport process to be frozen in time and add the missing terms for a complete description by conversely considering the dynamics of the transport process in isolation with the epidemic frozen in time. Since the transport process dynamically changes only the links in the transport layer of the epidemic network and consequently does not affect the state of health of the individual nodes, it will manifest itself exclusively in the evolution of the expected number of pairs and, in general, higher-order motifs of the transport layer. In order to simplify notation, we adopt the notation of Kiss et al. (2017)   Mean-field transition diagrams for SIRS-epidemic dynamics in a multiplex network up to pairmotifs. These diagrams capture the mean-field transition rates between individual-and pair-motifs for SIRS-epidemic dynamics adapted to the case of a multiplex network where in layer λ the infection rate along S-I-links is β λ (cf. Kiss et al. 2017) and where additionally links in a transport layer (ω = t.) are subject to changes due to the transport process studied in this work. The expected number of (directed) pairs of individuals in state of health h and h connected via a link in layer λ is denoted as [h λ ∼ h ]. Importantly, in this case links whose incident nodes have the same state of health are counted twice. Similarly, the expected number of (directed) triples of individuals in states h, h and h with the first two connected via a link in layer λ and the second two via a link in layer λ is denoted as The transport process manifests itself in the additional approximate transition terms (∂ t p t 2 )[h] t [h ] t to every pair of nodes in states h and h connected via a link in the transport layer (ω = t.), where ∂ t p t = −μ p t and = 1 − P the graph Laplacian of the transport network the transport layer, where ∂ t p t = −μ p t and = 1 − P the graph Laplacian of the transport network. The full transition diagrams for a mean-field description of SIRS-epidemic dynamics from which we will later deduce the mean-field differential equations are thus given by the diagrams in Fig. 2.
Indeed, denoting the number of individuals in a state of health h occupying a site x at time t as h t (x), we have that for any τ > 0 sufficiently small, Since by assumption, individuals occupying the same site generate links between them in a transitive way, denoting the number of links in the transport layer between individuals in state of health h and h at site x at time t as {h t.
∼ h }(x), we consequently also have that Indeed, suppose individual n makes a transition X n t → X n t+τ . Then, provided H n t = h, X n t = x, and X n t+τ = x , h t (x) → h t+τ (x) = h t (x) − 1 and h t x → h t+τ x = h t x + 1 while everything else remains constant. Thus the net change of h t (x) is (δ x,X n t+τ − δ x,X n t )δ h,H n t . Since at every point in time each individual can make such a transition, we sum this across all the individuals and obtain the overall net change as claimed in Eq. (1). For links, note that {h t.
and since we obtain the expressions as claimed in Eq.
(2). At this point, it might be worthwhile mentioning that we are counting the pairs as directed links following the convention of Kiss et al. (2017). In particular, this means that any h-h-link is counted twice. Now, recalling that for the random walk every individual performs independently through the transport network we have where = 1 − P is the graph Laplacian of the transport network, we compute, by the law of total expectation, where in the last step we used that x (x, x ) = 1 − x P(x, x ) = 0 since P is a stochastic matrix. Essentially along the same lines, we also find that and E ⎡ ⎣ n,n :n =n Thus, and so that after passing to the limit τ → 0, we finally arrive at and ( 1 2 ) Since we are only interested in the overall number of h-h -links in the transport layer, we obtain that from the latter by the simply taking the sum across all sites x. Hence, While this equation is exact, it depends on the precise number of individuals at every site and in every state of health which is not suitable for a complete meanfield description at the population level. Thus, assuming the site a given individual is occupying is independent of its state of health, we make the approximation where p t (x) is the probability that an individual performing a random walk through the transport network occupies site x at time t, and find using Eqs. (11) and (13) to describe the effect of the transport process on the epidemic network. This completes the argument and explains the additional transition terms in the diagrams in Fig. 2. From these diagrams, we can now deduce the mean-field equations, which are for ω ∈ {c., t.}, where we do not explicitly highlight that the transition terms due to the transport process are only approximations in view of other approximations that are going to be introduced down the line in the form of moment closures. These equations are valid for any topology of the epidemic network. However, as indicated earlier, these equations are not closed in the sense that the evolution equations for motifs of first or second order, in turn, depend on motifs of the next higher order. This gives rise to a hierarchy of models that will only terminate at motifs of system size and at this point defeating the point of mean-field equations as a low-dimensional description of the macroscopic behaviour of a system (House and Keeling 2011). This is where moment closures enter the scene (Kuehn 2016). By expressing higher-order motifs in terms of lower-order ones, one breaks the hierarchy and obtains a closed system of ordinary differential equations.
In this work, we will only focus on first-and second-order mean-field equations, i.e. equations that have been closed at the level of individuals and at the level of pairs, respectively (Kiss et al. 2017). In order to do so, we will from now on assume that the community layer of the epidemic network is a k-regular network, i.e. that all its nodes have degree k, which greatly simplifies the analysis (for the irregular case, see Appendix).

First-order mean-field equations
In order to derive first-order mean-field equations, we will use the well-known pairclosure (Kiss et al. 2017). For that, we first note that, in addition to the community layer, also the transport layer of the epidemic network is in expectation regular with degree p t 2 |N |. Indeed, one can verify that the probability for any individual n to have degree k in the transport layer is given as |N Thus, the pair-closure is given as where κ c. ( p) = k and κ t. ( p) = p 2 |N | are the (expected) degree of the networks in the community and transport layer, respectively. Applying this moment closure to the system of Eq. (15), we obtain as the first-order mean-field equations that describe the dynamics at the level of individuals with initial conditions p 0 some probability distribution on the transport network and t is a conserved quantity. Therefore, these first-order equations are not independent and the total number of individuals as a constant of motion can be used to derive a reduced system of differential equations that equivalently describe the system by eliminating one of the equations for the susceptible, infected, and recovered individuals.

Second-order mean-field equations
Since we have described the dynamics up to the level of pair motifs, we can go one step further and also derive the second-order mean-field equations. For that, we need to close the system of Eq. (15) at the level of pairs. In order to do so, we will use the closure which is frequently used for regular networks (Kiss et al. 2017) and where as before for ω ∈ {c., t.} are the second-order mean-field equations that describe the dynamics at the level of pairs with initial conditions p 0 some probability distribution on the transport network and As it has been for the first-order mean-field equations, these equations are not independent. Indeed, there are several conservation relations with which one can derive a reduced set of differential equations that equivalently describe the full system. Besides the total number of individuals also the (expected) total number of links in both layers of the epidemic network is conserved since In addition, we also have pair-conservation so that ∼ I] t from which the claim follows via standard arguments (Kiss et al. 2017, Proposition 4.2).

SIS-and SIR-epidemic mean-field equations in a singular limit
As mentioned in the beginning, we can derive the corresponding mean-field descriptions for SIS-and SIR-epidemic dynamics by taking the limits σ → ∞ and σ → 0, respectively. While in the case of the latter, one can simply set σ = 0, the procedure is more involved in the case of the former and to actually perform the limit one might draw on results from geometric singular perturbation theory. More specifically, as σ becomes large, we observe a time-scale separation in the system with the transition R → S occurring on a fast and the remaining transitions on a slow time scale, so that in the limit σ → ∞ only the slow dynamics remain. The mean-field equations describing them are then given as (see Appendix) The corresponding first-order mean-field equations can be obtained analogously or, alternatively, by simply applying the pair-closure from before in Eq. (16) to these Fig. 3 Comparison between stochastic trajectories of the model and first-and second-order mean-field solutions. Stochastic trajectories together with their corresponding mean-field solutions for the fractions of susceptibles, infected, and recovered are shown for SIS-, SIRS-and SIR-epidemic dynamics in a population of 1000 individuals in a 10-regular network as the community layer of the epidemic network and the Munich U-Bahn network that consists of approximately 100 sites as the transport network. The epidemic parameters are set to β c. = 1 6 , β t. = 1 20 β c. , γ = 1, and, in case of SIRS-dynamics, σ = 1 5 . The mobility rate is varied, with μ = 1 in A and μ = 10 in B. In each case, the mean initial prevalence is set to 5 % and, at t = 0, the individuals are all located at a single site ("Marienplatz") on transport network. The stochastic trajectories are shown without having been time-shifted equations. We remark that the reduced slow limiting systems are relatively straightforward in our context but far more complicated singular limits have recently emerged as well in epidemic dynamics depending on which singular perturbation parameters are considered (Jardón-Kojakhmetov et al. 2021;Li et al. 2016;Schecter 2021).

Numerical comparison of first-and second-order mean-field equations with the stochastic process
Having derived mean-field equations to describe the dynamics of the stochastic process at a macroscopic level, obviously raises the question as to how well they describe the full stochastic process. For that, we computed stochastic trajectories by simulating the stochastic process using a discrete-event simulation (Kiss et al. 2017;Law 2015). A numerical comparison between stochastic trajectories and the corresponding first-and second-order mean-field solutions for SIRS-, SIS-and SIR-epidemic dynamics in a population with 1000 individuals in a 10-regular network as the community layer of the epidemic network and the Munich U-Bahn network which consists of 96 sites as the transport network shows overall good agreement between the stochastic trajectories and the mean-field solutions (Fig. 3). However, as expected, the second-order meanfield solutions are generally in better agreement with the stochastic trajectories.

The epidemic threshold in the presence of transport
Due the transport process the individuals are transiently linked to other individuals that they would not have been otherwise and as such exposed to a higher risk of infection. The way the transport process increases the infection pressure can be seen most directly Fig. 4 The presence of the transport process lowers the epidemic threshold. The equilibrium prevalence near the epidemic threshold between disease-free and endemic state is shown as a function of the relative strength of the infection rate in the transport layer β t. β c. and the basic reproduction number in the absence of transport β c. k γ (in first-order approximation) for the first-and second-order mean-field description in A and B, respectively, as well as for the mean of stochastic simulations in C. For simplicity, this is done for SIS-epidemic dynamics (for SIRS-epidemic dynamics one can obtain qualitatively similar plots) in a population of 1000 individuals in a 10-regular network as the community layer of the epidemic network and the Munich U-Bahn network consisting of approximately 100 sites as the transport network. The recovery rate and the mobility rate are set to γ = 1 and μ = 10, respectively. As the infection rate in the transport layer β t. is increased, an increasing number of infections occurs via the transport layer which in turn has as a consequence that the epidemic threshold is lower in the first-order mean-field models. In fact, the effective infection rate in the presence of the transport process turns out to be λ β λ κ λ ( p t ) In other words, the transport process effectively amounts to β t. β c. p t 2 |N | ≥ β t. β c.

|N | |X |
additional contacts for the average individual. These additional, transient contacts can considerably alter the fate of the epidemic. If, e.g. in case of SIS-dynamics, the community layer alone is just sufficiently sparsely connected so that the epidemic cannot be sustained and will die out eventually, the presence of the transport process can lead to an endemic state (Fig. 4). More precisely, consider the first-order mean-field model in Eq. (17) and let ( 2 2 ) Then, the system undergoes a transcritical bifurcation when χ ( p ∞ ) = 1 where p ∞ is the unique solution to the equation corresponding to the disease-free and the endemic state, respectively. Indeed, due to the strongly connected transport network there is a unique solution p ∞ of ∂ t p t = −μ p t = 0 with p ∞ ≥ 0 and p ∞ = 1, the equilibrium solution of the random walk on the transport network. From there, a proof of the above statement proceeds entirely analogously to the case of standard SIRS-epidemic dynamics (Kiss et al. 2017).
In particular, since χ ( p ∞ ) ≥ χ ( p ∞ ) | β t. =0 , in the presence of the transport process and a non-vanishing infection rate in the transport layer, the epidemic threshold towards an endemic state is lower (Fig. 4).
In terms of mitigation strategies for an epidemic, this means that, without comprehensive testing and subsequent quarantining of infected individuals as well as additional hygiene measures on public transport, people that regularly participate in public transport, e.g. by commuting to work, would have to restrict their personal contacts even more drastically than people that do generally avoid public transport and by that keep their number of effective contacts low. Alternatively, the density of people in public transport has to be kept low, which would result in a lower infection rate and similarly keep the overall number of effective contacts low.
Finally, let us also remark that the corresponding results for SIR-and SIS-epidemic dynamics can again be obtained by considering the limits σ → 0 and σ → ∞, respectively. Unfortunately, a full analytical description of the equilibrium state of the second-order mean-field equations in terms of closed-form expressions is virtually impossible in most cases. Even in the simplest case of the SIS-epidemic dynamics and considering the fully reduced system after making use of all the conservation relations, a description of the endemic state turns out to be very complicated.

Non-local, fractional transport dynamics
So far, we have considered transport under the dynamics of the standard graph Laplacian. In this case, the individuals performing the random walk can only move to sites in the immediate neighbourhood of the graph. However, empirical studies suggest that among other things human mobility patterns are characterised by heavy-tailed distributed jump-lengths (Brockmann et al. 2006) reminiscent of Lévy flights on the network (Zaburdaev et al. 2015). In the following, we will consider a generalisation of our model towards fractional dynamics on the transport network, which are known to give rise to heavy-tailed jump-length distributions (Riascos and Mateos 2014;Michelitsch et al. 2017;Benzi et al. 2020;Michelitsch et al. 2019).
In defining the fractional dynamics for an exponent 0 < α ≤ 1, we follow the approach of Benzi et al. (2020). We consider the unnormalised Laplacian L := K − A which is symmetric and non-negative. As such, we can define its fractional power in the usual way: Given the spectral decomposition L = λ λ λ , we have L α = λ λ α λ . The fractional degree-matrix is then given as x) and we obtain the fractional Laplacian (α) = K (α)−1 L α so that the fractional transition matrix is given as P (α) = 1 − (α) . With the fractional transition matrix or Laplacian, we can define our model entirely analogously and also the results we have obtained so far carry over.
While for α = 1, we trivially recover the dynamics of a classical random walk, for 0 < α < 1, we obtain superdiffusive behaviour in the transport process as well as an algebraic decay in the jump-length distribution (Benzi et al. 2020;Riascos and Mateos 2014). One of the consequences that come with the introduction of the fractional dynamics is that it eventually raises the epidemic threshold (Fig. 5). Under fractional dynamics, individuals spread more evenly on the transport network facilitated by long- Fig. 5 Non-local, fractional dynamics in the transport process raise the epidemic threshold. The equilibrium prevalence near the epidemic threshold between disease-free and endemic state is shown as a function of the fractional exponent α and the basic reproduction number in the absence of transport β c. k γ (in first-order approximation) for the first-and second-order mean-field description in A and B, respectively, as well as for the mean of stochastic simulations in C. For simplicity, this is done for SIS-epidemic dynamics (for SIRS-epidemic dynamics one can obtain qualitatively similar plots) in a population of 1000 individuals in a 10-regular network as the community layer of the epidemic network and the Munich U-Bahn network consisting of approximately 100 sites as the transport network. The infection rate in the transport layer, the recovery rate, and the mobility rate are set to β t. = 1 20 β c. , γ = 1, and μ = 10, respectively. As the fractional exponent α is lowered, the individuals spread more evenly across the transport network and do not form clusters which in turn has a consequence that the epidemic threshold is higher distance jumps. As such, the equilibrium distribution of the transport process on the network approaches a uniform distribution as the fractional exponent becomes small.
Indeed, it is well-known that the equilibrium distribution for the transport dynamics is given as p (Benzi et al. 2020), where k (α) (x) = λ λ α λ (x, x) for any 0 < α ≤ 1. On one hand, considering the case α = 1, we immediately see that in equilibrium there is cluster formation at well-connected sites. On the other hand, if we let α tend to 0 these clusters dissolve. More specifically, for |X | where we have used that λ λ = 1 and that the eigenvector corresponding to the eigenvalue 0 is constant. Hence, as α approaches 0, the fractional degrees at every site converge to a constant and, consequently, the equilibrium distribution is continuously approaching a uniform distribution. In turn, since the term p 2 for any probability distribution p is minimised by the uniform distribution, we have that ∞ ) for any 0 < α ≤ 1. However, while one might assume that the dependence in α is monotonic, so that the epidemic threshold towards an endemic state is higher the more α is decreased, this turns out to be false in general. In fact, it depends on the topology of the network and there are counterexamples of networks where the epidemic threshold decreases before it finally increases (see Appendix).
Overall, the effect the fractional dynamics have on the epidemic threshold is rather subtle depending on the topology of the transport network. In fact, the more regular it is, the smaller it is. Conversely, the more irregular it is, the more apparent the difference in the epidemic threshold for fractional and standard random walk (α = 1) dynamics is.

Temporary instead of permanent participation in the transport-epidemic dynamics
The transport model we have considered so far has the apparent flaw that every individual is permanently moving through the transport network. In contrast, more realistically, individuals would enter the transport network at one site and leave it again at another site and, therefore, only temporarily move through the transport network participating in the transport-induced epidemic dynamics.
In order to incorporate that into the present model, we extend a given transport network with a set of sites, accessible by some or all other sites, that will be regarded as situated beyond the transport network in the sense that whenever individuals occupy these sites they are not aware of each other and as such do not generate a transient link in the transport layer. This change then shows up in the way the distribution on the extended transport network enters the mean-field description. Specifically, if (X , A) is the extended transport network with X = X ∪ and the aforementioned set of sites beyond the transport network, then any term p t 2 showing up in the equations is replaced by the term p t 1 c 2 = x∈X \ p t (x) 2 .
Indeed, if is now the graph Laplacian of the extended transport network and starting from Eq. (12), we obtain the overall number of h-h -links in the transport layer in this case by only summing across the sites in X as opposed to also those in , which yields an exact equation analogous to Eq. (13). After making the same approximation as before, which takes the place of Eq. (14) for the extended transport network. Thus, instead of p t 2 we have p t 1 c 2 for an extended transport network with void sites . Moreover, note that this reduces to the previous case where there are no void sites and = ∅. By a similar argument where one restricts the sum across the transport network's sites to exclude those in the void, one can show that in expectation the degree of an individual's node in the transport layer is analogously p t 1 c 2 |N | instead of p t 2 |N |. Overall, the introduction of void sites binds some of the density of individuals on the transport network, reduces the number of links in the transport layer and by that the infectious pressure. Consequently, this raises the epidemic threshold but otherwise does not qualitatively change the behaviour of the dynamics.

Numerical investigation of non-equilibrium dynamics
So far, we have exclusively focused on the equilibrium of the dynamics. In the following, we will discuss two other aspects regarding dynamics that do not necessarily Fig. 6 Comparison between stochastic trajectories of the model and first-and second-order mean-field solutions for a time-varying transport network. Stochastic trajectories together with their corresponding mean-field solutions for the fractions of susceptibles, infected, and recovered are shown for SIS-, SIRSand SIR-epidemic dynamics in a population of 1000 individuals in a 10-regular network as the community layer of the epidemic network and the Munich U-Bahn network that consists of approximately 100 sites as the transport network. The dynamics switch between the transport network as considered before and one with an attracting core that arises from the former after resolving any link as two directed links in either direction and then deleting any link that strictly increases the number of steps required to reach the core and by that creating a time-varying transport network. The epidemic parameters are set to β c. = 1 6 , β t. = 1 20 β c. , γ = 1, and, in case of SIRS-dynamics, σ = 1 5 . The mobility rate is varied, with μ = 1 in A and μ = 10 in B. In each case, the mean initial prevalence is set to 5 % and, at t = 0, the individuals are all located at a single site ("Marienplatz") on the transport network. The stochastic trajectories are shown without having been time-shifted. The periods when the transport network was attracting are shaded in grey admit an equilibrium on one hand and dynamics before reaching equilibrium on the other hand.

Time-dependent transport networks
Besides static, undirected transport networks that we have considered so far, it is also interesting to have the network dynamically changing. Inspired by the dynamics induced by the daily commute of people from their home to e.g. the central business district, one might consider a network that is temporarily attracting, drawing a random walker into a predefined core. One way to model this is to consider, in addition to a base network, a second (directed) transport network, derived from the former, where the strength of the links is biased so that a random walk will concentrate a large proportion of its mass at a small set of nodes, an attracting core. Given a function f : [0, ∞[→ [0, 1], we can then define transport dynamics on a dynamically changing network, governed by a timedependent Laplacian˜ t = (1 − f (t)) + f (t) with and the Laplacians of the base network and the attracting network, respectively. Depending on the modulation function f one can model different dynamics. In the case of periodic dynamics such as the initially mentioned daily commute, an obvious choice would be an oscillating function. In contrast, a single event temporarily drawing people to a particular site, can be realised by a function that switches only once between the two networks.
The mean-field description that we have derived readily applies to this case and for a numerical comparison between stochastic trajectories and first-and second-order mean-field solutions we consider again as a transport network the Munich U-Bahn network with the sites in the city centre as an attracting core. The attracting network arises from the base model after resolving every link between any two sites as two directed links in either direction and then deleting any link that strictly increases the number of steps required to reach the core. As before, we find overall good agreement between the stochastic trajectories and the corresponding first-and second-order meanfield solutions (Fig. 6).
However, it should come as no surprise that in general, the transport process and with it the whole system does not approach an equilibrium anymore. Moreover, it is not even possible to derive bounds on the epidemic dynamics in terms of the dynamics on one or the other transport network.

Surges in infections as a consequence of an accumulation of individuals at a single site
While the rate μ at which the individuals move through the transport network does not influence the equilibrium of the whole system, it can hugely influence the trajectory towards the equilibrium. Specifically, this is most profound when a considerable proportion of individuals accumulate at a single site in the transport network. In that case, the accumulation site has to be evacuated sufficiently fast, i.e. μ has to be large, in order for this to not lead to a surge in infections and a rise in prevalence. Starting with a relatively low prevalence of the epidemic and all the nodes accumulated at a particular site, we may consider the trajectories towards equilibrium under SIS-dynamics for different mobility rates. When the mobility rate is sufficiently high, the accumulation at a single site does not drastically alter the trajectory and the prevalence simply rises to the endemic equilibrium. However, for low mobility rates, one observes the prevalence rising beyond the endemic equilibrium and only slowly recovering to equilibrium afterwards (Fig. 7A). Moreover, the threshold mobility rate for these two extremes critically depends on the initial accumulation site in the transport network. In particular, it turns out that it is not necessarily the degree of the site but rather its centrality in the network that determines this threshold mobility rate. The more central a site is, the lower is this mobility rate (Fig. 7B).
Again, in terms of mitigation strategies, this suggests that for mass events there is a critical time for people to dwell at a particular site. In addition, this time depends on the location. Therefore, for such events to be as safe as possible and not imply a huge number of infections, they should be held at well-connected sites.

Discussion
In this work, we have introduced a network model that combines epidemic dynamics with a transport process in a multiplex network structure. The latter consists of two layers, a static one and one that is dynamically adapting in response to a transport process on a separate network. We have derived the model's mean-field description and analysed it from a dynamical systems' point of view, characterising its long-term behaviour and specifically focusing on the epidemic threshold between a disease- Fig. 7 Surge of infections due to accumulation of individuals at a single site. While the mobility rate does not affect the equilibrium prevalence, it can alter the epidemic in the short term. For that consider SISepidemic dynamics in a population of 1000 individuals in a 10-regular network as the community layer of the epidemic network. As before, the Munich U-Bahn network consisting of approximately 100 sites is taken as the transport network. The epidemic parameters are set to β c. = 1 6 , β t. = 1 20 β c. , γ = 1. In A, the trajectories towards equilibrium for different mobility rates are shown when every individual begins at the same initial site. For high mobility rates, the prevalence simply rises from the initial prevalence to equilibrium. However, there is a threshold mobility rate, μ * , under which prevalence rises beyond the equilibrium value. In B, it is shown that this threshold mobility rate depends on the initial site and its neighbourhood in the transport network. Specifically, for sites that are central and well-connected this threshold is lower than for more decentral ones free and an endemic equilibrium state. We have shown that the transport process induces additional ways for the contagion to spread and that way lowers the epidemic threshold. We then generalised the transport process to fractional and, in particular, non-local dynamics and have shown that this, conversely, raises the epidemic threshold. The extent of which depends on the topology of the transport network and is more pronounced the more irregular the transport network is. Although we have derived the mean-field description up to second order, our analysis rests upon only the firstorder equations. The two-layer multiplex structure of the epidemic network essentially doubles the number of equations required for a second-order mean-field description making it hard to find an analytical description of the endemic equilibrium. Yet, we have shown that first-and second-order numerical solutions are in good agreement so that even the former approximates the dynamics sufficiently well most of the time.
We have presented a generic, conceptual model for epidemic dynamics in the presence of a transport process. In the interest of analytical tractability, we have considered a continuous-time random walk that every individual performs independently on the transport network. However, as has been pointed out elsewhere, these processes do not capture all aspects of human mobility, in particular, since many of those are inherently non-Markovian. This includes behavioural features such as preferential return and recency (Song et al. 2010;Barbosa et al. 2015) as well as more statistical features such as heavy-tailed waiting times (Brockmann et al. 2006). Intuitively, though, we would expect that these features drive the epidemic by facilitating cluster formations and on average longer waiting times at any site in the transport network.
For the analytical treatment, we have considered static transport networks. However, in reality, one should expect them to dynamically change throughout the day. If we think e.g. about the daily commute, it is obvious that in the morning and evening the flux of individuals moving through the transport network is slightly directed inwards to and outwards from the city centre, respectively. Similarly, events can draw a proportion of individuals towards a certain site. In our model, as a straightforward generalisation, dynamic transport networks are reflected in a dynamic (fractional) Laplacian governing the transport dynamics. The derivation of the mean-field descriptions then carries through analogously. In general though, with a dynamic transport network, the system will not approach an equilibrium anymore and it seems very difficult to prove anything analytically about the epidemic threshold in this case. However, it is conceivable for any such dynamics to temporarily induce cluster formations and therefore again potentially lead to surges in infections. As we have shown, the latter crucially depends on the relation between the characteristic time scales of the transport process and epidemic spreading as well as the accumulation site and the topology of the transport network.
In the context of this work, one potentially interesting idea for future research would be to imagine in addition to the epidemic network having also a multiplex structure on the side of the transport network. Such a structure could then account for different modes of transport so that every individual can switch between the different layers at certain sites or for individual transport networks where every individual moves within its own layer. Another idea would be to consider an adaptive transport process, e.g. by restricting the mobility of individuals in a certain state of health or by avoiding sites with a high local prevalence of the epidemic.
Overall, given the importance transport processes have for the epidemic spreading, we do hope this work is going to inspire further mathematical modelling investigations into the intricate interplay between transport and contact processes that may then in turn inform policies for the mitigation of an epidemic and the design of transport networks in the future.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A.1 The singular limit → ∞ in the second-order mean-field equations of SIRS-dynamics
In this section, we will demonstrate how to derive the second-order mean-field equations of SIS-epidemic dynamics in the limit σ → ∞ from the corresponding second-order mean-field equations of SIRS-epidemic dynamics. As already mentioned, the techniques to perform this limit are provided by geometric singular perturbation theory; see e.g. Kuehn (2015); Wechselberger (2020) and references therein. When σ becomes large, we eventually observe a time-scale separation in the dynamics with transitions S → I and I → R on a comparatively slow and the transition R → S on a fast time-scale. As such, the idea is to consider the problem on the fast time-scale first, then to determine the critical manifold where the fast dynamics come to a halt, and finally deduce the residual slow dynamics on this critical manifold, which are the only ones remaining in the limit σ → ∞.
In order to do that, we start by rescaling time in (19) to the fast time-scale via the transformation τ : t → σ t. Then the SIRS-second-order mean-field equations read which can be written compactly as for some particular choice of G.
Since N has full column rank, the requirement N f ( τ ) = 0 is only fulfilled when f ( τ ) = 0. Thus, the critical manifold S of this system is given as the zero level-set of f , i.e. S = {ξ : f (ξ ) = 0}. In particular, since spec D f ( ) N = {−1, −2}, this critical manifold is normally hyperbolic and attracting.
In the singular limit σ → ∞, the slow dynamics of the entire system are contained entirely in the critical manifold. On their respective time-scale, they are known to be concisely given as (Wechselberger 2020, Lemma 3.4) Thus, finally, after performing the algebra, this can be expanded to and, except for the additional terms we have as a result of the transport process, these equations are essentially the same as the ones that have been derived also elsewhere for SIS-epidemic dynamics (Kiss et al. 2017).

A.2 Monotonicity of the epidemic threshold in the fractional exponent
In this section, we will discuss the question of monotonicity of the epidemic threshold as one varies the fractional exponent of the transport dynamics. This essentially comes down to the question of whether the function α → p ∞ the equilibrium distribution of the fractional random walk with exponent 0 < α ≤ 1 on some network, which we will assume to be connected, is monotonic. As mentioned in the main text, there is no guarantee for this to be monotonic in general. Specifically, as we will argue, it depends on the topology of the network. In the following, we will first consider networks with a strong block structure and show that these networks violate monotonicity. In contrast to that, we will then consider star networks for which one can explicitly prove monotonicity.
Before moving on, recall that, if L is the Laplacian of some graph with spectral decomposition L = λ λ λ where λ is the orthogonal rank-1 spectral projection onto the eigenspace for the eigenvalue λ, the fractional Laplacian is given as L α = λ λ α λ . The equilibrium distribution of the corresponding random walk is then is the fractional degree at site x. Moreover, as α tends to 0, this distribution becomes uniform and we will denote the limiting distribution as p where, in particular, p (α) ∞ 2 ≥ p (0) ∞ 2 since the uniform distribution among all discrete probability distributions is the unique distribution minimising the 2-norm. Therefore, if α → p (α) ∞ 2 is monotonic, it is necessarily non-decreasing. Using that ∂ α L α = 1 α L α ln L α , we have that where tr(L α ) > 0 so that monotonicity follows if

Non-monotonic behaviour for networks with a block structure
As it will turn out, networks with a strong block structure produce non-monotonic behaviour and as such violate inequality (A7). Before presenting an analytical argument, let us first consider a numerical approach that will motivate considering this particular network topology.
For that, we recall that for any Laplacian L and 0 < α ≤ 1, L α is again a Laplacian matrix (Michelitsch et al. 2019), so that for inequality (A7) to hold it is equivalent to show that for any Laplacian matrix L tr(L) In fact, without loss of generality, it is enough to consider Laplacian matrices L with tr(L) = 1 only. Indeed, suppose the above has been established for such ones, for arbitrary L, letL = 1 tr(L) L. In that case,L lnL = 1 tr(L) L ln L − ln tr(L) tr(L) L and tr(L lnL) = 1 tr(L) tr(L ln L) − ln tr(L). With that, on one hand and on the other hand which implies the assertion. Now, considering only Laplacian matrices L with tr(L) = 1, has the distinctive advantage that by that one can separate the spectra of L and L ln L. Indeed, if tr(L) = 1, spec(L) ⊂ [0, 1] whereas spec(L ln L) ⊂ [−e −1 , 0] so that L and L ln L are positive and negative semidefinite, respectively. Moreover, it allows for a systematical with the diagonal elements set appropriately so that column and row sums are 0. Hence, through sampling the respective standard simplex uniformly, we can numerically test the validity of inequality (A8) and consequently monotonicity for all possible networks. The first counterexample we found that way was on the standard 9-simplex, i.e. in dimension 5, while we were not able to find any counterexamples in lower dimensions. For the matrix Thus, inequality (A8) and therefore also inequality (A7) at α = 1 are not satisfied, so that for the corresponding network α → p This counterexample as well as the others we have subsequently found have in common that the underlying networks exhibit a block structure while the degree distribution is approximately uniform. In order to make this precise, we will construct a prototypical example of such a network and demonstrate that these networks do indeed produce non-monotonic behaviour.
Let d 1 , d 2 ≥ 2 and suppose that d 1 > d 2 . Furthermore, let δ = d 1 −d 2 d 1 +d 2 −2 and define the (d 1 , d 2 )-block adjacency matrix This constitutes an undirected, weighted network with d 1 +d 2 nodes. It is connected and the corresponding node degrees in the first an second block are (1 − δ) (d 1 − 1) + d 2 and (1 + δ) (d 2 − 1) + d 1 , respectively, so that they coincide up to an order of since (1 − δ) (d 1 − 1) = (1 + δ) (d 2 − 1). Moreover, we can write the corresponding Laplacian matrix as where L K d and L K d,d denote the Laplacian matrices corresponding to graphs K d and K d,d , respectively. The former denotes the complete graph on d vertices while the latter denotes the complete bipartite graph on two sets of vertices with size d and d . Now, for every exponent 0 < α ≤ 1 the map is well-defined and continuous for any self-adjoint matrix L and thus in particular for L a Laplacian matrix. Indeed, the mappings L → L α and L → L α ln L α via the functional calculus (λ → λ α and λ → λ ln λ can be continuously defined on R by extending them with 0 beyond their immediate domain of definition so that one can approximate them via polynomials on a sufficiently large compact set and proceed via a triangle-inequality estimate to show continuity), the trace, the projections L α → L α (x, x) and L α ln L α → (L α ln L α )(x, x), and, finally, the arithmetic operations are all continuous and therefore also their composition.
By the initial discussion, this shows that there exists an exponent α at which α → p (α) ∞ 2 is strictly decreasing contradicting the fact that if it were monotonic it is necessarily non-decreasing and thus demonstrates that indeed monotonicity is violated for networks with a strong block structure (Fig. 8B).

Monotonic behaviour for star networks
In contrast to the networks with a block structure, numerical evidence suggests that many other networks indeed produce monotonic behaviour. Here, we will explicitly consider star networks.
Let d ≥ 1 and consider a star-graph of d + 1 vertices with adjacency matrix The Laplacian spectrum of such a graph is given as {0, 1, d + 1} so that from the corresponding spectral decomposition of the Laplacian matrix we find that for any function f where 0 and d+1 are the orthogonal rank-1 projections onto the linear spaces spanned by (1, . . . 1) and (d, −1, . . . − 1), the eigenvectors for the eigenvalues 0 and d + 1, respectively.
, we then have that which establishes monotonicity, since, again by the initial discussion, this implies that α → p (α) ∞ 2 is non-decreasing.

A.3 The mean-field equations for irregular networks in the community layer
When deriving the mean-field equations up to second order, we have assumed that the community layer of the epidemic network is regular, i.e. its nodes all have the same degree. In this section, we will outline how to derive mean-field models in the case when this assumption fails and the community layer is irregular. This can be seen as a generalisation of what we presented so far, however, the analysis also proves slightly more difficult. As before, the mean-field description of the epidemic dynamics with the transport dynamics frozen in time can be deduced via standard techniques (Kiss et al. 2017). Rather than the overall expected number of individuals in a given state of health, one considers the expected number of individuals in a given state of health as well as with a certain degree (in the community layer of the epidemic network) together with the corresponding higher-order motifs. Then, since in the transport process we do not distinguish between the degree of the individuals and their movement is independent of each other, the coupling of the epidemic and the transport dynamics happens entirely analogously to the regular case. respectively. In addition, in order to simplify notation, we set [h k As can be shown along the same lines as before, the transport process amounts to transition terms (∂ t p t 2 )[h k ] t [h k ] t to the expected number of pairs of individuals in state of health h and h and degree k and k in the transport layer. Thus, the mean-field equations up to the level of pairs are given as for ω ∈ {c., t.} and k and k all possible degree values. For regular networks, this system of equations is the same as the one we have seen earlier in Eq. (15).

Moment-closures for first-and second-order mean-field equations
In order to close this system at the level of pairs and triples to arrive at the first-and second-order mean-field equations, one can apply the following closures. For the firstorder equations, the closure relation is given as [S k c. ∼ I] t ≈ p t 2 [S k ] t [I] t (A24) Fig. 9 Comparison between stochastic trajectories of the model and first-and second-order mean-field solutions. Stochastic trajectories together with their corresponding mean-field solutions for the fractions of susceptibles, infected, and recovered are shown for SIS-, SIRS-and SIR-epidemic dynamics in a population of 1000 individuals in a network with 900 nodes of degree 3 and 100 nodes of degree 73 as the community layer of the epidemic network and the Munich U-Bahn network that consists of approximately 100 sites as the transport network. The epidemic parameters are set to β c. = 1 6 , β t. = 1 20 β c. , γ = 1, and, in case of SIRS-dynamics, σ = 1 5 . The mobility rate is varied, with μ = 1 in A and μ = 10 in B. In each case, the mean initial prevalence is set to 5 % and, at t = 0, the individuals are all located at a single site ("Marienplatz") on transport network. The stochastic trajectories are shown without having been time-shifted and for the second-order equations, it is given as where κ c. k ( p) = k and κ t. k ( p) = p 2 |N |. For irregular epidemic networks, these closures are well-known (Kiss et al. 2017) and their generalisation to the multilayer network we are considering here follows along the same lines as before, using that the transport layer is in expectation regular. Note that we denote the moments of the degree distribution as K z = k |N k | |N | k z where N k denotes the set of individuals with degree k in the community layer.
Applying these closure relations, specifically the first-order mean-field equations are given as In order to compute the remaining eigenvalues, note that for an m ×m square-matrix W with W i,i = q i (a i a i + c) where q i ≥ 0 for every i, i q i = 1, and c > 0, which yields the assertion upon computing the roots of this polynomial. Now, observe that the matrix diag(|N k |) k B appearing in the Jacobian above is of a similar form as the one of the matrix W in the statement above. In fact, diag(|N k |) k B k,k = β c. K |N k | |N | kk + β t. β c. K |N | p ∞ 2 with k |N k | |N | = 1. With that and continuing from Eq. (A28), we find that These eigenvalues are all real-valued and they are strictly negative as long as