Epidemic threshold in pairwise models for clustered networks: closures and fast correlations

The epidemic threshold is probably the most studied quantity in the modelling of epidemics on networks. For a large class of networks and dynamics, it is well studied and understood. However, it is less so for clustered networks where theoretical results are mostly limited to idealised networks. In this paper we focus on a class of models known as pairwise models where, to our knowledge, no analytical result for the epidemic threshold exists. We show that by exploiting the presence of fast variables and using some standard techniques from perturbation theory we are able to obtain the epidemic threshold analytically. We validate this new threshold by comparing it to the threshold based on the numerical solution of the full system. The agreement is found to be excellent over a wide range of values of the clustering coefficient, transmission rate and average degree of the network. Interestingly, we find that the analytical form of the threshold depends on the choice of closure, highlighting the importance of model selection when dealing with real-world epidemics. Nevertheless, we expect that our method will extend to other systems in which fast variables are present.


Introduction
Epidemic dynamics on networks, being susceptible-infected-susceptible (SIS), susceptible-infected-recovered (SIR) or otherwise, are often modelled as continuous time Markov chains with discrete but extremely large state spaces of order m N , where m denotes the number of different disease statuses (e.g. m = 2 for SIS and m = 3 for SIR) and N stands for the number of nodes in the network. This makes the analysis of the resulting exact system almost impossible, except for some specific network topologies such as the fully connected network, networks with considerable structural symmetry or networks with few nodes (Kiss et al. 2017;Holme 2017).
Often, this problem is dealt with by focusing on mean-field models where the goal is to derive, often heuristically, a system of ordinary or integro-differential equations that describe (non-Markovian) epidemics for some average quantities, such as the expected number of nodes in various states, the expected number of links in various states or the expected number of star-like structures (focusing on a node and all of its neighbours). These methods usually rely on closures to break the dependency on higher-order moments (e.g. the expected number of nodes in a state depends on the expected number of links in certain states and so on). Such an approach has led to a number of models including heterogeneous or degree-based mean-field (Pastor-Satorras and Vespignani 2001;Pastor-Satorras et al. 2015), pairwise (Rand 1999;Keeling 1999), effective-degree (Lindquist et al. 2011), edge-based compartmental (Miller et al. 2012) and message passing (Karrer and Newman 2010a), to name a few. These models essentially differ in the choice of variables over which the averaging is done. Perhaps the most compact model with the fewest number of equations is the edge-based compartmental model (Miller and Volz 2013) which is valid for heterogeneous networks with Markovian SIR epidemics, although extensions of this model for arbitrary infection and recovery processes are possible (Sherborne et al. 2018).
Pairwise models have been extremely popular and the very first model for regular networks and SIR epidemics (Rand 1999;Keeling 1999) has been generalised to heterogeneous networks (Eames and Keeling 2002), preferentially mixing networks (Eames and Keeling 2002), directed (Sharkey et al. 2006) and weighted networks (Rattana et al. 2013), adaptive networks (Gross et al. 2006;Kiss et al. 2012;Szabó-Solticzky et al. 2016), and structured networks (House et al. 2009) among others. Perhaps this is due to the relative simplicity and transparency of the pairwise model, whereby variables have a straightforward interpretation and a basic understanding of the network and epidemic dynamics coupled with good bookkeeping leads to valid and analytically tractable model equations. Pairwise models have been successfully used to derive analytically the epidemic threshold and final epidemic size, with these results mostly limited to networks without clustering. The propensity of contacts to cluster, i.e. two neighbours of a node being neighbours of one another, is known to lead to many complications, and modelling epidemics on clustered networks using analytically tractable mean-field models is still limited to networks with specific structural features (House et al. 2009;Newman 2009;Miller 2009a, b;Karrer and Newman 2010b;Volz et al. 2011;Ritchie et al. 2016). However, using approaches borrowed from percolation theory (Miller 2009b) and focusing more on the stochastic process itself (Trapman 2007a), some results have been obtained. For example, Miller (2009b) showed that for the SIR epidemic on clustered networks with heterogeneous degree distributions, the basic reproduction number is given by where k i stands for the ith moment of the degree distribution, T is the probability of infection spreading across a link connecting an infected to a susceptible node and n denotes the average number of triangles that a node belongs to. The first positive term in Eq. (1.1) corresponds to the threshold for configuration-type networks without clustering. The second term in Eq. (1.1), which is negative, shows that clustering reduces the epidemic threshold when compared to the unclustered case, the contribution of the remaining terms being of a smaller order. For pairwise models, clustering first manifests itself by requiring a different and more complex closure, which makes the analysis of the resulting system, even for regular networks and SIR dynamics, challenging. Furthermore, it turns out that such a closure may in fact fail to conserve pair-level relations and may not accurately reflect the early growth of quantities such as closed loops of three with all nodes being infected (House and Keeling 2010). Such considerations have led to an improved closure being developed in an effort to keep as many true features of the exact epidemic process as possible (House and Keeling 2010). In this paper we focus on the classic pairwise model for regular networks with clustering, using both the simplest closure and a variant of the improved closure. We show that by working with two fast variables, corresponding to correlations between neigbouring nodes during the epidemic, we can analytically determine the epidemic threshold as an asymptotic expansion in terms of the global clustering coefficient φ, defined in Sect. 2.1.
The use of fast variables is not new (Keeling 1999;Juher et al. 2013;Llensa et al. 2014;Britton et al. 2016;Eames 2008). However, in many cases the epidemic threshold has only been obtained numerically and it was framed in terms of a growth-rate-based threshold, which is equivalent to the basic reproduction number at the critical point. Eames (2008) considered a hybrid pairwise model incorporating random and clustered contacts, with the analysis focusing on the growth-rate-based threshold. Eames (2008) derived a number of results, some analytic (the critical clustering coefficient for which an epidemic can take off) and some semi-analytic, and showed, in agreement with most studies, that clustering inhibits the spread of the epidemic when compared to an equivalent network without clustering but with equivalent parameter values governing the epidemic process. However, no analytic expression for the epidemic threshold was provided.
More recently, Li et al. (2018) calculated the epidemic threshold in a pairwise model for clustered networks with closures based on the number of links in a motif, rather than nodes. This led to where n is the average number of links per node, φ is the global clustering coefficient, and τ and γ are the infection and recovery rates, respectively. The expression above can be expanded in terms of the clustering coefficient φ to give which again demonstrates that clustering reduces the epidemic threshold. Building on these results, and effectively extending the work by Keeling (1999) and Eames (2008), our paper presents a method to determine the epidemic threshold analytically and applies it in the context of pairwise models with two different closures for clustered networks. The paper is structured as follows. In Sect. 2 we outline the model with closures for unclustered and clustered networks discussed in Sect. 3. In Sect. 4 we briefly review existing results and approaches for the pairwise model with the simple closure and then focus on the correlation structure in terms of fast variables, showing that the epidemic threshold can be expressed via the solution of a cubic polynomial. This key solution is determined numerically and analytically as an asymptotic expansion in terms of the clustering coefficient. In Sect. 5 we show that our approach extends to a compact version of the improved closure, thus validating and generalising our approach. Finally, we conclude with a discussion of the results, including comparing the threshold to other known results and touching upon a number of possible extensions.

