Epidemic Threshold in Continuous-Time Evolving Networks

TL;DR

This paper introduces a theoretical framework to accurately compute epidemic thresholds in continuous-time evolving networks, unifying various models and overcoming limitations of previous approximation methods.

Contribution

It presents the infection propagator approach and the weak commutation condition, enabling a unified analysis of epidemic thresholds in continuous-time network models.

Findings

Provides a coherent connection between discrete and continuous time models.

Introduces the infection propagator approach for epidemic threshold calculation.

Unifies different network models under a single formalism.

Abstract

Current understanding of the critical outbreak condition on temporal networks relies on approximations (time scale separation, discretization) that may bias the results. We propose a theoretical framework to compute the epidemic threshold in continuous time through the infection propagator approach. We introduce the {\em weak commutation} condition allowing the interpretation of annealed networks, activity-driven networks, and time scale separation into one formalism. Our work provides a coherent connection between discrete and continuous time representations applicable to realistic scenarios.

Equations

$$p_{h+1,i}=$$

$$+p_{h,i}\left(1-\mu \Delta t\right),$$

$$P(T_{step})=\prod _{k=1}^{T_{step}}\left[1-\mu \Delta t+\lambda \Delta t{A}_k\right].$$

$$\frac{P(t+\Delta t)-P(t)}{\Delta t}=P(t)\left[-\mu +\lambda A(t+\Delta t)\right],$$

$$\dot{P}(t)=P(t)\left[-\mu +\lambda A(t)\right],$$

$${\dot{P}}_{ij}(t)=\sum _k\left[\lambda A_{kj}(t)-\mu {\delta }_{kj}\right]P_{ik}(t).$$

$${\dot{p}}_i(t)=\lambda \left[1-p_i(t)\right]\sum _j{A}_{ij}(t)p_j(t)-\mu p_i(t).$$

$${\dot{p}}_i(t)=\sum _j\left[\lambda A_{ij}(t)-\mu {\delta }_{ij}\right]p_j(t).$$

$$P(t)=1+\sum _{j>0}{\mu }^j{P}^{(j)}(t),$$

$$P^{(j)}(T)=$$

$$\cdot [\gamma A(x_{j-1})-1]\cdots [\gamma A(x_1)-1].$$

$$P(t)=\mathcal{T}\exp {\int }_0^{t}dx\left[-\mu +\lambda A(x)\right].$$

$$[A(t),{\int }_0^{t}dxA(x)]=0,\forall t\in [0,T],$$

$$P(T)=e^{T\left[-\mu +\lambda ⟨A⟩\right]},$$

$${\lambda }_c=\frac{\mu }{\rho [⟨A⟩]}.$$

$${\lambda }_c\delta =\frac{\mu }{m\left(⟨a⟩+\sqrt{⟨a^2⟩}\right)},$$

Full text

**Epidemic Threshold in

Continuous-Time Evolving Networks

** Eugenio Valdano1,$\dagger$, Michele Re Fiorentin2, Chiara Poletto1,

and Vittoria Colizza1,3

1Sorbonne Universités, UPMC Univ Paris 06, INSERM, Institut Pierre Louis d’Épidémiologie et de Santé Publique (IPLESP UMRS 1136), F75012, Paris, France. 2Center for Sustainable Future Technologies, CSFT@PoliTo, Istituto Italiano di Tecnologia, corso Trento 21, 10129 Torino, Italy. 3ISI Foundation, Torino, Italy. $\dagger$Current affiliation: Department d’Enginyeria Informàtica i Matemàtiques, Universitat Rovira i Virgili, 43007 Tarragona, Spain.

**Abstract

Current understanding of the critical outbreak condition on temporal networks relies on approximations (time scale separation, discretization) that may bias the results. We propose a theoretical framework to compute the epidemic threshold in continuous time through the infection propagator approach. We introduce the weak commutation condition allowing the interpretation of annealed networks, activity-driven networks, and time scale separation into one formalism. Our work provides a coherent connection between discrete and continuous time representations applicable to realistic scenarios. **

