Introduction

In recent years epidemiological modelling, along with many other fields, has seen renewed activity thanks to the emergence of network science (Newman 2018; Barabasi 2016; Zhan et al. 2020; Masuda and Holme 2017). Approaching these models from the view of complex coupled systems has shed new light onto spreading processes where the early black-box ordinary differential equation (ODE) models from Kermack and McKendrick had its limitations (Keeling and Eames 2005; Humphries et al. 2020). These ODE models assume homogeneous mixing of the entire population, which may be an appropriate approximation for small communities. However, when attempting to model the spread of disease at a national or international level, they fail to capture how heterogeneities in both travel patterns and population distributions contribute to and affect the spread of disease. Epidemiological models on complex networks aim to solve this problem by moving away from averaged dynamics of populations and mean-field descriptions. Instead, the focus is on interactions between individuals or metapopulations, where the spreading process is driven by contacts in the network (Yang et al. 2003; Sharkey et al. 2015; Colizza et al. 2007).

There have been many improvements made in regards to network models, e.g., generalised multi-layer network structures or more specifically temporal networks that allow for the network structure to change with time (Masuda and Holme 2017; Iannelli et al. 2017; Koher et al. 2016; Lentz et al. 2016). Temporal networks are a natural way of representing contacts and lead to an insightful interplay between the disease dynamics and the evolving network topology (Karrer and Newman 2010; Shrestha et al. 2015; Lentz et al. 2013). With the ever growing availability of mobility and contact data it has become easier to provide accurate and high-resolution data to inform network models. The results can be extremely useful tools for public-health bodies and other stakeholders (Gethmann et al. 2019; Tratalos et al. 2020; Génois and Barrat 2018).

In previous works, a widely used epidemiological concept is the individual-based model (Newman 2018; Valdano et al. 2015; Sharkey 2011). It assumes statistical independence in the state of each vertex. A major problem associated with such a model is that it suffers quite badly from an echo chamber effect due to the fact that there is no memory of past interactions due to statistical independence. There have been efforts to ameliorate this problem by introducing memory at the level of each vertex’s direct neighbours. These models referred to as contact-based (Koher et al. 2019) or pair-based (Frasca and Sharkey 2016) and have been shown to significantly reduce the echo chamber effect, depending on the underlying network structure. These two models differ in their initial approach. The contact-based model takes an edge-based perspective, which extends the message passing approach (Karrer and Newman 2010; Kirkwood 1935), and all dynamic equations are formulated in terms of edges. By contrast, the pair-based model keeps the vertex-based approach of the individual-based model and dynamic equations are in terms of vertices.

In this paper, we extend the Pair-Based (PB) model to a temporal setting giving a Temporal Pair-Based (TPB) Model. We show how it can be drastically reduced and simplified under a certain dynamical assumption (Sharkey and Wilkinson 2015). We deal specifically with susceptible–infected–recovered (SIR) models. Once the TPB model is written in concise form, it is then possible to show that the contact-based model is equivalent to a linearised version of the TPB model. We then establish the conditions for an epidemic to occur according to the TPB model, also known as the epidemic threshold. We investigate how the TPB model performs on a number of synthetic and empirical networks and investigate what kind of network topologies work best with the TPB model.

The remainder of this paper is structured as follows: In "SIR network model" section, we summarize the theoretical framework. Then, we address the calculation of an epidemic threshold in "Epidemic threshold" section, before presenting the main results in "Results" section. Finally, we offer some conclusions in "Conclusions" section.

SIR network model

Let us consider a temporal network \({\mathcal {G}}=\{G_1,\dots ,G_T\}\) to be a series of networks \(G_t = (V, E_t)\), which all share the same vertex set V but differ in their edge set \(E_t\). The adjacency matrix for the network at time t will be denoted by \(A^{[t]}\), and \(A^{[t]}_{ij}=1\) implies an edge between vertices i and j at time t, whereas \(A^{[t]}_{ij}=0\) implies no edge. 

Reduced master equations

Let \(\Omega\) be the set of compartments in the model, that is, in the SIR model: \(\Omega = \{S,I,R\}\). Let \({\mathbf {X}}^n=\left( X_1^n, X_2^n,\dots ,X_N^n \right) ^T\) be the vector whose ith element refers to the state of the ith vertex in the network at time step n, thus \({\mathbf {X}}^n\in \Omega ^{N}\) where \(|V|=N\). The evolution of the disease is then described by the master equations (Gardiner and Gardiner 2009),

$$\begin{aligned} P\left( {\mathbf {X}}^{n+1}\right) = \sum _{{\mathbf {X}}^{n}\in \Omega ^{N}}P\left( {\mathbf {X}}^{n+1}|{\mathbf {X}}^{n}\right) P\left( {\mathbf {X}}^{n}\right) . \end{aligned}$$
(1)

In other words, we assume that the infection process is Markovian. \(P({\mathbf {X}}^{n+1})\) is the probability of the network being in the particular configuration of states given by \({\mathbf {X}}^{n+1}\) and \(P({\mathbf {X}}^{n+1} | {\mathbf {X}}^{n})\) is probability of the network moving from the configuration of states \({\mathbf {X}}^{n}\) to \({\mathbf {X}}^{n+1}\) between their respective time steps. These equations describe the entire process on the network. However, in order to progress the system forward one step in time, the probabilities of all combinations of state vectors must be found. This usually is not feasible for network processes with potentially billions of vertices as for the SIR process the total combination of states is given by \(3^N\).

