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 , where m denotes the number of different disease statuses (e.g. for SIS and 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 bywhere 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 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 towhere 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 givewhich 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.

Model formulation

The network

We begin by considering a population of N individuals with its contact structure described by an undirected network with adjacency matrix where if nodes i and j are connected and zero otherwise. Self-loops are excluded, so and for all . 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 . Hence, the number of doubly counted links in the network is . More importantly, using simple matrix operations on G, we can calculate the global clustering coefficient of the networkwhere yields six times the number of closed triples or loops of length three (uniquely counted) and 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 equal 1 if the individual at node i is of type A and equal zero otherwise. Then single nodes (singles) of type A can be counted as , pairs of nodes (pairs) of type can be counted as and triples of nodes (triples) of type can be counted as . This method of counting means that pairs are counted once in each direction, so , and [AA] is even. Using this notation to keep track of singles, pairs and triples leads to the following system of pairwise equations describing the SIR epidemic on networks: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.

Fig. 1

Flow diagrams showing the flux between compartments of singles (left) and compartments of pairs (right) 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. ), so we do not need to directly calculate both quantities

Closures

A quick inspection of the unclosed pairwise system (2.2)–(2.6) reveals that only triples of type [ASI] need closing, with . These triples, as well as triples of type [RSI], are illustrated in Fig. 2 for unclustered and clustered networks.

Fig. 2

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 take values from the set

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 which take values from the set . Several observations can be made. The expected number of 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 , with , are independent and identically Bernoulli distributed random variables with probability , where 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 closureIt 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. . Furthermore, the number of pairs of different types must sum to nN. More subtle conditions refer to the conservation of link types at node level () and pair level (), respectively. It turns out that the closure for unclustered networks (3.1) conserves these relations (Kiss et al. 2017). Finally, the validity of closures can be empirically assessed by looking at the initial growth rate of the number of open and closed triples, where the number of open triples comprised of three infectious nodes should grow differently to the number of such closed triples. Of course such subtle tests are usually preceded by direct comparisons between the numerical solution of the closed pairwise system and explicit stochastic network simulations for a range of parameters. Such tests initially focus on prevalence of infection and final epidemic size but may include expected number of pairs.

Closures for clustered networks

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 . Here, represents the probability that nodes i and j, connected by a link, are both infected, while 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 SIS 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 inequalitieshold, where [AB] and [A] with represent the expected counts of pairs and singles of the corresponding types taken from the exact model, i.e., the continuous time full Markov chain.

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 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, , that is able to capture the propensity that two nodes connected by a link are in states A and B, respectively. This is given bywhere . This effectively compares the expected number of edges of type [AB] to what its value would be if nodes were labelled at random with [A] nodes of type A and [B] nodes of type B. If , then nodes of type A and B are positively correlated, whereas if nodes of type A and B are negatively correlated, . As expected, means that nodes are effectively labelled as type A or B at random. Equation (3.2) implies thatWe can modify to reflect these observations, leading to . 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 towhere 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 toThis same closure has been derived by Keeling et al. (1997) and Keeling (1999). Framing and more generally and independently of the network type, i.e. simply considering , the following statement holds:

Proposition 1

Consider a closure of the following form . If , where A is taken over all possible states, then .

Proof

.

Improved closure

We note that while satisfies the above proposition, does not. In particular, we findHowever, for the clustered part of the network this is not the case. We find thatwhich does not result in the desired . This can be corrected in a straightforward way by definingHence we can now writeas 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.

Proposition 2

For closures applied across the clustered part of the network and assuming that the number of nodes in state R is negligible, it follows thatand

All three probabilities follow from their definitions and assuming that . Since links are negatively correlated (3.2), it follows that and as a resultWhile 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 is less obvious. A close inspection reveals that can be rewritten asHowever, combining with , as given in Eq. (3.2), leads to . Finally, using the relation in Eq. (3.12) yieldsEquation (3.13) illustrates that as expected . Again, this simply shows that for clustered networks and for the setup in 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 , the improved closure yieldsWe 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.

Results for the pairwise model with the simple closure

Background

Using the simple closure for clustered networks (3.7), and writing , we obtain the following closed pairwise model equations describing an SIR epidemic on a clustered regular network of N individuals with degree n:For model equations (4.1)–(4.5), the basic reproductive ratio () is considered by Keeling (1999). Starting from the evolution equation of the expected number of infectious nodes leads towhere is defined in Eq. (3.3). Taking into account that and that initially , Keeling (1999) claimed that . It is important to note that this is not the classical 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 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 and found that this variable is given by the ordinary differential equation (ODE)However, the ODE above depends on the behaviour of and Keeling (1999) showed thatConsidering the quasi-equilibrium of , referred to as , in Eq. (4.6) together with the expression for in Eq. (4.7), one finds that is given byHence, R can be calculated as , at least numerically. Variables such as and 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 (), a further simplification of Eq. (4.8) can be achieved (Keeling 1999). To determine in this case, one can simply solveto find and thus .