Contagion processes, such as the spread of diseases, information, or innovations [1, 2, 3, 4, 5], share a common theoretical framework coupling the underlying population contact structure with contagion features to provide an understanding of the resulting spectrum of emerging collective behaviors [6]. A common keystone property is the presence of a threshold behavior defining the transition between a macroscopic-level spreading regime and one characterized by a null or negligibly small contagion of individuals. Known as the epidemic threshold in the realm of infectious disease dynamics [1], the concept is analogous to the phase transition in non-equilibrium physical systems [7, 8], and is also central in social contagion processes [5, 9, 10, 11, 12, 13].

A vast array of theoretical results characterize the epidemic threshold [14], mainly under the limiting assumptions of quenched and annealed networks [15, 4, 16, 17, 18], i.e., when the time scale of the network evolution is much slower or much faster, respectively, than the dynamical process. The recent availability of data on time-resolved contacts of epidemic relevance [19] has, however challenged the time scale separation, showing it may introduce important biases in the description of the epidemic spread [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 19] and in the characterization of the transition behavior [34, 35, 31, 36, 37]. Departing from traditional approximations, few novel approaches are now available that derive the epidemic threshold constrained to specific contexts of generative models of temporal networks [38, 39, 22, 40, 35, 41, 32] or considering generic discrete-time evolving contact patterns [42, 43, 44]. In particular, the recently introduced infection propagator approach [43, 44] is based on a matrix encoding the probabilities of transmission of the infective agent along time-respecting paths in the network. Its spectrum allows the computation of the epidemic threshold at any given time scale and for an arbitrary discrete-time temporal network. Leveraging an original mapping of the temporal network and epidemic spread in terms of a multilayer structure, the approach is valid in the discrete representation only, similarly to previous methods [17, 18, 35].

Meanwhile, a large interest in the study of continuously evolving temporal networks has developed, introducing novel representations [45, 20, 27, 19] and proposing optimal discretization schemes [46, 47, 44] that may however be inaccurate close to the critical conditions  [48]. Most importantly, the two representations – continuous and discrete – of a temporal network remain disjointed in current network epidemiology. A discrete-time evolving network is indeed a multilayer object interpretable as a tensor in a linear algebraic representation [49]. This is clearly no longer applicable when time is continuous, as it cannot be expressed in the form of successive layers. Hence, a coherent theoretical framework to bridge the gap between the two representations is still missing.

In this Letter, we address this issue by analytically deriving the infection propagator in continuous time. Formally, we show that the dichotomy discrete time – continuous time translates into the separation between a linear algebraic approach and a differential one, and that the latter can be derived as the structural limit of the former. Our approach yields a solution for the threshold of epidemics spreading on generic continuously evolving networks, and a closed form under a specific condition that is then validated through numerical simulations. In addition, the proposed novel perspective allows us to cast an important set of network classes into one single rigorous and comprehensive mathematical definition, including annealed [50, 4, 51] and activity-driven [35, 52] networks, widely used in both methodological and applied research.

Let us consider a susceptible-infected-susceptible (SIS) epidemic model unfolding on a continuously evolving temporal network of $N$ nodes. The SIS model constitutes a basic paradigm for the description of epidemics with reinfection [1]. Infectious individuals (I) can propagate the contagion to susceptible neighbors (S) with rate $\lambda$, and recover to the S state with rate $\mu$. The temporal network is described by the adjacency matrix $A(t)$, with $t\in\left[0,T\right]$. We consider a discretized version of the system by sampling $A(t)$ at discrete time steps of length $\Delta t$ (Fig. 1). This yields a finite sequence of adjacency matrices $\{A_{1},A_{2},\cdots,A_{T_{step}}\}$, where $T_{step}=\left\lfloor T/\Delta t\right\rfloor$, and $A_{h}=A(h\Delta t)$. The sequence approximates the original continuous-time network with increasing accuracy as $\Delta t$ decreases. We describe the SIS dynamics on this discrete sequence of static networks as a discrete-time Markov chain [17, 18]:

[TABLE]

where $p_{h,\,i}$ is the probability that a node $i$ is in the infectious state at time step $h$, and $\mu\Delta t$ ($\lambda\Delta t$) is the probability that a node recovers (transmits the infection) during a time step $\Delta t$, for sufficiently small $\Delta t$.