Instead, it is possible to describe the evolution of the disease using a system of Reduced Master Equations (RME) (Sharkey 2011), which describes the evolution of subsystems within the network, such as individual vertices, removing the need to obtain every possible combination of states. An important note is that these RMEs are in fact not themselves true master equations as they are not necessarily linear due to the fact that the transition rates of the subsystems are non-linear combinations of the transitions rates of the original system. However, we shall continue to use the term RME introduced by the author of Sharkey (2011).

For notational convenience we use the following notation for the joint marginal probabilities,

$$\begin{aligned} P\left( X^{n}_{i_1},X^{n}_{i_2},\cdots ,X^{n}_{i_m}\right)&= \langle X^{n}_{i_1}X^{n}_{i_2}\cdots X^{n}_{i_m}\rangle . \end{aligned}$$
(2)

When we wish to specify a particular realisation of \(X_i^n\), we denote it by \(S_i^n\), \(I_i^n\) or \(R_i^n\) to imply \(X_i^n=S\), \(X_i^n=I\) or \(X_i^n=R\), respectively. Employing this new notation we start with the RME, which describes the evolution of individual vertices,

$$\begin{aligned} \langle X^{n+1}_{i} \rangle = \sum _{X^n_i\in \Omega }\langle X^{n+1}_i|X_i^n \rangle \langle X_i^n \rangle . \end{aligned}$$
(3)

For SIR dynamics, the evolution of each vertex in each compartment is given as the following,

$$\begin{aligned} \langle S^{n+1}_i\rangle&= \langle S^{n}_i\rangle -\langle I_i^{n+1}|S^{n}_i \rangle \langle S^{n}_i\rangle \end{aligned}$$
(4a)
$$\begin{aligned} \langle I^{n+1}_i\rangle&= \langle I^{n}_i\rangle -\langle R_i^{n+1}|I^{n}_i \rangle \langle I^{n}_i\rangle + \langle I_i^{n+1}|S^{n}_i \rangle \langle S_i^n\rangle \end{aligned}$$
(4b)
$$\begin{aligned} \langle R^{n+1}_i\rangle&=\langle R^{n}_i\rangle + \langle R_i^{n+1}|I^{n}_i \rangle \langle I_i^n\rangle , \end{aligned}$$
(4c)

where \(\langle I_i^{n+1}|S^{n}_i \rangle\) reads the probability vertex i is infected at time \(n+1\) given it was susceptible at time n and similarly for \(\langle R_i^{n+1}|I^{n}_i \rangle\). Note that \(\langle R_i^n\rangle\) can be recovered using the conservation of the probabilities \(\langle S_i^n\rangle + \langle I_i^n\rangle + \langle R_i^n\rangle = 1\). In order to compute the transition rates we define the following quantities, the probability of infection on contact \(\beta\), the rate of recovery \(\mu\). \(A^{[n]}\) is the temporal adjacency matrix of the network on which the process is occurring. Following directly from Frasca and Sharkey (2016), the transition rates of moving from S to I, and I to R are given by

$$\begin{aligned} \langle I_i^{n+1}|S^{n}_i \rangle&= \frac{1}{\langle S_i^n\rangle }\left[ \beta \sum _{j_1\in V}A_{ij_1}^{[n]} \langle S_i^n I_{j_1}^n\rangle \right. - \beta ^2\sum _{j_1, j_2\in V} A_{ij_1}^{[n]}A_{ij_2}^{[n]}\langle S_i^n I_{j_1}^n I_{j_2}^n\rangle \nonumber \\&\quad +\dots \left. -\, (-\beta )^{N-1}{\sum _{j_1,\dots ,j_{N-1}\in V}} A_{ij_1}^{[n]}\dots A_{ij_{N -1}}^{[n]}\langle S_i^n I_{j_1}^n \dots I_{j_{N-1}}^n\rangle \right] \end{aligned}$$
(5a)
$$\begin{aligned} \langle R_i^{n+1}|I^{n}_i \rangle&= \mu . \end{aligned}$$
(5b)

For the expression within the square brackets of Eq. (5a), the first term is the probability that vertex i is infected by some other vertex in the network and double-counts events. Each subsequent term accounts for double-counting and over-correcting in the previous. These equations describe the probabilistic SIR process on temporal networks. Note that the system of equations is not closed as it lacks a description for their joint probabilities. There are a number of ways in which this problem can be tackled, usually by making a number of numerical or dynamical approximations (Yang et al. 2003; Valdano et al. 2015; Koher et al. 2019; Gómez et al. 2010). In the next sections we attempt to improve on and unify many existing approaches showing how they are derived from the system of RMEs given by Eqs. (4) and (5).

Temporal individual-based model

One of the most commonly used epidemiological models on networks is the individual-based (IB) model, which has been extended to the temporal setting in Valdano et al. (2015). We refer to this extension as the temporal individual-based (TIB) model. The key idea is the assumption of statistical independence of vertices or the mean field approximation, i.e., the factorisation \(\langle X_{i_1}^nX_{i_2}^n\dots X_{i_M}^n\rangle =\langle X_{i_1}^n\rangle \langle X_{i_2}^n\rangle \dots \langle X_{i_M}^n\rangle\). By assuming this independence of vertices, Eq. (5a) simplifies to

$$\begin{aligned} \langle I_i^{n+1}|S^{n}_i \rangle&= \beta \sum _{j_1\in V}A_{ij_1}^{[n]} \langle I_{j_1}^n\rangle \nonumber \\&\quad -\, \beta ^2\sum _{j_1, j_2\in V} A_{ij_1}^{[n]}A_{ij_2}^{[n]}\langle I_{j_1}^n\rangle \langle I_{j_2}^n\rangle +\dots \nonumber \\&\quad -\,(-\beta )^{N-1}{\sum _{j_1,\dots ,j_{N-1}\in V}} A_{ij_1}^{[n]}\dots A_{ij_{N -1}}^{[n]}\langle I_{j_1}^n\rangle \dots \langle I_{j_{N-1}}^n\rangle , \end{aligned}$$
(6a)