The network
We begin by considering a population of N individuals with its contact structure described by an undirected network with adjacency matrix G = (g i j ) i, j=1,2,...,N where g i j = 1 if nodes i and j are connected and zero otherwise. Self-loops are excluded, so g ii = 0 and g i j = g ji for all i, j = 1, 2, . . . N . The network is static and regular, such that each individual has exactly n edges or links. The sum over all elements of G is defined as ||G|| = i, j g i j . Hence, the number of doubly counted links in the network is ||G|| = n N . More importantly, using simple matrix operations on G, we can calculate the global clustering coefficient of the network where trace(G 3 ) yields six times the number of closed triples or loops of length three (uniquely counted) and ||G 2 || − trace(G 2 ) is twice the number of triples (open and closed, also uniquely counted).

SIR dynamics
The standard SIR epidemic dynamics on a network is considered. The dynamics are driven by two processes: (a) infection and (b) recovery from infection. Infection can spread from an infected and infectious node to any of its susceptible neighbours and this is modelled as a Poisson point process with per-link infection rate τ . Infectious nodes recover from infection at constant rate γ .

The unclosed pairwise model
Let (2.5) (2.6) We note that Eqs. (2.4)-(2.6) contain triples but evolution equations for these are not given. To determine solutions of the system, we must find a way to account for these triples in terms of pairs and singles, a method referred to as closing the system. The system above is exact before a closure is applied. This means that it can be derived directly from the exact stochastic epidemic model on the network, given by a continuous time Markov Chain, without making any approximations [a precise proof for the SIS epidemic was given by Taylor et al. (2012)]. The flow between compartments and the rates of the SIR pairwise model are illustrated in Fig. 1. The system given above only contains dynamically relevant variables, i.e. those that emerge naturally but following a strict bookkeeping rule, and those that appear when a chosen closure for the triples is considered.

Closures
A quick inspection of the unclosed pairwise system (2.2)-(2.6) reveals that only triples of type [AS I ] need closing, with A ∈ {S, I }. These triples, as well as triples of type [RS I ], are illustrated in Fig. 2   for the SIR pairwise model. In the compartments of pairs, straight arrows denote infections coming from within the pair (with a rate depending on a pair) or from outside the pair (with a rate depending on a triple), and curved arrows denote a recovery. The colour indicates the status of the "first" node in the pair. Symmetry allows us to conclude that some of the variables (see lighter shaded variables on the right hand side of the pairs diagram) must equal their symmetric version (e.g. [RS] = [S R]), so we do not need to directly calculate both quantities General setup for a central susceptible node with a given infected neighbour for a unclustered and b clustered regular networks with degree n. Dashed arrows indicate that the infected node may be connected to the other neighbours of the central susceptible node. Random variables X 1 , X 2 , . . . , X n−1 take values from the set {S, I , R}

Closure for unclustered networks
First, we consider the situation depicted in Fig. 2a. We aim to find an approximation for the distribution of the random variables X i which take values from the set {S, I , R}. Several observations can be made. The expected number of A − S type links is [AS] and the total number of links emanating from susceptible nodes counted across the whole network is n [S]. Hence, the most straightforward approximation would be to assume that X i , with i = 1, 2, . . . , n − 1, are independent and identically Bernoulli distributed random variables with probability p uc A|S−I = [AS] n [S] , where p uc A|S−I stands for the probability that a neighbour of a susceptible node already connected to an infected node will be in state A, provided that the network is unclustered. Averaging across the whole network leads to the closure It is important to note that the new closed system, obtained upon using Eq. (3.1) in system (2.2)-(2.6), is effectively an approximation of the exact pairwise model (2.2)-(2.6) and one should question if closure (3.1) conserves the properties of the stochastic process and those of the counting on the network. For example, it is expected that in the closed system the number of nodes is conserved, i.e.