By mapping the system into a multilayer structure encoding both network evolution and diffusion dynamics, the infection propagator approach derives the epidemic threshold as the solution of the equation $\rho\left(P(T_{step})\right)=1$ [43, 44], where $\rho$ is the spectral radius of the following matrix:

[TABLE]

The generic element $P_{ij}(T_{step})$ represents the probability that the infection can propagate from node $i$ at time step $1$ to node $j$ at time step $T_{step}$, when $\lambda$ is close to $\lambda_{c}$ and within the quenched mean-field approximation (locally tree-like network [53]). For this reason, $P$ is denoted as the infection propagator.

To compute the continuous-time limit of the infection propagator, we observe that $P$ obeys the recursive relation $P(h+1)=P(h)\left[1-\mu\Delta t+\lambda\Delta tA_{h+1}\right]$. Expressed in continuous time and dividing both sides by $\Delta t$, the relation becomes:

[TABLE]

that in the limit $\Delta t\rightarrow 0$ yields

[TABLE]

a system of $N^{2}$ coupled differential equations whose components are

[TABLE]

The lhs of Eqs. (4) is the derivative of $P$ that is well behaved if all entries are continuous functions of time. $A_{ij}(t)$ are, however often binary, so that their evolution is a sequence of discontinuous steps. To overcome this, it is possible to approximate these steps with one-parameter families of continuous functions, compute the threshold, and then perform the limit of the parameter that recovers the discontinuity. More formally, this is equivalent to interpret derivatives in the sense of tempered distributions [54].

In order to check that our limit process correctly connects the discrete-time framework to the continuous time one, let us now consider the standard Markov chain formulation of the continuous dynamics:

[TABLE]

Performing a linear stability analysis of the disease-free state (i.e., around $p_{i}(t)=0$) in the quenched mean-field approximation [17, 18], we obtain

[TABLE]

We note that this expression is formally equivalent to Eq. (5). In particular, each row of $P_{ij}$ of Eq. (5) satisfies Eq. (7). Furthermore, the initial condition $P_{ij}(0)=\delta_{ij}$ guarantees that in varying the row $i$ we consider all vectors of the space basis as initial condition. Every solution $p(T)$ of Eq. (7) can therefore be expressed as a linear combination of the rows of $P(T)$. Any fundamental matrix solution of Eq. (7) obeys Eq. (5) within the framework of Floquet theory of nonautonomous linear systems [55].

The equivalence of the two equations shows that our limit of the discrete-time propagator encodes the dynamics of the continuous process. It is important to note that the limit process leading to Eq. (4) entails a fundamental change of paradigm on the representation of the network structure and contagion process, where the linear algebraic representation suitable in discrete time turns into a differential geometrical description of the continuous-time flow. While network and spreading dynamics in discrete time are encoded in a multilayer adjacency tensor, the continuous time description proposed in Eq. (5) rests on a representation of the dynamical process in terms of a manifold whose points are adjacency matrices (or rank-2 tensor in the sense of [49]) corresponding to possible network and contagion states. The dynamics of Eq. (5) is then a curve on such manifold, indicating which adjacency matrices to visit and in which order. In practice, we recover that the contagion process on a discrete temporal network corresponding to an ordered subset of the full multilayer structure of [49] becomes in the limit $\Delta t\rightarrow 0$ a spreading on a continuous temporal network represented through a 1-dimensional ordered subset of a tensor field (formally the pull-back on the evolution curve). The two frameworks, so far considered independently and mutually exclusive, thus merge coherently through a smooth transition in this novel perspective.

We now turn to solving Eq. (4) to derive an analytic expression of the infection propagator. By defining the rescaled transmissibility $\gamma=\lambda/\mu$, we can solve Eq. (4) in terms of a series in $\mu$ [56]

[TABLE]

with $P^{(0)}=1$ and under the assumption that $\gamma$ remains finite around the epidemic threshold for varying recovery rates. The recursion relation from which we derived Eq. (4) provides the full propagator for $t=T$. Eq. (8) computed in $T$ therefore yields the infection propagator for the continuous-time adjacency matrix $A(t)$, and is defined by the sum of the following terms

[TABLE]