which upon factorising can be concisely written as

$$\begin{aligned} \langle I_i^{n+1}|S^{n}_i \rangle = 1 - \prod _{k\in V}\left( 1-\beta A_{ki}^{[n]}\langle I_k^n\rangle \right) . \end{aligned}$$
(7)

Upon substituting the transition rates \(\langle I_i^{n+1}|S^{n}_i \rangle\) and \(\langle R_i^{n+1}|I^{n}_i \rangle\) under the assumption of statistical independence, the full TIB model is written as

$$\begin{aligned} \langle S^{n+1}_i\rangle&= \langle S^{n}_i\rangle \prod _{k\in V}\left( 1-\beta A_{ki}^{[n]}\langle I_k^n\rangle \right) \end{aligned}$$
(8a)
$$\begin{aligned} \langle I^{n+1}_i\rangle&= \langle I^{n}_i\rangle (1-\mu )+\langle S^{n}_i\rangle \left( 1- \prod _{k\in V}\left( 1-\beta A_{ki}^{[n]}\langle I_k^n\rangle \right) \right) . \end{aligned}$$
(8b)

This model closes Eq. (5a) at the level of vertices, thus ignoring all correlations with other vertices at previous times. However, ignoring all past correlations causes the model to suffer quite badly from an echo chamber effect (Shrestha et al. 2015). This echo chamber has the effect of vertices artificially amplifying each other’s probability of being infected \(\langle I_i^n \rangle\) at each new time step, as the marginal probability of each vertex is highly correlated with the rest of the network and the factorisation of Eq. (5a) means each vertex forgets its past interactions. As demonstrated in Shrestha et al. (2015), it is possible to show that in the absence of a recovered compartment, a static network of two linked vertices for non-zero initial conditions has probabilities of being infected which converge according to \(\lim _{n\rightarrow \infty }\langle I_0^n\rangle =\lim _{n\rightarrow \infty }\langle I_1^n\rangle =1\) for the TIB model.

Temporal pair-based model

In contrast to the TIB model, instead of assuming independence of vertices we can approximate the marginal probabilities in terms of combinations of lower order marginals using some form of moment closure (Frasca and Sharkey 2016; Sharkey and Wilkinson 2015). Here, we make an equivalent assumption to that of the message passing approaches (Karrer and Newman 2010; Shrestha et al. 2015). We assume the network contains no time-respecting non-backtracking cycles. In other words, starting at some initial vertex i that leaves via vertex j, there is no way to find a time-respecting path returning to this vertex that does not return via j. This is equivalent to a tree network when the temporal network is viewed in its static embedding of the supra-adjacency representation (Bianconi 2018). This allows us to write all higher order moments in Eq. (5a) as a combination of pairs \(\langle S_i^n I_k^n\rangle\). To show why this is possible, consider the three vertices ijk connected by two edges through i. If conditional independence of these vertices is assumed given we have the state of i, then one can make the following assumption,

$$\begin{aligned} \langle X_i^nX_j^nX_k^n\rangle = \langle X_j^nX_k^n|X_i^n\rangle \langle X_i^n\rangle = \frac{\langle X_i^nX_j^n\rangle \langle X_i^nX_k^n\rangle }{\langle X_i^n\rangle }. \end{aligned}$$
(9)

This has the effect of assuming the network is tree-like in structure as it implies any interaction between vertices j and k must occur through vertex i. Thus, the process is exact on networks that contain no time-respecting non-backtracking cycles and otherwise provides an improved approximation of varying degree, which depends on the true network structure. The result obtained in Eq. (9) is often referred to as the Kirkwood closure (Kirkwood 1935). Under the assumption that the network is tree like, the following simplification is obtained for Eq. (5a),

$$\begin{aligned} \langle I_i^{n+1}|S^{n}_i \rangle =1 - \prod _{k\in V}\left( 1 - \beta A^{[n]}_{ki}\frac{ \langle S_i^nI_k^n\rangle }{\langle S_i^n\rangle }\right) . \end{aligned}$$
(10)

However, we run into the problem that we have no description for pairs of vertices. Thus, we derive expressions for their evolution from the RMEs for pairs of vertices which is given by,

$$\begin{aligned} \langle X_i^{n+1}X_j^{n+1} \rangle = \sum _{(X_{i}^n,X_{j}^n)\in \Omega ^2}\langle X_i^{n+1}X_j^{n+1}|X_i^nX_j^n \rangle \langle X_i^nX_j^n \rangle , \end{aligned}$$
(11)

For \(\langle S_i^{n+1}I_j^{n+1} \rangle\), we obtain

$$\begin{aligned} \langle S_i^{n+1}I_j^{n+1} \rangle&= \langle S_i^{n}I_j^{n}\rangle + \langle S_i^{n+1}I_j^{n+1}|S_i^nS_j^n \rangle \langle S_i^nS_j^n \rangle \\&\quad -\, \langle I_i^{n+1}I_j^{n+1}|S_i^nI_j^n \rangle \langle S_i^nI_j^n \rangle \\&\quad -\, \langle S_i^{n+1}R_j^{n+1}|S_i^nI_j^n \rangle \langle S_i^nI_j^n \rangle \\&\quad -\, \langle I_i^{n+1}R_j^{n+1}|S_i^nI_j^n \rangle \langle S_i^nI_j^n \rangle . \end{aligned}$$
(12)