Unfortunately when , 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, and , 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 thatWe know the quantity remains non-negative regardless of time in the epidemic process, and we choose to consider the epidemic threshold in terms of . This leads to . When an epidemic will occur, and when 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 at some initial time close to . At or at the disease-free steady state, both [SI] and [I] are zero and hence their ratio is ill-defined. Furthermore, gaining knowledge about will involve and this term suffers from the same problem, being ill-defined at . While this is similar to the approach taken by Keeling (1999), we focus on the variables and , 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 [SI] / [I] and [II] / [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 and 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.

Fig. 3

Illustration of the dynamics of prevalence, [I] / N, over time (a, b), compared to that of (c, d) and (e, f) for the pairwise model with the simple (left column) and the improved (right column) closures. Parameter values are , , and

We continue by deriving differential equations for the variables and . Differentiating and and using Eqs. (4.1)–(4.5) leads toWe 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 links to 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, i.e. . First we show that the trajectories of the system remain in D for all appropriate initial conditions and all positive times. When , then . When , then . However, by condition (4.15), . Finally, we need to show that if then . By substituting , and after some algebra we obtain that . The observations prove that D is invariant. A typical picture of the phase diagram is given in Fig. 4.

Fig. 4

Illustration of the typical phase plane of system (4.11)–(4.12). The null clines (dashed) and (dash-dotted), and the (continuous) line are plotted together with a typical trajectory () that is attracted to the unique steady state of the system. Parameter values are , , and

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 and are given bySeveral observations can be made. If , then will be a decreasing function in and its intersection with the horizontal axis is at , which happens to be less than n. Furthermore, it is straightforward to show that , which means that is an increasing function in . Given the behaviour of the null clines at , it follows that their intersection gives rise to a unique steady state. Requiring that is equivalent toThis is the same as found by Keeling (1999) in the limit of large and when assuming that at the threshold . This condition can also be derived directly from Eq. (4.21) by replacing (which corresponds to the threshold ) and then taking the limit of large . In fact, when 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 null cline lies to the left of then the null cline may cross the boundary at a smaller value of than the null cline does. However, this cannot happen because, 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 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 then and thus no limit cycle can occur. Now setting and neglecting the term defines the following curveThis intersects the horizontal axis at . Given that the slope of is positive, the divergence will remain negative in D as long as the intersection with the horizontal axis is beyond n. This requires thatRearranging this, we obtainHence, if the above holds then the unique steady state is globally stable. It is worth noting that ifthen the global stability also holds for all , and as long as the above inequality holds.

Numerical tests suggest that global stability holds for all reasonable parameter values. For example, if the point of intersection of with the horizontal axis is in , then the non-existence of the limit cycle can be shown as follows. To the left of the divergence is negative and the lower right quadrant of D is repellent.

Fast variables without clustering

When clustering is negligible and hence , we find thatwhere . The steady states of the system (4.18)–(4.19) are given by and . Based on Eq. (4.10), it follows that

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 , findingPlugging 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 . We note that the steady state in has to be biologically plausible. 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.

Asymptotic expansion of the epidemic threshold

The case of can be regarded as a perturbation of the case without clustering and we thus set out to find using a perturbation method. More precisely, we seek to find the roots of the cubic polynomial, given in Eq. (4.21), in terms of an asymptotic expansion in powers of , that isPlugging (4.22) into Eq. (4.21) leads toCollecting terms of order in (4.23) and after some algebra we find that satisfies:Hence, . The other solution, is not biologically feasible since by definition is positive. As expected, corresponds to the unclustered case. Collecting terms of order in (4.23), we find a polynomial in terms of and :Equation (4.25) leads towhich, after substituting and yieldsTo summarise, we have determined the first two coefficients and 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 . Second, the coefficient of the first order correction of can be rearranged in terms of , the threshold for the case of unclustered networks, leading towhere 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 thatThe 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 . From the same figure it is clear that higher values of clustering push the critical 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.

Fig. 5

Assessing the validity of the epidemic threshold based on the asymptotic approximation (4.27) (dashed line and markers—) by comparing it to the epidemic threshold based on the numerical solution of the cubic equation (4.21) (continuous lines). In the right hand column we compare both threshold curves in the plane. In the left hand column both curves are compared to the final epidemic size based on numerical integration of the pairwise model equations with the simple closure. Parameter values are , and from top to bottom the clustering coefficients are

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 () 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 . The initial conditions for the closed pairwise systems were set in the following way: , , , and . The ODEs were run for a sufficiently long time () to ensure that the epidemic died out. It is worth noting that the correct numerical solution of the cubic equation can be chosen by keeping in mind that .

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 closurewhich 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 . In Appendix A we show that this requires knowledge about the behaviour of variable and indeed a system of differential equations involving these two variables can be derived. This system is given belowAs previously, the steady states of this system are of interest and apart from the trivial 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 toPlugging Eq. (5.4) into Eq. (5.3) and collecting powers of leads to the following cubic equationwhere and . 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.

Global stability of the steady state

It is worth considering whether the trajectories of the system governed by Eqs. (5.2)–(5.3) remain in for all appropriate initial conditions and all positive times. When , then , so the line is stationary and solutions remain in D. Moreover, on this line, which is greater than zero when . This is a condition which will resurface later when the intersection of the null clines is analysed. If , then meaning that the solution cannot leave D along the line. Finally, we need to show that if then . By substituting , and after some algebra we obtain that . These findings prove that D is invariant.

To continue we focus on showing that (5.2)–(5.3) admits a unique steady state which is biologically meaningful, i.e. . The null cline corresponding to can be rewritten to giveIt is straightforward to check thatmeaning that the function is decreasing for all . Setting in (5.6) leads to , which can be both negative or positive. On the other hand setting in (5.6) yields . This null cline has a singularity at , with . If then the branch on the left of the vertical asymptote will not intersect D. This happens exactly when . So in this case the branch of the null cline to the right of the asymptote intersects the -axis at and the -axis at , where the intersection with the -axis happens at a positive value, namely , 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 . This is illustrated in Fig. 6 (left panel). When the singularity point is positive, , that is when , 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.

Fig. 6

Illustration of the typical phase plane of system (5.2)–(5.3). The null clines (dashed) and (dash-dotted), and the (continuous) line are plotted together with a typical trajectory () that is attracted to the unique steady state of the system. Parameter values are , , (left panel), (right panel) and

The null cline corresponding to is given byThis null cline passes through and the derivative of 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 and their product is . Therefore, two singularity points exist, one for negative and the other for positive . is such thatHence, 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 , 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. (5.2)–(5.3), and computing , the divergence of the system yieldsIt is easy to show that this is negative. Even if is neglected, we have thatsince is greater than both n and .

As in Sect. 4.4, we require the roots of the cubic polynomial given in Eq. (5.5). To do so, we express as an asymptotic expansion in powers of . We substitutePlugging the expansion for (5.11) into Eq. (5.5) leads toAlternatively, substituting (5.4) into the differential equation for (5.3), setting the expression equal to zero and rearranging leads toSubstituting (5.11) into (5.13) and collecting terms of order yieldsFollowing the same process to collect terms of order , we findwhich can be rearranged to yieldwith defined in (5.17). In summary, we have determined the first two coefficients and 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 obtainwhich, upon neglecting terms in of order larger than one, can be rearranged to findThe expression for (5.22) can be used to determine the epidemic threshold as followsIt is straightforward to see that again , with clustering making the spread of the epidemic less likely.

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 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 , with only small deviations even for . 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.

Fig. 7

Assessing the validity of the epidemic threshold based on the asymptotic expansion (5.20) (dashed line and markers—) by comparing it to the epidemic threshold based on the numerical solution of the cubic equation (5.5) (continuous lines). In the right hand column we compare both threshold curves in the plane. In the left hand column both curves are compared to the final epidemic size based on numerical integration of the pairwise model equations with the compact improved closure. Parameter values are , and from top to bottom the clustering coefficients are

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 toEquation (6.1) can be expanded in terms of to givewhich again reflects our finding that clustering reduces the epidemic threshold.

Similarly but for clustered networks with heterogeneous degree distributions, Miller (2009b) found thatwhere 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 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 , above reduces towhere we have used the fact that a global clustering coefficient of translates to a node on average being part of 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 referred to as , 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 does not have the same biological interpretation as . 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 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],[SI], the generation of new infectious cases at the the disease-free steady state is given bywhere the lower right term is obtained from Eq. (4.3) by looking at the rate of growth of [SI] in terms of [SI] itself and evaluating it at the disease-free equilibrium, that isNow all other transfers between compartments are summarised in the V matrix, which is given belowwhere the lower right term describes the rate at which [SI] pairs are depleted. This is obtained from Eq. (4.3) as followswhere again all expressions were evaluated at the disease free steady state. Now is given by the leading eigenvalue of , which isObviously, 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 , , which agrees perfectly with the two results quoted above. Considering the case, we write , and . Plugging these into Eq. (6.7), leads toWhile the first term in the expansion for 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 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 . 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 iswhich can be inverted to give in terms of clusteringAssuming that both infection and recovery are Markovian with rates (infection through global links), (infection within households) and , and following the calculations by Ball et al. (2010) it is easy to show that the epidemic threshold iswherePlugging in the expression for , as in Eq. (6.9), leads toIt is now obvious that decreases as increases, but to keep in the spirit of this section we expand the above in terms of . Given that around the following expansion holds , we can rewrite to giveTwo important remarks can be made. First, even though 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 () precisely when any other threshold parameter (such as ) 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 , the threshold from the household model decreases as 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 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, 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 (Kiss et al. 2018). 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 [SR] / [R] and [IR] / [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 effective-degree, 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 or . For ease of treatment, let be the number of nodes with degree with . 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 and , is simply . Then, the closure could be considered as followswhere denotes the class of susceptible nodes of degree . 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 [ASI], 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.