Equations (8) and (9) can be put in a compact form by using Dyson’s time-ordering operator $\mathcal{T}$ [57]. It is defined as $\mathcal{T}A(t_{1})A(t_{2})=A(t_{1})A(t_{2})\theta(t_{1}-t_{2})+A(t_{2})A(t_{1})\theta(t_{2}-t_{1})$, with $\theta$ being Heaviside’s step function. The expression of the propagator is thus

[TABLE]

Eq. (10) represents an explicit general solution for Eq. (4) that can be computed numerically to arbitrary precision [56]. The epidemic threshold in the continuous-time limit is then given by $\rho\left(P(T)\right)=1$.

We now discuss a special case where we can recover a closed-form solution of Eq. (10), and thus of the epidemic threshold. We consider continuously evolving temporal networks satisfying the following condition (weak commutation):

[TABLE]

i.e. the adjacency matrix at a certain time $A(t)$ commutes with the aggregated matrix up to that time. In the introduced tensor field formalism, the weak commutation condition represents a constraint on the temporal trajectory, or equivalently, an equation of motion for $A(t)$.

Eq. (11) implies that the order of factors in Eq. (9) no longer matters. Hence, we can simply remove the time-ordering operator $\mathcal{T}$ in Eq. (10), yielding

[TABLE]

where $\left\langle A\right\rangle=\int_{0}^{T}\!\operatorname{d}\!{t}A(t)/{T}$ is the adjacency matrix averaged over time. The resulting expression for the epidemic threshold for weakly commuting networks is then

[TABLE]

This closed-form solution proves to be extremely useful as a wide range of network classes satisfies the weak commutation condition of Eq. (11). An important class is constituted by annealed networks [50, 4, 51]. In the absence of dynamical correlations, the annealed regime leads to $\langle[A(x),A(y)]\rangle=0$, as the time ordering of contacts becomes irrelevant. Eq. (11) can thus be reinterpreted as $\langle[A(t),A(x)]\rangle_{x}=0$, where the average is carried out over $x\in\left[0,t\right)$. For long enough $t$, $\int_{0}^{t}\!\operatorname{d}\!{x}A(x)/t$ approximates well the expected adjacency matrix $\left\langle A\right\rangle$ of the annealed model, leading the annealed regime to satisfy Eq. (13). This result thus provides an alternative mathematical framework for the conceptual interpretation of annealed networks in terms of weak commutation. Originally introduced to describe disorder on quenched networks [58, 59], annealed networks were mathematically described in probabilistic terms, with the probability of establishing a contact depending on the degree distribution $P(k)$ and the two-node degree correlations $P(k^{\prime}|k)$ [50]. Here we show that temporal networks whose adjacency matrix $A(t)$ asymptotically commutes with the expected adjacency matrix are found to be in the annealed regime.

Eq. (13) can also be used to test the limits of the time scale separation approach, by considering a generic temporal network not satisfying the weak commutation condition. If $\mu$ is small, we can truncate the series of the infection propagator (Eq. (8)) at the first order, $P=1+\mu P^{(1)}+{\mathcal{O}}(\mu^{2})$, where $P^{(1)}(T)=T[\gamma\left\langle A\right\rangle-1]$, to recover indeed Eq. (13). The truncation thus provides a mathematical expression of the range of validity of the time-separation scheme for spreading processes on temporal networks, since temporal correlations can be disregarded when the network evolves much faster than the spreading process.

Extending the result of the annealed networks, we show that the weak commutation condition also holds for networks whose expected adjacency matrix depends on time as a scalar function (instead of being constant as in the annealed case), $\left\langle A(t)\right\rangle=c(t)\left\langle A(0)\right\rangle$. Also in this case we have $\langle[A(x),A(y)]\rangle=0$, so that the same treatment performed for annealed networks applies. Examples are provided by global trends in activation patterns, as often considered in infectious disease epidemiology to model seasonal variations of human contact patterns (e.g. due to the school calendar) [60].