Note that the RME for \(\langle S_i^{n}S_j^{n} \rangle\) is also required, which we find to be the following

$$\begin{aligned} \langle S_i^{n+1}S_j^{n+1}\rangle&= \langle S_i^nS_j^n\rangle -\langle S_i^{n+1}I_j^{n+1}|S_i^nS_j^n \rangle \langle S_i^nS_j^n\rangle \\&\quad -\,\langle I_i^{n+1}I_j^{n+1}|S_i^nS_j^n \rangle \langle S_i^nS_j^n\rangle \\&\quad -\,\langle I_i^{n+1}I_j^{n+1}|S_i^nS_j^n \rangle \langle S_i^nS_j^n\rangle . \end{aligned}$$
(13)

Since only the probabilities \(\langle S_i^nI_j^n \rangle\) and \(\langle S_i^nS_j^n \rangle\) are needed in order to describe the RMEs in Eq. (10), we consider those two combinations of states. From Frasca and Sharkey (2016), we obtain the exact transition rates for pairs of vertices and find that we can factorise the pair-wise transition rates similar to Eq. (5a). Here, we give the expression for \(\langle S_i^{n+1}I_j^{n+1}|S_i^nS_j^n \rangle\) only while the rest of the pair-wise transition rates are given in “Appendix”:

$$\begin{aligned}{}&\langle S_i^{n+1}I_j^{n+1}|S_i^nS_j^n \rangle \nonumber \\&\quad = \frac{1}{\langle S_i^nS_j^n\rangle }\left[ \beta \sum _{k_1\in V}A_{ik_1}^{[n]} \langle S_i^n S_j^nI_{k_1}^n\rangle \right.- \beta ^2\sum _{k_1, k_2\in V} A_{ik_1}^{[n]}A_{ik_2}^{[n]}\langle S_i^nS_j^n I_{k_1}^n I_{k_2}^n\rangle \nonumber \\&\qquad +\dots \nonumber \left. - (-\beta )^{N-2}{\sum _{k_1,\dots ,k_{N-2}\in V}} A_{ik_1}^{[n]} \dots A_{ik_{N -2}}^{[n]}\langle S_i^n S_{j}^n \dots I_{k_{N-2}}^n\rangle \right] \nonumber \\&\qquad \times \left[ 1-\beta \sum _{k_1\in V}A_{ik_1}^{[n]} \langle S_i^n S_j^nI_{k_1}^n\rangle \right. + \beta ^2\sum _{k_1, k_2\in V} A_{ik_1}^{[n]}A_{ik_2}^{[n]}\langle S_i^nS_j^n I_{k_1}^n I_{k_2}^n\rangle \nonumber \\&\qquad -\dots \left. + (-\beta )^{N-2}{\sum _{k_1,\dots ,k_{N-2}\in V}} A_{ik_1}^{[n]}\dots A_{ik_{N -2}}^{[n]}\langle S_i^n S_{j}^n \dots I_{k_{N-2}}^n\rangle \right] . \end{aligned}$$
(14)

In the above equation, the term in the first pair of square brackets corresponds to the probability that vertex i does not become infected and the term in the second pair of square brackets corresponds to the probability that vertex j becomes infected. Upon applying our moment closure technique Eq. (14) may be written as

$$\begin{aligned} \langle S_i^{n+1}&I_j^{n+1}|S_i^nS_j^n \rangle \nonumber \\&=\prod _{\begin{array}{c} k\in V\\ k\ne j \end{array}} \left( 1- \beta A^{[n]}_{ki}\frac{\langle S_i^nI_k^n \rangle }{\langle S_i^n\rangle }\right) \left[ 1 - \prod _{\begin{array}{c} k\in V\\ k\ne i \end{array} }\left( 1- \beta A^{[n]}_{kj}\frac{\langle S_j^nI_k^n \rangle }{\langle S_j^n\rangle }\right) \right] . \end{aligned}$$
(15)

By introducing the following functions, the RMEs for pairs as well as the individual vertices can be written more concisely. The probability that vertex i does not become infected at time step \(n+1\), given that i is not infected at time step n is denoted by

$$\begin{aligned} \Psi _i^n=\prod _{k\in V}\left( 1 - \beta A^{[n]}_{ki}\frac{ \langle S_i^nI_k^n\rangle }{\langle S_i^n\rangle }\right) . \end{aligned}$$
(16)

Similarly, the probability that vertex i does not become infected at time step \(n+1\), given that i is not infected at time step n while excluding any interaction with j, is given by

$$\begin{aligned} \Phi _{ij}^n=\prod _{\begin{array}{c} k\in V\\ k\ne j \end{array}}\left( 1 - \beta A^{[n]}_{ki}\frac{\langle S_i^nI_k^n\rangle }{\langle S_i^n\rangle }\right) . \end{aligned}$$
(17)

Then, the evolution of the state of every vertex in the network is determined by the following closed set of equations,

$$\begin{aligned} \langle S_i^{n+1} \rangle&= \Psi _i^n\langle S_i^n \rangle \end{aligned}$$
(18a)
$$\begin{aligned} \langle I_i^{n+1}\rangle&= (1- \mu )\langle I_i^n\rangle + [1 - \Psi _i^n]\langle S_i^n\rangle \end{aligned}$$
(18b)
$$\begin{aligned} \langle S_i^{n+1}I_j^{n+1}\rangle&= (1-\mu ) (1-\beta A_{ji})\Phi _{ij}^n\langle S_i^nI_j^n\rangle \nonumber \\&\quad +\, \Phi _{ij}^n\left( 1-\Phi _{ji}^n\right) \langle S_i^nS_j^n\rangle \end{aligned}$$
(18c)
$$\begin{aligned} \langle S_i^{n+1}S_j^{n+1}\rangle&=\Phi _{ij}^n\Phi _{ji}^n \langle S_i^nS_j^n\rangle . \end{aligned}$$
(18d)