Simple closure
The presence of closed loops of length three, as illustrated in Fig. 2b, introduces some complications. Namely, a neighbour of the central susceptible node that is itself connected to an infected neighbour of the central node is less likely to be susceptible due to the added pressure from the infected neighbour, when compared to the case when the force of infection is distributed evenly, as it is the case for the closure for unclustered networks (3.1). More precisely, the epidemic process on the network displays clear correlations. Cator and Van Mieghem (2014) have shown that the exact SIS and SIR epidemics on networks are non-negatively correlated in the sense that P(I i I j ) ≥ P(I i )P(I j ). Here, P(I i I j ) represents the probability that nodes i and j, connected by a link, are both infected, while P(I i ) stands for the probability of node i being infected. For this result to hold, all processes must be Markovian and infection rates across all links and recovery rates of all nodes have to be fixed a priori. Using the pairwise model for an S I S epidemic on an unclustered network with closure (3.1), it has been shown that the same correlation is preserved when averaging at the population level (Kiss et al. 2017). While the proof has not been extended to the pairwise SIR model, intuitively we expect to find the same correlation structure. Based on these observations, we assume that the correlation structure in exact SIS and SIR epidemics on networks averaged at the population level is maintained. Hence, the inequalities Intuitively, this means that as the epidemic spreads on the network, infected nodes are more likely to have neighbours which are themselves infected (either those that infected the node or were infected by it), and at the 'front' of the epidemic we would expect to observe a 'sea' of susceptible nodes alongside a 'front' of links between susceptible and infected nodes that drives the epidemic. Hence, clustering and correlations need to be accounted for and a new p c A|S−I for clustered networks needs to be defined. This has been done by Keeling (1999) [see also work by Rand (1999) and Keeling et al. (1997)] and relies on a correlation factor, C AB , that is able to capture the propensity that two nodes connected by a link are in states A and B, respectively. This is given by n[S] C AI . However, before the closure can be written, open and closed loops need to be treated separately. In order to do this, we split the closure based on whether the neighbour whose state is to be determined is part of a closed loop of three nodes and thus in direct contact with an infectious node, or not. This leads to where φ is defined in Eq. (2.1). With this in mind, the closure can be derived by averaging equation (3.1) over the unclustered and clustered parts of the network. This leads to (3.7) This same closure has been derived by Keeling et al. (1997) and Keeling (1999). Framing p uc A|S−I and p c A|S−I more generally and independently of the network type, i.e. simply considering p A , the following statement holds:

Improved closure
We note that while p uc A|S−I satisfies the above proposition, p c A|S−I does not. In particular, we find However, for the clustered part of the network this is not the case. We find that which does not result in the desired (n − 1)[S I ]. This can be corrected in a straightforward way by defining Hence we can now write as required. It is informative to investigate the relationship between the various probability models that lead to different closures. This is summarised in the following proposition. While p c S|S−I has a natural interpretation (it is a simple discounted variant of the probability from the unclustered network case and takes into account the observation that if the neighbour of a central susceptible node is connected to one of the infected neighbours of the same node then it is less likely that the node in question is susceptible), the interpretation of p  Fig. 2, it is less likely to find neighbours who are susceptible compared with the unclustered network case.
Taking into account the new way of defining p c new A|S−I , the improved closure yields . (3.14) We finally note that the closures rely heavily on the assumption of how the states of the neighbours are distributed, and the assumption of independent and identically Bernoulli-distributed variables is a strong one. For clustered networks in particular, we have illustrated different ways of incorporating correlations induced by closed cycles of length three. Despite these seemingly strong assumptions, it is known that the pairwise model for unclustered networks is equivalent to the edge-based compartmental equivalent on configuration networks (Miller and Kiss 2014;Kiss et al. 2017) and the latter has been shown to be the limiting system of the stochastic network epidemic model (Decreusefond et al. 2012;Janson et al. 2014). For clustered networks we are not aware of such results.