When the time scale separation approach is not applicable, we find another class of weakly commuting temporal networks that are used as a paradigmatic network example for the study of contagion processes occurring on the same time scale of contacts evolution – the activity-driven model [35]. It considers heterogeneous populations where each node $i$ activates according to an activity rate $a_{i}$, drawn from a distribution $f(a)$. When active, the node establishes $m$ connections with randomly chosen nodes lasting a short time $\delta$ ($\delta\ll 1/a_{i}$). Since the dynamics lacks time correlations, the weak commutation condition holds, and the epidemic threshold can be computed from Eq. (13). In the limit of large network size, it is possible to write the average adjacency matrix as $\left\langle A\right\rangle_{ij}=\frac{m\delta}{N}(a_{i}+a_{j})+{\mathcal{O}}(\frac{1}{N^{2}})$. Through row operations we find that the matrix has $\mbox{rank}(\left\langle A\right\rangle)=2$, and thus only two non-zero eigenvalues, $\alpha,\sigma$, with $\alpha>\sigma$. We compute them through the traces of $\left\langle A\right\rangle$ ($\mbox{tr}\!\left[\left\langle A\right\rangle\right]=\alpha+\sigma$ and $\mbox{tr}[\left\langle A\right\rangle^{2}]=\alpha^{2}+\sigma^{2}$) to obtain the expression of $\rho[\left\langle A\right\rangle]$ for Eq. (13): $\rho[\left\langle A\right\rangle]=\alpha=m\delta\left(\left\langle a\right\rangle+\sqrt{\left\langle a^{2}\right\rangle}\right)$. The epidemic threshold becomes

[TABLE]

yielding the same result of Ref. [35], provided here that the transmission rate $\lambda$ is multiplied by $\delta$ to make it a probability, as in [35].

Finally, we verify that for the trivial example of static networks, with an adjacency matrix constant in time, Eq. (13) reduces immediately to the result of Refs. [17, 18].

We now validate our analytical prediction against numerical simulations on two synthetic models. The first is the activity-driven model with activation rate $a_{i}=a$, $m=1$, and average inter-activation time $\tau=1/a=1$, fixed as the time unit of the simulations. The transmission parameter is the probability upon contact $\lambda\delta$ and the model is implemented in continuous time. The second model is based on a bursty inter-activation time distribution $P(\Delta t)\sim(\epsilon+\Delta t)^{-\beta}$ [31], with $\beta=2.5$ and $\epsilon$ tuned to obtain the same average inter-activation time as before, $\tau=1$. We simulate a SIS spreading process on the two networks with four different recovery rates, $\mu\in\{10^{-3},10^{-2},10^{-1},1\}$, i.e. ranging from a value that is 3 orders of magnitude larger than the time scale $\tau$ of the networks (slow disease), to a value equal to $\tau$ (fast disease). We compute the average simulated endemic prevalence for specific values of $\lambda,\mu$ using the quasi-stationary method [61] and compare the threshold computed with Eq. (13) with the simulated critical transition from extinction to endemic state. As expected, we find Eq. (13) to hold for the activity-driven model at all time scales of the epidemic process (Fig. 2), as the network lacks temporal correlations. The agreement with the transition observed in the bursty model however is recovered only for slow diseases, as at those time scales the network is found in the annealed regime. When network and disease time scales become comparable, the weakly commuting approximation of Eq. (13) no longer holds, as burstiness results in dynamical correlations in the network evolution [31].

Our theory offers a novel mathematical framework that rigorously connects discrete-time and continuous-time critical behaviors of spreading processes on temporal networks. It uncovers a coherent transition from an adjacency tensor to a tensor field resulting from a limit performed on the structural representation of the network and contagion process. We derive an analytic expression of the infection propagator in the general case that assumes a closed-form solution in the introduced class of weakly commuting networks. This allows us to provide a rigorous mathematical interpretation of annealed networks, encompassing the different definitions historically introduced in the literature. This work also provides the basis for important theoretical extensions, assessing, for example, the impact of bursty activation patterns or of the adaptive dynamics in response to the circulating epidemic. Finally, our approach offers a tool for applicative studies on the estimation of the vulnerability of temporal networks to contagion processes in many real-world scenarios, for which the discrete-time assumption would be inadequate.

Acknowledgments We thank Luca Ferreri and Mason Porter for fruitful discussions. This work is partially sponsored by the EC-Health contract no. 278433 (PREDEMICS) and the ANR contract no. ANR-12-MONU-0018 (HARMSFLU) to VC; the EC-ANIHWA contract no. ANR-13-ANWA-0007-03 (LIVEepi) to EV, CP, VC.