This approximation allows a large increase in accuracy compared to TIB model while only adding two equations to the final model. All past dynamic correlations are now tracked by the model and so the echo chamber effect is eliminated, but only with direct neighbours, that is, vertices which share an edge. A major benefit of this particular TPB model over other existing iterations (Koher et al. 2019; Gómez et al. 2010) is that this model can be implemented as an element-wise sparse matrix multiplication rather than having to iterate through all edges for every time step, making it extremely computationally efficient and fast on even large networks. It also benefits from a low conceptual cost by not deviating from a vertex-based perspective, like the contact-based models, which move to the perspective of edges and thus define the model in terms of the line-graphs and non-backtracking matrices (Koher et al. 2019).

Similar to Shrestha et al. (2015), we can compare pair-based models to the TIB model using the two vertex example. In that illustrative configuration, we consider two vertices connected by an undirected static edge and give the two vertices some initial non-zero probability \(\langle I_0^0\rangle = \langle I_1^0\rangle = z\) of being infected. We then run the TIB and TPB models for some given parameters \(\beta\) and \(\mu\) and compare it to the ground truth, which is the average of a number of Monte-Carlo(MC) realisations.

Fig. 1
figure 1

Running 25 time steps of the TIB (blue dashed) and TPB (green solid) SIR model as well as the average over 10\(^5\) Monte-Carlo (MC) simulations (red dash-dotted) for the two vertex example. Parameters: \(\beta = 0.4\), \(\mu = 0.2\), and \(\langle I^0_0 \rangle = \langle I^0_1 \rangle =0.2\)

From Fig. 1, it becomes apparent how the TIB model fails to capture the true SIR process on the network due to the previously discussed echo chamber induced by assuming statistical independence of vertices. It becomes clear that the TPB model accurately describes the underlying SIR process for this simple example as each vertex is able to recover the dynamic correlations of past interactions with direct neighbours.

Equivalence between the contact-based and pair-based models

In the contact-based model (Koher et al. 2019), the central component is \(\theta _{ij}^n\), which is the probability that node j has not passed infection to node i up to time step n. From \(\theta _{ij}^n\), the quantity \(\langle S^{n}_i\rangle\) may be computed as

$$\begin{aligned} \langle S_i^{n+1}\rangle = \langle S_i^0\rangle \prod _{j\in V}\theta _{ij}^{n+1}. \end{aligned}$$
(19)

This equation is the basis for the contact-based model and allows us to easily compare with the pair-based model as it describes the same quantity as our Eq. (18a). The authors also assume that the evolution of \(\theta _{ij}^n\) satisfies the following relation

$$\begin{aligned} \theta _{ij}^{n+1}&= \theta _{ij}^n - \beta A^{[n]}_{ji}\frac{\langle S_i^nI_j^n\rangle }{\langle S_i^n\rangle }\nonumber \\ \theta _{ij}^{0}&=1. \end{aligned}$$
(20)

In the pair-based model, the evolution of the susceptible probability, given by Eq. (18a), can similarly be rewritten in terms of its initial conditions,

$$\begin{aligned} \langle S_i^{n+1}\rangle&= \Psi _i^n\langle S_i^n\rangle \end{aligned}$$
(21a)
$$\begin{aligned}{}&= \langle S_i^0\rangle \prod ^{n}_{m=0}\Psi _i^m \end{aligned}$$
(21b)
$$\begin{aligned}{}&= \langle S_i^0\rangle \prod _{j\in V}\prod ^{n}_{m=0} \left( 1 - \beta A^{[m]}_{ji}\frac{\langle S_i^mI_j^m\rangle }{\langle S_i^m\rangle }\right) . \end{aligned}$$
(21c)

From equating (19) and (21) it is clear that if the models are exactly equivalent then \(\theta _{ij}\) is defined by

$$\begin{aligned} \theta _{ij}^{n+1} = \prod ^{n}_{m=0}\left( 1 -\beta A^{[m]}_{ji}\frac{\langle S_i^mI_j^m\rangle }{\langle S_i^m\rangle }\right) . \end{aligned}$$
(22)

However, this contradicts the assumption made by Eq. (20). Thus the pair-based and contact-based models are only equivalent if the following linearisation is assumed:

$$\begin{aligned} \prod ^{n}_{m=0}\left( 1 - \beta A^{[m]}_{ji}\frac{\langle S_i^mI_j^m\rangle }{\langle S_i^m\rangle }\right) \approx 1 - \sum ^n_{m = 0}\beta A^{[m]}_{ji}\frac{\langle S_i^mI_j^m\rangle }{\langle S_i^m\rangle }, \end{aligned}$$
(23)

which then implies (22) can be written as

$$\begin{aligned} \theta _{ij}^{n+1}&= 1 - \sum ^{n-1}_{m = 0}\beta A^{[m]}_{ji}\frac{\langle S_i^mI_j^m\rangle }{\langle S_i^m\rangle } - \beta A^{[n]}_{ji}\frac{\langle S_i^nI_j^n\rangle }{\langle S_i^n\rangle } \end{aligned}$$
(24a)
$$\begin{aligned}{}&= \theta _{ij}^n - \beta A^{[n]}_{ji}\frac{\langle S_i^nI_j^n\rangle }{\langle S_i^n\rangle }. \end{aligned}$$
(24b)