Background
Using the simple closure for clustered networks (3.7), and writing ξ = (n−1) n , we obtain the following closed pairwise model equations describing an SIR epidemic on a clustered regular network of N individuals with degree n: (4.1) (4.4) For model equations (4.1)-(4.5), the basic reproductive ratio (R 0 ) is considered by Keeling (1999). Starting from the evolution equation of the expected number of infectious nodes leads to[ where C S I is defined in Eq. (3.3). Taking into account that τ n = β and that initially Keeling (1999) claimed that R 0 = C S I β/γ . It is important to note that this R 0 is not the classical R 0 in the sense of describing the expected number of new infections produced by a typical infectious individual when introduced into a fully susceptible population. Rather it can be thought of as a growth-rate-based threshold, and has the same properties as the classical R 0 when both are exactly one. In what follows, we will simply refer to it as R (Eames 2008;Kiss et al. 2012).
In order to determine R explicitly, Keeling (1999) considered the early behaviour of C S I and found that this variable is given by the ordinary differential equation (ODE) However, the ODE above depends on the behaviour of [I ]C I I /N and Keeling (1999) showed that (4.7) Considering the quasi-equilibrium of C S I , referred to as C * S I , in Eq.(4.6) together with the expression for [I ]C I I /N in Eq. (4.7), one finds that C * S I is given by Hence, R can be calculated as C * S I β/γ , at least numerically. Variables such as C S I and C I I describe the correlations between the states of neighbouring nodes on the network as the epidemic unfolds and these have been studied numerically by Keeling (1999).
For model equations (4.1)-(4.5) and when there is no clustering in the network (φ = 0), a further simplification of Eq. (4.8) can be achieved (Keeling 1999). To determine R = C * S I β/γ in this case, one can simply solve to find C * S I = n−2 n and thus R = (n−2)τ γ .
Unfortunately when φ = 0, according to our knowledge, the quasi-equilibrium values can only be determined numerically via Eq. (4.8). In what follows, we show that by working with two new variables, α = [S I ]/[I ] and δ = [I I ]/[I ], which are still closely linked to the correlations formed during the spreading process, it is possible to obtain the epidemic threshold as the solution of a cubic equation and, more importantly, we show that this solution can be approximated by an asymptotic expansion in powers of φ.

Epidemic threshold
Consider the initial phase of an infection invading an entirely susceptible population in the pairwise model, described by Eqs. (4.1)-(4.5). We find thaṫ We know the quantity γ [I ] remains non-negative regardless of time in the epidemic process, and we choose to consider the epidemic threshold in terms of When R > 1 an epidemic will occur, and when R < 1 the epidemic will die out. Although we know the values of τ and γ , to determine if an epidemic will occur a priori, we require further knowledge about the quantity [S I ] [I ] at some initial time close to t = 0. At t = 0 or at the disease-free steady state, both [S I ] and [I ] are zero and hence their ratio is ill-defined. Furthermore, gaining knowledge about [I ] and this term suffers from the same problem, being ill-defined at t = 0. While this is similar to the approach taken by Keeling (1999), we focus on the variables [S I ] [I ] and [I I ] [I ] , and we motivate our choice below. The problem of finding the epidemic threshold can be dealt with in at least two other different but equivalent ways. First, one can carry out a simple linear stability analysis of the disease-free steady state as is shown in Appendices B and C. Second, the threshold can also be computed as the largest eigenvalue of the next generation matrix, see Sect. 6. However, in both cases, the variables [S I ]/[I ] and [I I ]/[I ] turn out to play a key role and their values for small times need to be determined.

Fast variables with the simple closure
To circumvent the problem of the ill-defined variables above, we exploit the fact that α := [S I ] [I ] and δ := [I I ] [I ] are fast variables when compared to the time course of the epidemic. Figure 3 shows clearly that α and δ are fast compared to the epidemic process and that they quickly converge to a quasi-equilibrium. Hence, at early times, α and δ attain their quasi-equilibrium values, and these are the values that can be used to compute the epidemic threshold.
We continue by deriving differential equations for the variables α = We note that this approach has already been exploited by Juher et al. (2013), Llensa et al. (2014) and Britton et al. (2016), with the authors focusing on combinations of SIS, SIR and SEIR models without demography and rewiring of S − I links to S − S links. In all these papers, systems of fast variables are derived and analysed in detail to gain information about the epidemic threshold and the stability of the disease-free or endemic steady states.

Global stability of the steady state
The analysis of system (4.11)-(4.12) is carried out in detail by Trapman (2007b) (see Appendix A of this paper). The only caveat there is that the global stability of the unique steady state was not confirmed, leaving the possibility of the existence of a limit cycle. Below, we sketch the main ideas of the proof and provide some extra results by using the Bendixson criterion. The starting point is to show that system (4.11)-(4.12) admits a unique steady state which is biologically meaningful First we show that the trajectories of the system remain in D for all appropriate initial conditions and all positive times. When δ = 0, then dδ/dt = 2τ α > 0. When α = 0, then dα/dt = 0. However, by condition (4.15), Finally, we need to show that if α + δ = n then d(α + δ)/dt < 0. By substituting δ = n − α, and after some algebra we obtain that d(α + δ)/dt = γ (α − n) − τ (n − 1)φα(1 − α/n) 2 < 0. The observations prove that D is invariant. A typical picture of the phase diagram is given in Fig. 4.
It turns out that both null clines can be written conveniently with α being the independent and δ being the dependent variable. The null clines corresponding to dα/dt and dδ/dt are given by (4.14) Several observations can be made. If ξ n(1−φ)−1 > 0, then δ 1 (α) will be a decreasing function in α and its intersection with the horizontal axis is at α 1 = ξ n(1−φ)−1 1−ξφ , which happens to be less than n. Furthermore, it is straightforward to show that dδ 2 (α)/dα > 0, which means that δ 2 (α) is an increasing function in α. Given the behaviour of the null clines at α = 0, it follows that their intersection gives rise to a unique steady state.
This is the same as found by Keeling (1999) in the limit of β = τ n large and when assuming that at the threshold C S I = γ /β. This condition can also be derived directly from Eq. (4.21) by replacing α = τ/γ (which corresponds to the threshold R = τ α γ ) and then taking the limit of large τ . In fact, when φ > (n − 2)/(n − 1) the disease dies out. Hence, the two null clines define a unique point of intersection as long as the condition above, (4.15), is met. The same argument holds even if the singularity of the second null cline happens to be in (0, n). However, we must also exclude the possibility that the intersection point will lie outside D. For example if the δ 2 null cline lies to the left of δ 1 then the δ 2 null cline may cross the α + δ = n boundary at a smaller value of α than the δ 1 null cline does. However, this cannot happen because, Fig. 4 Illustration of the typical phase plane of system (4.11)-(4.12). The null clines δ 1 (dashed) and δ 2 (dash-dotted), and the α + δ = n (continuous) line are plotted together with a typical trajectory ( ) that is attracted to the unique steady state of the system. Parameter values are N = 10,000, n = 5, φ = 0.5 and τ = γ = 1 0 1 2 3 4 5 0 1 2 3 4 5 in such a case, D would not be invariant since the solutions would leave D across the region of this boundary limited by the two null clines, which contradicts that fact that d(α + δ)/dt < 0 on this boundary.
Provided that condition (4.15) holds, Fig. 4 shows that a unique steady state exists. Trapman (2007b) showed that this steady state is locally stable but global stability was not confirmed. Here, in addition to the results by Trapman (2007b) we show that under certain assumptions the existence of a limit cycle can be ruled out by applying the Bendixson criterion. This also ensures the global stability of the unique steady state. Dividing the equations by α, the divergence of the system takes the form: We separated the above expression into the non-clustered and clustered parts of the network. It is obvious that when φ = 0 then B(α, δ) < 0 and thus no limit cycle can occur. Now setting B(α, δ) = 0 and neglecting the −γ /α term defines the following curve This intersects the horizontal axis at α B = n ξφ − n/2. Given that the slope of δ B is positive, the divergence will remain negative in D as long as the intersection with the horizontal axis is beyond n. This requires that Rearranging this, we obtain φ < 2n 3(n − 1) .
Hence, if the above holds then the unique steady state is globally stable. It is worth noting that if then the global stability also holds for all φ < (n − 2)/(n − 1), and as long as n < 6 the above inequality holds. Numerical tests suggest that global stability holds for all reasonable parameter values. For example, if the point of intersection of δ B with the horizontal axis is in (α β , n), then the non-existence of the limit cycle can be shown as follows. To the left of δ B the divergence is negative and the lower right quadrant of D is repellent.

Fast variables with clustering
When clustering is present in the network, the differential equations for α and δ are more complex and thus steady states are harder to compute. Firstly, we set Eq. (4.11) equal to zero and rearrange to isolate δ, finding (4.20) Plugging Eq. (4.20) into Eq. (4.12) leads to the following cubic equation in α: The solution of the cubic equation (4.21) provides the steady state(s) of system (4.11)-(4.12), and allows the computation of the threshold via the formula R c = τ α * γ . We note that the steady state in α has to be biologically plausible. α = [S I ] [I ] restricts the steady state to be positive and to be less than n, since the average number of susceptible neighbours averaged over all infected nodes needs to be less than the average degree.
(4.25) Equation (4.25) leads to which, after substituting α 0 = (n − 2) and ξ = (n−1) n yields α 1 = −2(n − 1) n 2 2τ (n − 1)(n − 2) + γ n τ (n − 2) + γ . (4.26) To summarise, we have determined the first two coefficients α 0 and α 1 of the asymptotic expansion (4.22) which solves the cubic equation (4.21). Hence, the true solution is approximated by: We make several remarks. First, the epidemic threshold will be given by R c = τ α/γ . Second, the coefficient of the first order correction of α can be rearranged in terms of R = τ (n−2) γ , the threshold for the case of unclustered networks, leading to where a = 2(n − 1)/n and where terms in φ of order larger than one have been neglected. Finally, it is clear that due to the first order correction being negative, we have that (4.29) The goodness of the estimate for α (4.27) is tested by comparing it to the numerical solution of the cubic equation (4.21). This is done in Fig. 5 for five different values of the clustering coefficient. The asymptotic approximation performs well and only breaks down for values of clustering larger than φ 0.3. From the same figure it is clear that higher values of clustering push the critical R c = 1 curve to higher values of τ and n. Hence, in the presence of clustering a viable epidemic requires either a denser network or a higher transmission rate, noting that the transmission rate and the recovery rate γ are not strictly independent.

Numerical examples
In the previous section we have demonstrated that for the pairwise model with the simplest closure for clustered networks, the determination of the epidemic threshold involves the solution of a cubic equation. While this solution can be obtained numerically, we presented an asymptotic approximation of the solution in terms of powers of the clustering coefficient φ. In Fig. 5 we present a systematic test of the newly determined threshold by comparing the threshold based on the numerical solution of the cubic equation (4.21) (continuous line in the (τ, n, 0) plane), the asymptotic approximation of the solution to the cubic equation (4.27) (dashed line and markers-•) and the numerical solution of the full ODE system corresponding to the closed pairwise model (4.1)-(4.5).
The agreement between the explicit numerical solution of the closed pairwise system and threshold based on the numerical solution of the cubic equation is excellent for all clustering values and other parameter combinations. Moreover, the agreement of these results with the threshold based on the asymptotic approximation is also excellent and remains valid for values of 0 ≤ φ ≤ 0.3. The initial conditions for the closed pairwise systems were set in the following way:

Results for the pairwise model with the compact improved closure
Starting from the improved closure (3.14) but in line with Proposition 2, we adapt the closure so that the term responsible for the approximation on the clustered part of the network does not consider variables, singles or pairs involving the R class. This leads to the new closure which we refer to as the compact improved closure. Plugging Eq. (5.1) into the exact system (2.2)-(2.6) leads to a self-consistent system that is written out in full in Appendix A.
In line with our procedure so far, we aim to find the epidemic threshold of this new pairwise system with the compact improved closure. It turns out that the approach used for the pairwise system with the simple closure is applicable to this case, and the steps and results are summarised below.

Fast variables with the compact improved closure
As we have shown before, finding the threshold relies on finding the quasi-equilibrium of α = [S I ] [I ] . In Appendix A we show that this requires knowledge about the behaviour of δ = [I I ] [I ] variable and indeed a system of differential equations involving these two variables can be derived. This system is given below As previously, the steady states of this system are of interest and apart from the trivial (α * , δ * ) = (0, 0) steady state, the quasi-equilibrium can be found by first expressing δ as a function of α. This can be done by setting Eq. (5.2) equal to zero and rearranging, leading to Plugging Eq. (5.4) into Eq. (5.3) and collecting powers of δ leads to the following cubic equation where A = (n − 2) − 2φ(n − 1) and B = γ /τ . It is worth noting that in this case it is easier to work with δ, but any results can be converted in terms of α which is the main variable of interest. However, before we proceed with the asymptotic expansion of the solution, we show that there is a unique biologically feasible steady state.
To continue we focus on showing that (5.2)-(5.3) admits a unique steady state which is biologically meaningful, i.e. (α * , δ * ) ∈ D. The null cline corresponding to dα/dt can be rewritten to give It is straightforward to check that meaning that the function is decreasing for all α. Setting α = 0 in (5.6) leads to δ = n(n − 2)/(2φ(n − 1) − (n − 2)), which can be both negative or positive. On the other hand setting δ = 0 in (5.6) yields α = (n − 2). This null cline has a singularity at α * = (n − 2) − 2φ(n − 1), with α * < (n − 2) < n. If α * < 0 then the branch on the left of the vertical asymptote will not intersect D. This happens exactly when 2φ > (n − 2)/(n − 1). So in this case the branch of the null cline to the right of the asymptote intersects the α-axis at ((n − 2), 0) and the δ-axis at (0, n(n − 2)/(2φ(n − 1) − (n − 2))), where the intersection with the δ-axis happens at a positive value, namely n(n − 2)/(2φ(n − 1) − (n − 2)) > 0, and this inequality holds true due to requiring that α * is negative. This point may be greater than n but also intersects the horizontal axis at (n − 2, 0). This is illustrated in Fig. 6 (left panel). When the singularity point is positive, α * > 0, that is when 2φ < (n − 2)/(n − 1), then the intersection with the δ-axis happens at a negative value of δ. This is also illustrated in Fig. 6 (right panel), where the positive singularity is clearly visible with the intersection with the δ-axis being out of the range of the plot.
The null cline corresponding to dδ/dt is given by This null cline passes through (α, δ) = (0, 0) and the derivative of α n (δ) is always positive, namely, The denominator is a quadratic polynomial in δ with the discriminant being always positive and thus leading to two distinct real roots. From the equation it follows that sum of the roots is (n − 2) − 2φ(n − 1) and their product is −2n < 0. Therefore, two singularity points exist, one for negative and the other for positive δ. α n (δ) is such that Hence, D happens to lie, at least partly, in the area defined by the two singularity points (i.e. the region between the two vertical asymptotes if considered in the (δ, α) plane). In this area the null cline increases with δ starting from (α, δ) = (0, 0), see both panels in Fig. 6. Summarising, we have shown that the null clines will intersect at a unique point, and this point cannot be outside D due to the orientation of the vector fields, see also the argument presented in Sect. 4.3.1. Finally, we show that the existence of a limit cycle can be ruled out by applying the Bendixson criterion. This also ensures the global stability of the unique steady state. Dividing Eqs.
It is easy to show that this is negative. Even if − γ α is neglected, we have that since (n + δ) is greater than both n and (n − 1).
(5.13) Substituting (5.11) into (5.13) and collecting terms of order φ 0 yields . (5.17) Following the same process to collect terms of order φ 1 , we find which can be rearranged to yield with δ 0 defined in (5.17). In summary, we have determined the first two coefficients δ 0 and δ 1 of the asymptotic expansion for δ given in Eq. (5.11). Hence, the true solution is approximated by the following expression: Finally, we are able to plug (5.20) into the quasi-equilibrium point for α, given in Eq. (5.4), to obtain which, upon neglecting terms in φ of order larger than one, can be rearranged to find The expression for α (5.22) can be used to determine the epidemic threshold as follows It is straightforward to see that again R cci ≤ R, with clustering making the spread of the epidemic less likely.

Numerical examples
In Fig. 7 we repeat the systematic test of comparing the epidemic threshold generated via the numerical solution of the cubic equation (5.5), the epidemic threshold generated by the asymptotic expansion (5.23) and the numerical value of the final epidemic size predicted by the pairwise model with the compact improved closure, over a wide range of (τ, n) values. Several observations can be made. First, it is clear that higher values of clustering push the location of the threshold to higher τ and n values, meaning that the limiting effect of clustering on the epidemic spread can only be overcome if either the value of the transmission rate or average degree increases. Second, the agreement between the threshold based on the numerical solution of the cubic equation (5.5) and the asymptotic expansion (5.20) is excellent over a wide range of φ values. In fact, in this case the agreement is excellent for 0 ≤ φ ≤ 0.45, with only small deviations even for φ = 0.6. The agreement between the numerical solution of the pairwise model and the threshold based on the numerical solution of the cubic equation (5.5) remains excellent across all parameter values.

Comparing epidemic thresholds based on different models
Exploiting the presence of fast variables and combining this with elements of perturbation theory allowed us to compute the epidemic threshold for the pairwise model with two different closures that take clustering into account. Our results are in line with the findings by Li et al. (2018) and Miller (2009b). Li et al. (2018) calculated the epidemic threshold in a pairwise model for clustered networks with a closure based on the number of links in a motif, rather than nodes. This led to (6.1) Equation (6.1) can be expanded in terms of φ to give which again reflects our finding that clustering reduces the epidemic threshold.
Similarly but for clustered networks with heterogeneous degree distributions, Miller (2009b) found that where k i stands for the ith moment of the degree distribution, T is the probability of infection spreading across a link connecting an infected to a susceptible node and n denotes the average number of triangles that a node belongs to. The expression above again shows that clustering reduces the epidemic threshold when compared to the unclustered case. Furthermore, if the network is regular and we assume that infections and recoveries are Markovian processes with rates τ and γ respectively, giving T = τ/(τ + γ ), R 0 above reduces to where we have used the fact that a global clustering coefficient of φ translates to a node on average being part of 1 2 n(n − 1)φ uniquely counted triangles. This in turn coincides with Eq. (6.2). This is perhaps unexpected since the first expression was obtained based on a new type of closure for pairwise models while the other expression was based on percolation theory type arguments. Trapman (2007a) considered specific networks with household structure to investigate the effects of clustering and infectious period distribution on a modified version of R 0 referred to as R * , and lower and upper bounds for the value of this quantity were found. Similarly Ball et al. (2010) considered a random network incorporating household structure and provided the basic reproduction number which takes into account within household and global contacts.
However, as elaborated upon in Sect. 4.1, the R threshold that we compute is a growth-rate-based threshold and whilst at the threshold R = 1 ⇐⇒ R 0 = 1, R does not have the same biological interpretation as R 0 . Despite this, our analysis confirms that clustering starves the spreading epidemic of susceptible neighbours such that the epidemic is less likely to spread if the networks are clustered, all other parameters being equal. More importantly, the epidemic threshold is model-dependent and the pairwise model with the compact improved closure leads more readily to epidemic outbreaks when compared to the pairwise model with the simple closure, see Figs. 5 and 7. While this ordering is true for the parameters used in this paper, we cannot conclude that this ordering holds true for all parameter values. Further research may focus on the ordering of these thresholds and gaining a better understanding of the impact of model choice on the values of the epidemic threshold.
The computation of the true R 0 for pairwise models can be attempted by considering the next generation matrix approach (Van den Driessche and Watmough 2002). Looking at the pairwise model with the simplest closure and ordering the variables involved in the spreading process as: [I ],[S I ], the generation of new infectious cases at the the disease-free steady state is given by where the lower right term is obtained from Eq. Now all other transfers between compartments are summarised in the V matrix, which is given below where the lower right term describes the rate at which [S I ] pairs are depleted. This is obtained from Eq. (4.3) as followṡ where again all expressions were evaluated at the disease free steady state. Now R 0 is given by the leading eigenvalue of F V −1 , which is Obviously, this seems like a rather complicated expression since the quasi-equilibrium values for α and δ are needed. These are only available as asymptotic expansions in powers of φ. Nevertheless, for φ = 0, R 0 = τ (n−1) τ +γ , which agrees perfectly with the two results quoted above. Considering the φ > 0 case, we write R 0 = r 0 + φr 1 , α = α 0 + φα 1 and δ = δ 0 + φδ 1 . Plugging these into Eq. (6.7), leads to While the first term in the expansion for R 0 agrees with the results quoted above, the second term seems less likely to be equivalent to those shown above. This same approach can be used to compute R 0 when the compact improved closure is used. We believe that comparing these different ways of computing the epidemic threshold can contribute to reconciling different methods and will lead to more clarity and transparency between various modelling approaches. Finally, we report some results concerning networks composed of two layers, local within household and global contacts, where epidemic threshold-like quantities have been proposed (Ball et al. 2010). Taking the infection rates over global/network and local/household edges to be the same means that households in the model can be viewed as a device for introducing clustering into the network. This observation motivates our short analysis below. We consider the simple example of a network with all households of size three with additional global contacts assigned to nodes according to a configuration-like network with a regular degree, say μ D . This is to keep in line with our assumption of regular random networks. Based on results by Ball et al. (2010), the clustering in such a network is , (6.8) which can be inverted to give μ D in terms of clustering Assuming that both infection and recovery are Markovian with rates λ G (infection through global links), λ L (infection within households) and γ , and following the calculations by Ball et al. (2010) it is easy to show that the epidemic threshold is where Plugging in the expression for μ D , as in Eq. (6.9), leads to It is now obvious that R * decreases as φ increases, but to keep in the spirit of this section we expand the above in terms of φ. Given that around x = 0 the following expansion holds x − · · · , we can rewrite R * to give (6.13) Two important remarks can be made. First, even though R * defines an epidemic threshold, it does not have the same interpretation as the basic reproduction number: it is the household reproduction number. However, it is a threshold parameter so it takes a value below/at/above its threshold value (= 1) precisely when any other threshold parameter (such as R 0 ) is below/at/above its threshold value. Secondly, the dependency on φ for the various epidemic thresholds differs. While for most thresholds considered here this dependency is via a negative term of O(φ), the threshold from the household model decreases as O((φ) −1/2 ) as φ increases away from zero. This may indicate a clear difference in the underlying models but all models may be correct as long as their individual assumptions are met. Therefore, further exploration may focus on understanding which assumptions lead to this discrepancy and what the implications of the various modelling approaches are when applying such models in reality.

Discussion
In this paper we derived an analytic epidemic threshold using pairwise models but for clustered networks. For the unclustered case this problem has been solved previously (Keeling 1999). Here, however, by exploiting the presence of fast variables and using elements of perturbation theory, we were able to find the epidemic threshold as an asymptotic expansion in powers of the clustering coefficient. Our analysis confirms that clustering starves the spreading epidemic of susceptible nodes such that the epidemic is less likely to spread if the networks are clustered, all other parameters being equal. More importantly, the epidemic threshold is modeldependent and the pairwise model with the compact improved closure leads more readily to epidemic outbreaks when compared to the pairwise model with the simple closure, see Figs. 5 and 7. While this ordering is true for the parameters used in this paper, it is easy to show that this relation can change if parameters are tuned accordingly.
We carried out a full analysis of two systems of fast variables (one corresponding to the simplest closure with no clustering, the other corresponding to the compact improved closure for clustered networks). Both systems exhibit similar behaviours but, surprisingly, the more complicated one (that with the compact improved closure) yields results with virtually no constraints on the parameter values.
It is obvious that the complexity of the closure has a bearing on the complexity of the resulting model. As shown in the paper, using the compact improved closure leads to a more complex model whose analysis is slightly more complicated. After submitting the present paper and while waiting for the reviews, we analysed the system with the full improved closure . However, our analysis only included the asymptotic expansion of the epidemic threshold without considering the detailed analysis of the system of fast variables (e.g. existence and uniqueness of a biologically feasible steady-state). This system is now four dimensional with not two but four fast variables (the extra variables being [S R]/[R] and [I R]/[I ]). In doing so, we were able to confirm the effectiveness and generality of the approach presented in the paper.
It will also be worthwhile to compare different models in order to identify the impact of clustering on epidemics by mapping out regions in the parameter space where its effect is strongest. It is known that when the network is dense the effect of clustering is limited and the same holds when the transmission/recovery rates are high/low, respectively. Moreover, as we have shown in Sect. 6 there is scope for reconciling epidemic thresholds computed from different mean-field or stochastic models where the network is a key ingredient. More importantly, while there is some agreement between the different epidemic threshold expressions, especially in some limits or particular cases, it is clear that the epidemic threshold is model dependent. Hence, the biology of the disease and the contact pattern has to be carefully analysed and taken into account when choosing models that are to be used in relation to actual epidemics.
Of course there remains the issue of accounting for degree heterogeneity in the network and this has been addressed to some extent by using percolation type approaches. The approach that we presented in this paper may be extended to degree-heterogeneous clustered networks, but this will require more sophisticated models such as effectivedegree, or compact/super-compact pairwise models (Simon and Kiss 2015). These will no doubt lead to more complex systems which are more challenging to analyse. The simplest starting point could be to consider a network with nodes having either degree k 1 or k 2 . For ease of treatment, let N i be the number of nodes with degree k i with i ∈ {1, 2}. Now one can assume that clustering in the network is introduced at random so it is going to be proportional to the degree and the mixing between the two groups of nodes. One can assume the simplest case of proportional mixing, where the number of links between nodes of degree k i and k j , n i, j is simply n i j = k i k j N i N j l k l N l . Then, the closure could be considered as follows where S i denotes the class of susceptible nodes of degree k i . Now appropriately scaled closures for the triples are needed, which will depend on the degree of the nodes and how clustering is apportioned over nodes of different degrees. The viability of such a model will then rely on whether such closures are compact and compatible enough to derive a reasonably simple overall expression for [AS I ], ideally one where the closure no longer depends on degree, but rather such information appears as some factor in the closure.
Finally, it would be worthwhile to test our findings against explicit stochastic network simulations. Since our focus was on exploiting the presence of fast variables and the use of perturbation analysis to determine the epidemic threshold analytically, such empirical validation was thought to be beyond the scope of the present work. We hope that the results of this paper may encourage other researchers to consider and tackle the challenges posed by modelling epidemic dynamics on clustered networks with heterogeneous degree distributions. = 2τ α − γ δ + 2τ (n − 1) φαδ 1 n + δ − τ αδ. When the disease-free steady state is stable, the system will always end up at the disease-free steady state and thus no epidemic will occur. When the disease-free steady state becomes unstable, there exists (at least) a second steady state whereby an epidemic will occur and [S] will no longer be equal to N . To determine a stability condition for the disease-free steady state (B.1), we must compute the Jacobian matrix J of the system (4.1)-(4.5), evaluated at the disease-free steady state, and solve to find its eigenvalues. By computing partial derivatives of each differential equation ( [I ] is also required in Eq. (4.10), and the threshold cannot be computed without it.

C Standard linear-stability analysis for the case of the compact improved closure
To determine an epidemic threshold, we consider conditions for stability of the diseasefree steady state (B.1). To do so, we compute the Jacobian matrix evaluated at the disease-free steady state as  [I ] . The Jacobian above becomes useful once analytical expressions for α and δ are obtained (or it could be an asymptotic expansion or even numerical values). Plugging these in the Jacobian will allow to either numerically or analytically compute the threshold. We note that using the linear-stability analysis or focusing on the initial growth rate should lead to the same threshold value, as was already shown for the case of the system with the simple closure in Sect. B.