This shows that the contact-based model is a linearised version of the pair-based model.

Epidemic threshold

One of the most important metrics used in epidemiology is the epidemic threshold, which determines the critical values of the model parameters at which a transition in qualitative behaviour occurs and the disease-free equilibrium (DFE), which we define as

$$\begin{aligned} \langle S_i^{n} \rangle = S_i^*,\quad \langle I_i^n\rangle = 0,\quad \langle R_i^n \rangle = R_i^*\quad \forall i, \end{aligned}$$
(25)

where \(\langle S_i^*\rangle + \langle R_i^*\rangle = 1\). From this definition it is clear that there exists a whole class of DFEs, which must considered. At this critical point, the DFE becomes unstable and on average an epidemic occurs. Calculating the epidemic threshold is a bit more difficult with the SIR model compared with the SIS due to the fact that the flow of probability is in only one direction between compartments \(S\rightarrow I \rightarrow R\). Therefore, the class of DFE solutions are always asymptotically stable. Thus, we will look at classifying the initial stability of the SIR model as we perturb it from the state,

$$\begin{aligned} \langle S_i^{n} \rangle = 1,\quad \langle I_i^n\rangle = 0,\quad \langle R_i^n \rangle = 0\quad \forall i, \end{aligned}$$
(26)

which we shall define as the pre-disease equilibrium (PDE). If it is unstable that means the disease has a chance to take hold and will spread through the network causing an epidemic before dying out. We now look at small perturbations from the PDE, if they vanish then the disease will die out and will not have a chance to propagate through the network. We shall define an epidemic in the SIR model as instability of the PDE under such perturbations. First, we look to linearise the difference equation for \(\langle I_i^{n+1}\rangle\) near the PDE, this translates to linearising the non-linear function \(\Psi _i^n\). Under the assumption \(\langle I_i^n\rangle = \epsilon _i\) for every vertex i such that \(0<\epsilon _i\ll 1,\) we find

$$\begin{aligned} 0\le \frac{\langle S_i^nI_j^n\rangle }{\langle S_i^n\rangle } \le \epsilon _j \end{aligned}$$
(27)

by the fact that for the joint probability \(\langle S_i^nI_j^n\rangle \le \min \{\langle S_i^n\rangle , \langle I_j^n\rangle \}\). Thus for \(\epsilon _i\ll 1\) we may assume \(\frac{\langle S_i^nI_j^n\rangle }{\langle S_i^n\rangle }\approx \epsilon _j\). Upon substituting this into \(\Psi _i^n\) we find that

$$\begin{aligned} \Psi _i^n \approx&\prod _{k\in V}\left( 1 - \beta A^{[n]}_{ki}\epsilon _k\right) \nonumber \\ \approx&\sum _{k\in V}\beta A^{[n]}_{ki}\epsilon _k \end{aligned}$$
(28)

We can then use this to linearise \(\langle I_j^{n+1}\rangle\) from (18). While \(\epsilon _i\ll 1\) holds, so does the approximation,

$$\begin{aligned} \langle I_i^{n+1}\rangle&\approx \langle I_i^{n}\rangle (1-\mu ) + \sum _{k\in V} \beta A^{[n]}_{ik}\langle I_k^{n}\rangle . \end{aligned}$$
(29)

This linearisation eliminates \(\langle S_i^n\rangle\) from the equation. Interestingly, this is exactly the form of the SIS model in the TIB framework for which the epidemic threshold is easily found (Valdano et al. 2015). Therefore, we find that the SIS and SIR models share the same epidemic threshold condition. We introduce the matrix \(M^{[n]}\), called the infection propagator, which is a linear map that describes the evolution of the SIR model close to the DFE:

$$\begin{aligned} M^{[n]}_{ij} = \beta A^{[n]}_{ij} + \delta _{ij}(1-\mu ). \end{aligned}$$
(30)

Following Valdano et al. (2015), we find that the condition required for an epidemic to occur is given by

$$\begin{aligned} \rho \left( \prod _{m=0}^{n}{M}^{[m]} \right) > 1. \end{aligned}$$
(31)

For the values of \(\beta\) and \(\mu\) which the above is satisfied, implies that when a disease is introduced into the network the PDE is unstable for a period of time. What this means is that in the equivalent SIS model with the same parameters, the proportion of infected vertices never settles on the PDE.

Results

In this section, we compare the accuracy of the TIB model and the TPB model against the ground truth MC average, that is, direct stochastic simulations. In short, we show how the TPB model can offer a massive increase in accuracy and also discuss when it fails to accurately capture the true dynamics of the stochastic SIR process. Furthermore, we validate the analytical epidemic threshold.

Tree network

The assumption in the TPB model is conditional independence between vertices with a neighbour in common given the common neighbours state. This is equivalent to assuming the network contains no time-respecting non-backtracking cycles (Hashimoto 1989; Bordenave et al. 2015). To illustrate this reasoning, we consider a static tree network, that is, tree networks contain no cycles of length 3 or greater, made up of 100 vertices. All vertices start from some initial non-zero probability \(\langle I_i^0\rangle = z\) of being infected. We then run the TIB model and the TPB model for some given parameters \(\beta\) and \(\mu\) and compare it to the ground truth, which is the average of a number of MC realisations.

Fig. 2
figure 2

a A random tree network made up of 100 vertices. b Time series of the TIB (blue dashed) and TPB (green solid) SIR model as well as 10\(^5\) Monte-Carlo (MC) simulations (red dash-dotted) for the tree network shown in a. Parameters: \(\beta = 0.4\), \(\mu = 0.02\), and \(\langle I^0_0 \rangle = \langle I^0_1 \rangle =0.03\)

Figure 2b shows how the TIB model fails to capture the true SIR process on the network due to the previously discussed echo chamber induced by assuming statistical independence of vertices. It becomes clear that the TPB model accurately describes the underlying SIR process for this simple example as each vertex is able to recover the dynamic correlations of past interactions with direct neighbours. As we will see from the next section, temporal networks that are well approximated by tree networks are also well approximated by the TPB model.

Empirical networks

In the following section, we consider two empirical temporal networks that all vary in both structure and temporal activity. For each of the empirical networks we wish to test our our findings that the TPB model offers an increase in accuracy over the TIB model. We observe a change in behaviour as the model parameters cross the epidemic threshold. We run the TIB and TPB SIR models for all our networks for different values of \(\beta\) and \(\mu\) and then compare them to the average of a sufficiently large number of MC simulations. This allows us to quantify how well the different models approximate the dynamics of the true SIR process. At each time step, the average prevalence of states within the network are collected and denoted as \(\langle S^n_{\text {avg}}\rangle\), \(\langle I^n_{\text {avg}}\rangle\) and \(\langle R^n_{\text {avg}}\rangle\) with the cumulative prevalence of infection being taken as \(\langle I^n_{\text {avg}}\rangle + \langle R^n_{\text {avg}}\rangle\). Then, we validate our analytical finding for the epidemic threshold of the TPB SIR process. For this purpose, we fix a value for \(\mu\) and then, for increasing values of \(\beta\), perform a number of MC simulations for long times in order to get a distribution of the final out break size, which is given by \(\lim _{n\rightarrow \infty }\langle I^n_{\text {avg}}\rangle + \langle R^n_{\text {avg}}\rangle\). In the long-term dynamics of the SIR process, \(\lim _{n\rightarrow \infty }\langle I^n_{\text {avg}}\rangle + \langle R^n_{\text {avg}}\rangle\) will usually exceed the observation time of the network. Therefore, periodicity of the networks is assumed in a similar way to Valdano et al. (2015) when computing the final outbreak sizes.

Table 1 List of empirical networks

Irish Cattle Trade

The Cattle Trades network consists of all trades between herds within the Republic of Ireland during the year 2017 with a temporal resolution of one day (Tratalos et al. 2020) (cf. Table 1). Due to the nature of the trade data, interactions are directional. Thus, this data set is modelled by a directed network, where each vertex represents a herd and each edge represents a trade weighted by the number of animals traded. The aggregated degree distribution of the network as shown in Fig. 3a indicates a scale-free behaviour often seen in empirical networks. The network appears to be quite sparse as is evident from Fig. 3b, with an average of only 347 edges per day while having an aggregated 1,041,054 edges over the entire year. The data also displays a strong bi-modal seasonal trend with there being two distinct peaks while there tends to be very little trades occurring on Sundays when the data points lie near zero. Although we ignore external drivers of the disease, this model still offers insight into how susceptible to epidemics the network is, as trade is the main vector of non-local transmission. There are a number of infectious diseases that affect cattle, such as Foot and Mouth Disease and Bovine Tuberculosis, the latter of which is still a major problem in Ireland, thus effective models for the spread of infectious diseases among herds are incredibly useful tools. In the present study, we focus on the SIR dynamics, but the TPB model framework can easily be extended to other models,

Fig. 3
figure 3

a In- and out-degree of the network aggregated over the entire year worth of data. b Time series of number of active edges per day in the network. c Trades at the level of counties aggregated over the observation period. The number of trades is indicated by the edge width and colour

Fig. 4
figure 4

TIB (blue dashed) and TPB (green solid) SIR models on the Irish cattle trade network together with the average of \(10^3\) Monte-Carlo simulations (red dash-dotted). a, c The time series for the prevalence of infection in the network and b, d the cumulative prevalence of infection within the network. The probability of infection is set to \(\beta = 0.3\) in a, b and \(\beta = 0.5\) in c, d. The initial chance of infection is 0.03. Other parameters: \(\mu =0.005\)

From Fig. 4 we can compare the performance of the TIB and TPB models on the cattle trades network. The figure shows a year worth of simulations of both models plotted against the average of \(10^3\) MC realisations for the same choice of parameters. As is evident from the figure, the TPB model offers a significant improvement over the TIB model as there is a clear agreement with the MC average in the both the average and cumulative average disease prevalence for both choices of parameters. The reason for such a significant improvement can be explained by the fact that the TPB model is exact on networks with no non-backtracking cycles. However, because the cattle trade network is a production network, there exist very few non-backtracking cycles making the network structure highly tree-like in its supra-adjacency embedding. This can be explained by the fact that the existence of such cycles are inefficient and cost prohibitive in the trade process. As a result the network is well approximated by a tree network and contains very few non-backtracking cycles. Therefore, the SIR process is well approximated by the TPB model for such a network. As shown in Fig. 5 the number of non-backtracking cycles as a fraction of the total number of paths in the network is very small, peaking at just over 0.025%. Hence, the network is unlikely to suffer from the echo chamber effect in the TPB model. However, there are very many reciprocal (bi-directional) links in the network, meaning many farms trade in both directions. The reason for the difference in prediction between the TIB and TPB models is that the TPB model is able to account for all these reciprocal edges.

Fig. 5
figure 5

a Proportion (bars) and number of all non-backtracking (NBT) paths (blue solid) and cycles (blue dashed) that close at each respective time step. b Final outbreak sizes for varying values of \(\beta\) with \(\mu =0.005\). The vertical line shows the critical probability \(\beta _{\text {crit}}\) obtained from the TPB model. The average of the Monte-Carlo (MC) is shown as green curve

Next, we test our analytical findings for the epidemic threshold on the cattle trade network. Figure 5b depicts the average final outbreak sizes of a number of realisations against increasing values for \(\beta\) while keeping \(\mu\) fixed at 0.005. The critical \(\beta\) which gives rise to an epidemic according to the analytically computed epidemic threshold for such a fixed \(\mu\) in the TPB model is given as \(\beta _{\text {crit}}\approx 0.025\). For values of \(\beta\) that are greater than the computed epidemic threshold, there is an obvious but gradual change in dynamics as local outbreaks no longer die out, but now propagate throughout the network leading to larger final outbreak sizes as the value for \(\beta\) gets larger, thus showing agreement with the analytical result for the epidemic threshold. Overall, we find that such trade networks are a good candidate for the TPB model as they avoid many non-backtracking cycles.

Conference contacts

The second data set (cf. Table 1) is the Conference network described in Génois and Barrat (2018). It includes the face-to-face interactions of 405 participants at the SFHH conference held in Nice, France 2009. Each snapshot of the network represents the aggregated contacts in windows of 20s. Since this data set describes face-to-face interactions, each contact is bi-directional and so an undirected network is the natural choice to model these interactions.

Fig. 6
figure 6

a Degree of the conference network aggregated over the entire observation period. b Time series of number of active edges per time step in the network

Because of the small number of nodes in the network, it is difficult to draw detailed conclusions from the degree distribution. As shown in Fig. 6a, there is a clear heavy tail with most vertices having a relatively small aggregated degree. In Fig. 6b we see that the network activity in this case shows a number of peaks occurring then quickly dying out. These are explained by breaks between sessions at the conference during which the participants converse and interact. Because of the time scale and observation period of this particular temporal network, it is not feasible to model the spread of disease as infection and recovery is unlikely to occur within the observation period, which is approximately 20 h. However, we can use our model to simulate the spread of viral information or “gossip” using the same dynamics as the SIR model. Infection is equivalent to receiving some information in such a way that it becomes interesting enough to for the individual to try and spread to those they contact in the future and recovery is equivalent to growing tired of the information and no longer inform others they meet.

Fig. 7
figure 7

TIB (blue dashed) and TPB (green solid) SIR models on the conference network together with the average of \(10^3\) Monte-Carlo (MC) simulations (red dash-dotted). a, c The time series for the prevalence of infection in the network. b, d Depict the cumulative prevalence of infection within the network. The probability of infection is set to \(\beta = 0.3\) in a, b and \(\beta = 0.5\) in c, d. The initial chance of infection is 0.03. Other parameters: \(\mu =0.005\)

Figure 7 shows the time series of the different models for two probabilities of infection: \(\beta = 0.3\) in (a),(b) and \(\beta = 0.5\) in (c),(d). Again, one can observe that in every case the TPB approximation offers an improvement over the TIB. However, compared to the MC simulations the TPB model does not offer perfect agreement in contrast to the cattle trade network. An interesting observation is that in panel (a) we see the average prevalence of the infection, according to the TPB model, has a higher peak than the TIB model and it may appear as though the TIB model performs better. However, when we look at the corresponding plot in panel (b), which shows the cumulative average prevalence of the disease, the TPB model is closer to the MC average at every time step. The reason the peak is bigger for the TPB model is that the TIB model has a sustained higher first peak, which uses up the pool of susceptible individuals leaving a smaller population available to catch the disease.

Fig. 8
figure 8

a Proportion (bars) and number of all non-backtracking (NBT) paths (blue solid) and cycles (blue dashed) that close at each respective time step. b Final outbreak sizes for varying values of \(\beta\) with \(\mu =0.005\)

The reason we do not see a good agreement with the MC average for this particular data set is due to the underlying topology of the network that is a physical social interaction network where individuals congregate in groups and most or all in the group interact with one another. This leads to large clusters that give rise to many non-backtracking cycles. The more time-respecting non-backtracking cycles that occur, the worse the TPB model will perform. It is for this reason that we see a relatively large deviation from the MC simulations for the TPB model. This can be explained by Fig. 8a, which in contrast to the cattle network shows that the number of NBT cycles is relatively dense at many points in time, reaching as high as 12%.

In Fig. 8b we see the distribution of final outbreak proportions against the critical \(\beta\) computed for the epidemic threshold of the TPB model. Again, for values of \(\beta\) that are greater than the computed epidemic threshold, there is an obvious but gradual change in dynamics as local outbreaks no longer die out, but now propagate throughout the network leading to larger final outbreak sizes as the value for \(\beta\) gets larger, thus showing agreement with the analytical result for the epidemic threshold. However, the agreement with the final outbreak sizes only remains consistent with the MC average for values of \(\beta\) below and slightly above the epidemic threshold.

Conclusions

In this paper, we have presented work done on SIR pair-based models by systematically extending them to a temporal setting and investigating the effect of non-backtracking cycles on the accuracy of the model on arbitrary network structures. We have found that the existence of many such non-backtracking cycles leads to a deviation in the pair-based model from the true SIR process due to the echo chamber effect they induce. Thus, the pair-based model is best suited to network structures which do not contain many cycles, such as production networks. We also find that our analytical finding for the epidemic threshold holds up when compared to numerical simulations, by showing a qualitative change in the final outbreak proportion.