Beyond clustering: mean-field dynamics on networks with arbitrary subgraph composition.

Clustering is the propensity of nodes that share a common neighbour to be connected. It is ubiquitous in many networks but poses many modelling challenges. Clustering typically manifests itself by a higher than expected frequency of triangles, and this has led to the principle of constructing networks from such building blocks. This approach has been generalised to networks being constructed from a set of more exotic subgraphs. As long as these are fully connected, it is then possible to derive mean-field models that approximate epidemic dynamics well. However, there are virtually no results for non-fully connected subgraphs. In this paper, we provide a general and automated approach to deriving a set of ordinary differential equations, or mean-field model, that describes, to a high degree of accuracy, the expected values of system-level quantities, such as the prevalence of infection. Our approach offers a previously unattainable degree of control over the arrangement of subgraphs and network characteristics such as classical node degree, variance and clustering. The combination of these features makes it possible to generate families of networks with different subgraph compositions while keeping classical network metrics constant. Using our approach, we show that higher-order structure realised either through the introduction of loops of different sizes or by generating networks based on different subgraphs but with identical degree distribution and clustering, leads to non-negligible differences in epidemic dynamics.

Introduction

Network models have revolutionised our way of thinking about complex phenomena such as the spreading of disease, information transmission and processing in the brain, and the formation and interaction of social groups. Mathematical epidemiology, in particular, has embraced and benefited greatly from the use of networks as a modelling paradigm, with examples ranging from data-driven models (Tildesley et al. 2010; Barabási and Albert 1999; Kiss et al. 2006) to theoretical models (Newman 2002; Pastor-Satorras and Vespignani 2001; Keeling 1999a; Kiss et al. 2013; Miller et al. 2011). These have been used to study the impact of different network properties on how diseases break out and spread. Network models have led to greater clarity in understanding and quantifying the impact of contact heterogeneity, preferential mixing and community structure, including households (Ball and Lyne 2001; Ball et al. 2010). Although clustering of contacts or transitivity, i.e., the propensity of nodes with a common neighbour to be connected, is pervasive in many real-world networks, it continues to pose many significant challenges to the community, both from the viewpoint of network generation and, even more so, from that of deriving well-performing approximate models.

To investigate the impact of network properties, one can either use empirical networks or synthetic ones that have been generated from theoretical network models with tunable properties (Newman 2002; Karrer and Newman 2010; Molloy and Reed 1995). Many algorithms exist for clustering, but it is generally the case that focusing on achieving a particular clustering leads to changes above and beyond those controlled by the algorithm. This can preclude correct analysis of the impact of clustering (Karrer and Newman 2010; Ritchie et al. 2014; House et al. 2009; House and Keeling 2010; Keeling 1999a; Green and Kiss 2010; Milo et al. 2002; Colomer-de Simón et al. 2013; Miller 2009; Gleeson et al. 2010; Kiss and Green 2008). When looking at the impact of higher-order structure, for example, it is important that the degree distribution and, if possible, the degree correlations remain the same between networks with different clustering. Some algorithms in this direction have been proposed (Karrer and Newman 2010; Volz et al. 2011; Miller 2009; Newman 2009) and are based on the notion of subgraphs, where clustering is achieved by mixing fully-connected subgraph types, such as fully-connected triples or quadruples, and non-fully connected subgraphs, such as overlapping triangles. Using such networks, Volz et al. (2011) have developed a low-dimensional ODE model that approximates well the expected value of a number of system-level quantities, and Karrer and Newman (2010) have provided final epidemic size results for networks built by using different mixtures of subgraphs. Furthermore, House and colleagues (2009, 2010) generalised the pairwise approach to closure at the level of all possible subgraphs involving four nodes. However, a number of outstanding issues remain. The Volz et al. model, which provides time evolution, can handle well only fully-connected subgraphs. Karrer and Newman’s approach, which combines a wider variety of subgraphs, can only characterise large-time limits. Finally, to our knowledge, House et al.’s (2009) approach has not been compared to stochastic simulations and it will perform poorly for heterogeneous networks.

In this paper, we provide a general and automated approach to deriving a set of ODEs that describe, to a high degree of accuracy, the expected values of prevalence or number of recovered individuals for networks that are generated based on an arbitrary set of subgraphs. This is achieved by a rigorous separation of the role of nodes within the subgraphs and by using the probability generating function (PGF) formalism to correctly track: (a) the distribution of subgraphs to which nodes belong and (b) the excess degree that is generalised from the classical notion of a stub of a single edge to different corner types given by subgraphs. This is a significant step forward as it allows us to: (a) accurately model and analyse dynamical processes on networks with higher-order structure, thus increasing model realism, (b) map out the impact of clustering in the classic sense, and more importantly, its impact at a higher level involving four or more nodes (Ritchie et al. 2014), and (c) provide much needed insights into the role of small subgraphs or network motifs/units in epidemiology and systems biology.

The paper is organised as follows. We first review how the probability generating function (PGFs) can be used to derive ODEs that capture epidemic dynamics on configuration model (CM) networks. Such PGF-based models operate by using the versatile properties of the PGF whereby it allows us to keep track of the fraction of susceptible individuals, their degree and excess degree. Next, we generalise the CM to the hyperstub configuration model (HCM). The HCM is a network construction algorithm that selects and connects hyperstubs as prescribed by the building blocks or subgraphs of the network, rather than at random. With a basic understanding of both the network and epidemic models, we then generalise the PGF formalism to HCM networks. This section includes a step-by-step explanation of the model derivation with examples for a particular network and a detailed presentation of the code-generating algorithm. A key component of the generalised model is to label and track the position of each and every node in all subgraphs in order to avoid any ambiguity as to the role of nodes in non-fully-connected subgraphs. We then compare our approach to state-of-the-art models that can, in principle, capture the system’s expected behaviour. Where fair comparisons are possible we show that our model displays excellent agreement with existing models, otherwise we show our model to either outperform existing models or to produce accurate results where other models fail. Finally, we use the generalised model to investigate the effect of loops/cycles as well as the impact of higher-order stucture, where global clustering is kept constant, on epidemic dynamics.

Materials and methods

In this section we consolidate and generalise existing work centred around deriving low dimensional, deterministic and approximate ODEs that capture the time evolution of epidemic dynamics on configuration model networks. First, we re-introduce the basic susceptible-infected-recovered (SIR) epidemic model on random graphs following Volz’s original PGF-based derivation (Volz 2008; Miller et al. 2011). This is followed by a rigorous formalisation of the hyperstub configuration model that was first presented by Karrer and Newman (2010). We then demonstrate how this model may be used to generate networks of differing subgraph compositions whilst keeping traditional network metrics such as first and second moments of the degree distribution, clustering and where possible the entire degree distribution, equal. Section 2.3 provides a derivation of the PGF-based approximate ODE model that accurately captures SIR dynamics on hyperstub configuration networks. This derivation is similar to Volz et al.’s (2011) PGF-based extension from configuration and unclustered to clustered networks, but generalised to incorporate arbitrary subgraphs. Finally, Sect. 3 provides an algorithm that automatically generates and solves ODEs presented in Sect. 2.3 for SIR epidemics on networks constructed using a user-specified set of subgraphs.

SIR epidemics on random graphs

The SIR compartmental model involves a population with three types of individuals—susceptible, infected or recovered—whose interactions are modelled by a network. Infection travels across edges at a per-edge rate of and individuals recover, independently, at rate . To account for the heterogeneous contact patterns, the model is centred around the PGF induced by the network’s degree distribution,where is the probability that a randomly chosen node has links. Before we can demonstrate the usefulness of storing the network in this compact way, we need to define the survivor function, . First, we define infectious contact to be the event whereby an infected node transmits to its neighbour , regardless of its state, i.e., irrespective of whether or not it is susceptible (Miller 2011). Next, we select an edge uniformly at random, with nodes and at its ends, and define a direction from node to node . Let be the probability that there has never been infectious contact from node to node by time . Since an infectious contact does not depend on the state of the receiving node, we used the same simplifying assumption as proposed in (Miller 2011), namely, “we disallow infectious contact from node to node ”. Otherwise, may be infected by some other source, and in turn, infect , thus increasing the probability of infectious contact from to . This definition effectively implies that is independent across all edges. For example the probability that a degree two node is susceptible at time is given by , or more generallywhere is the fraction of susceptibles at time . To analytically describe , we need to consider the rate at which a node with degree one becomes infected. This yieldswhere and denote the expected degree of a susceptible node and the expected number of edges per node at time . Hence, denotes the probability that a susceptible and infected node are connected at time . In other words, a node which up to time is susceptible will, on average, become infected at rate . It turns out that can be computed using the PGF and is given bywhich can be interpreted as the expected degree conditional on nodes being susceptible. To compute additional information from the PGF must be extracted, namely the excess degree. This involves selecting an edge at random and following it to its originating node. The observed degree of this node, excluding the edge by which it was selected, is known as the excess degree and has a distribution that is generated byAs before it is possible to condition this on susceptible nodes and thus to compute the expected excess degree of susceptible nodesBy assuming that the expected degree of a newly infected node is equal to the expected degree of a susceptible node, Volz uses the above, multiplied by , to model the expected number of edges the disease can spread across upon infection of a susceptible node. This can be used to derive the equations that describe the flux between edges in different states. Namely, these are given bywhere , and denote the infecting the or the recovering, being created by a node in a edge being infected by an external source to that edge and, finally, the susceptible in a edge being infected by an external source, respectively. Summarising all the above yields the complete system of equations,This concludes the derivation for PGF-based epidemic dynamics on random networks. Volz et al. (2011) extended this methodology to clustered networks by defining a joint probability distribution which describes the typical number of lines and triangles allocated to nodes. This particular derivation has been omitted from this paper. However, in the following section, we will outline a further generalisation of this whereby the joint probability specifies the distribution of subgraphs of various types around nodes. This then leads to more complex PGFs. In Appendix 6.3, we show how the PGF used in the main result of this paper can be made equivalent to the PGF resulting from Volz et al.’s original edge-triangle model.

Hyperstub configuration model

In this paper we generalise the configuration model (Bollobás 1980) to the hyperstub configuration model. Before we specify the model we need to establish how to classify hyperstubs, the set of stubs that connect a node to a subgraph, depending on their parent subgraph and their role within that subgraph.

To generate a hyperstub configuration network model one needs to first decide on a set of subgraphs or building blocks that will form the network. This is then followed by the identification of the number of different hyperstubs induced by the subgraphs: hyperstubs must be uniquely associated with both their parent subgraph and the orbit of their incident nodes  (Karrer and Newman 2010) where the orbit of a node is the set nodes with which it may be permuted such that no edges are created or destroyed. For example, in Fig. 1, subgraph contains two distinct orbits and .

Fig. 1

Subgraph notation and position labelling. Subgraphs are labelled by followed by a symbolic subscript for ease of reference

Once all hyperstubs have been identified it is possible to define a joint probability distribution that specifies the probability of a node having a certain combination of these. For example may denote the probability of a node having and . Using this distribution it is possible to generate hyperstub degree sequences. For network generation these sequences will be subject to cardinality constraints. For example, the sum of the degree sequence of must be divisible by three. Otherwise, the sequence needs to be re-generated. For asymmetric subgraphs, e.g., , the sum of the degree sequences of both types of hyperedge must also be equal. In practice, this can be achieved by generating a suitable degree sequence for one type of hyperedge and then randomly permute it to obtain a second sequence for the second hyperstub. has two degree sequences, one for each hyperstub, and both must be even since we select pairs of nodes from each to form the subgraph.

The network generating algorithm will then form a dynamic list for each hyperstub, where a node with hyperstub degree appears times. This is followed by selecting nodes from the lists, at random and without replacement, and by following the subgraphs’ hyperstub composition in order to construct subgraphs and the network. It is possible that self or multi-edges form in which case the selection is discarded and new samples chosen until a valid selection is obtained. This is repeated until all lists are empty.

In this paper we wish to both computationally generate networks and theoretically describe dynamics on such networks. The PGF of the hyperstub degree distribution provides the link between theory and simulation. The construction of the PGF induced by the hyperstub distribution can be achieved by encoding different levels of detail. At the simplest level nodes may belong to a number of subgraphs without further specifying their orbit or position within the subgraph (Volz et al. 2011). The PGF could be constructed at the level of hyperstubs but would not differentiate between topologically equivalent positions in the subgraph, and this is what we use in our network generating algorithm (nodes may now be allocated asymmetric subgraphs) (Karrer and Newman 2010). Finally, the PGF can be specified by accounting for all details described above with the addition of the precise position of nodes within the subgraph (used in the ODE derivation, Sect. 2.3). For network generation the PGF takes the general form,where is a placeholder and denotes the number of hyperstubs assigned to a node. The symbolic form of the PGF provides more flexibility for computation. Let us consider independently-distributed subgraphs as follows: , and (both hyperstubs of are Poisson distributed with parameter ). The PGF of such a network isFrom this PGF, the average number of subgraphs a node belongs to may be computedBy replacing with , where is the number of stubs contained within the hyperstub , the PGF of the classical degree distribution can be recoveredThe term accounts for the fact that is counted twice, once for each of its hyperstubs. The first and second moments of the degree distribution are directly computed using the linearity of expectation and the fact that . As well as recovering the degree distribution, it is possible to determine the expected number of triangles per node: , since on average each node in is incident to triangles. To summarise, we haveBy including a fourth subgraph in the above example, the equivalent of system Eq. (1) will be underdetermined with 3 equations and 4 unknowns. This allows the first and second moments and the expected number of triangles (and therefore clustering) to be fixed whilst varying the subgraph composition. For example, fixing , and , we can form the underdetermined system,where the columns of the LHS matrix correspond to contributions to , and respectively and denotes a complete subgraph of nodes. From this system it is possible to obtain two valid solutions: (1) and (2) , . Moreover, by replacing with other types of subgraph and updating the L.H.S matrix, several differing network models with the same first and second moments and clustering may be obtained. A selection of such networks used in the results section is listed below:While the three most basic network metrics for the networks above are identical, their degree distributions are not. However, it is also possible to generate classes of networks where the degree distribution is equal between networks but the subgraph composition is not. Let us consider networks constructed purely out of cycles, where, regardless of the length of the cycle, cycle hyperstubs are composed of only pairs of stubs. It is then possible to increase the size of cycles whilst maintaining identical classical degree distributions between different networks. This is implemented in the following way: first, allocate to each node, on average, a pair of cycle hyperstubs, then for each type of network allow the hyperstubs to form increasingly large cycles, starting with then and so on. If the hyperstubs are distributed such that then the classical degree distribution for each network will be such that only even degrees are possible, i.e., denoted for convenience. It is also possible to include a null, random, model for comparison, i.e., a network with degree distribution given by but connected at random. In our investigation we shall be using the following cycle-based networks:where and denote cycles of 5 and 6 nodes (pentagons and hexagons), respectively. Having thus created two classes of networks, the former will be used to show how conventional network metrics may not entirely capture the structure of the network as far as dynamics are concerned; the latter to investigate the effect of cycles of increasing length on dynamics.

SIR epidemics on hyperstub configuration model networks

This section presents the derivation of a general epidemic model for a network built from an arbitrary number of subgraph types. Conceptually, this model uses the node labelling approach of Karrer and Newman (2010) and generalises the PGF-type framework of Volz et al. (2011), Volz (2008). By taking this approach it is possible to derive ODEs that accurately predict the epidemic prevalence on networks that exhibit a variety of exotic subgraphs, both fully- and non-fully connected.

The first step is to choose the set of subgraphs to be included in the network. Let an arbitrary set of subgraphs be labelled by . For example, Fig. 1 shows different subgraphs, which result in distinct node positions, where stands for the total number of nodes over all subgraphs. For clarity, we recall that a hyperstub is the set of half-links connecting a node to a subgraph. This example highlights the key component of the model, namely to distinguish between all nodes of a subgraph even those that are topologically equivalent. This distinction makes it possible to deal with the added complexity of having to account for labelled subgraphs. Each node/position of a subgraph is labelled. This is reflected in a PGF that accounts for each and every node in each and every subgraph. This gives rise to a PGF of the following form:where is a placeholder and is such that is the number of times a node appears in position , .

For each subgraph its state at time is denoted by . This not only describes a subgraph and its state but also the expected number of the given subgraph in the given state at time , i.e., when appended with a state this notation has numerical meaning. Since accounts for the state of node, it will always explicitly depends on . To describe the flux between different subgraph states, infectious events within and between subgraphs need to be considered. This requires a generalisation of which was first given in Sect. 2.1. Accordingly, we now first select a hyperstub at random and then define a direction, from its parent subgraph to its incident node. An infectious contact is now the event that , regardless of its state, becomes infected by one of its adjacent nodes within that subgraph. now needs to reflect a node’s position in the subgraph. Hence, we define to be the probability that the group of edges connecting a node in position to the parent subgraph have not allowed for infectious contact from any infectious node in the subgraph to by time . Again, we impose that cannot transmit infection to the subgraph in question. Under these assumptions, the infectious contact through hyperstubs to position is now independent. A node that appears only times in position remains susceptible with probability . By geometrically compounding all into a PGF, it is possible compute the fraction of the susceptible population. This is given byThis probability is equal to the fraction of susceptible nodes in the population at time (Volz 2008). is referred to as a survivor function. It is dependant on time and may by computed from first principles using the definition of the Poisson process. However, in our formulation, it is computed from variables that denote the expected rate, , at which infection is transmitted to a node in position through the corresponding subgraph. We note that while is commonly used to denote the cumulative probability that infection may occur, we keep it as defined above to be consistent with the current literature on such models (Volz et al. 2011). Each position label has a variable associated with it. The following examples show these rates for positions , and , see Fig. 1:To generate the above identities, we consider a susceptible node in position and list all possible corresponding subgraph states that allow this node to be exposed to infection. can now be used to determine the probability that a susceptible node has an infectious neighbour within a certain subgraph type. This is done by dividing by the number of states that involve a susceptible at position :The expected degree of a susceptible node at position is given bywhere . To compute the expected degree for every position of every subgraph, one can take the Jacobian of evaluated at ,The th entry of this vector evaluated at shall be denoted . A susceptible node in position will have remained susceptible up to time , with probability after which infection may be transmitted at rate . This information may be used to form the following equation: decays at the rate at which a subgraph transmits infection to its node in position , conditional on that node being susceptible.

Once a node is newly infected it is important to determine what, if any, subgraph states are created or destroyed. To do this, we use the susceptible nodes’ excess degree prior to the infection. For the full derivation of susceptibles’ excess degree refer to Appendix 6.1. In this derivation, the excess degree must be generalised to account for the degree of the different positions a node may be in, i.e., . The expected excess degree for susceptible nodes is given bywhere is the Hessian of the PGF. denotes the expected number of positions associated with a node that has been selected at random, but proportionally to the number of positions associated with that node. It is now possible to formulate ODEs describing the evolution subgraph states. We denote the time derivative of a subgraph’s state by . This quantity is dimensionless but not normalised. For example, the number of unique links in a network of size is given by . To form the ODE for the subgraph state , we consider all possible ways in which this state may be created or destroyed, namelywhere denotes the first entry of the vector that is the product of the matrix multiplied from the left by vector . Conceptually denotes the expected number of nodes in position an infection will encounter upon infecting a susceptible node through any possible route, see Fig. 2. The first term on the RHS of Eq. (7) describes this state being destroyed by the infecting the or the recovering. The second term stands for this state being destroyed by the being infected by an outside source. Finally, the last term corresponds to this state being created by the second node of being infected by a source external to the subgraph. To further illustrate this, the equations for and are given,Equations for every state of every subgraph must be derived. In general, we first describe any infection and recovery events of nodes within a subgraph. Next we list all possibilities for susceptible nodes to be infected from sources external to that subgraph using the appropriate terms.

Fig. 2

Graphical representation of . and denote the excess degree of a susceptible node and rate of infection, respectively. We note that newly infected nodes are modelled as previously susceptible nodes so the product is being used to model the expected number of edges infection will be able to spread along upon infecting a susceptible node. This product implicitly considers all possible routes of infection into the node. The left hand side of the figure shows example subgraphs that are the source of infection for the central node. The right hand side of the figure graphically represents the expected excess degree of subgraphs for the central node

To compute network-level prevalences, we recall that can be computed at any time by Eq. (2). is computed directly by differentiating . Namely, since susceptibles become infected and infected nodes recover at rate , we haveThe total number of equations is given by , where denotes the number of nodes in a subgraph. In Appendix 6.2 we give more example ODEs and in Appendix 6.3 we show how our model is equivalent to previous systems developed for complete subgraphs (Volz et al. 2011).

Initial conditions

Let be the fraction of initially infected nodes. Hence, , where is the number of initially infected nodes and is the network size. Initial conditions for the and populations are given byAt time no hyperstub has transmitted infection, therefore, . For a subgraph that contains a single infected node, where is the expected hyperstub degree. For the subgraph with every node susceptible we set . By assuming that only a small fraction of the population, i.e., a single node, is initially infected, we do not allow non-zero initial conditions for subgraphs with more than one infectious node.

Automated code-generation of the mean-field model

We now present our methodology for computationally generating a complete system of equations for a network constructed from subgraphs following a configuration model. This procedure requires the PGF of a hyperstub degree distribution (HDD), the adjacency matrices of corresponding subgraphs, and epidemiological parameters as inputs. The algorithm will output the system of ODEs that will predict the network-level prevalence. Table 1 gives a brief summary of the variables that need to be generated, listed in the order they are generated in this section.

Table 1

Summary of the key system variables and their generation

Variable	Description	Generation
ψ	PGF of the HDD given as a function, not as a series	A symbolic software package can be used to compute the Jacobian and Hessian
θ_i(t)	Survivor functions with their evolution equations given by ODEs	These ODEs can be defined within a single for loop, see Eq. (6)
(S,I,R)	The prevalences of S, I and R, with the latter two given by numerical solutions of ODEs	From Eq. (8), it follows that S=ψ(θ)
T_i	Total rate of infection experienced by an S in position x_i	For a subgraph with m nodes, T_i may be generated by m nested for loops cycling through the possible states that a subgraph can be in, see Eq. (3 4 5)
G_x(S,I,…,R)	Expected prevalence of a subgraph in a given state	The equation for this is computed based on the rate matrix, Z, see Eq. (10)

Let denote the vector of states of a subgraph with denoting a specific state of . For the model, has elements. To generate from , the following steps are needed: (1) cycle through , (2) for each infectious contact to node in state , update to . Using the survivor functions can be computed, see Eq. (6), which are then used to compute the fraction of the population which is susceptible, infected or recovered, see Eq. (8).

The ODEs corresponding to subgraphs need to be represented with a rate matrix, . This matrix encodes all information relating to the given subgraph, namely the excess degrees, rates of infection over subgraphs , epidemiological parameters and , and implicitly encodes the subgraph’s adjacency matrix . To compute , we use Eq. (2) and a symbolic software package to calculate the Jacobian and Hessian of the PGF.

For each subgraph, we initialise the matrix as a square matrix with all entries set to zero. The th column and row of correspond to state . Once populated, the entry contains the rate at which state transitions to state .

To illustrate how to generate , we consider the subgraph, see Fig. 1, with states . We associate the state with the first row and column of . Moving along the top row, when a column index is reached that corresponds to a state that may transition to, we update the entry with the appropriate rate. The first row of is all zero except for and . The second row, corresponding to state , has entries and , see Eq. (7). Fill every row of the matrix in this way, refer to Appendix 6.4 the full matrix corresponding to . The algorithm for this process is given for an arbitrary subgraph in Appendix 6.6, and the corresponding Matlab code is provided as supplemental material but is also available on GitHub at https://github.com/martinritchie/PGF-ODEs

Using the rate matrix, the ODE for the subgraph state yieldsThe final step to generating the full system is to set the initial conditions. Only the initial conditions for subgraph states need computing as and are fixed as per the previous section. This can be done by cycling through each element of . If (a) is a purely susceptible state then we set , and if (b) contains a single infectious individual and is otherwise susceptible, we set . All other states are set to zero, as we assume that with a sufficiently small infectious seed, the probability of having two infectious individuals in a subgraph is zero.

Results

To validate the proposed mean-field model and to assess the goodness of the approximation, we compare results from the ODEs to output from stochastic simulations. Networks were generated following the configuration algorithm, please refer to Appendix 6.5. Typically we generated 500 networks of size and computed a single realisation of the epidemic, according to the Gillespie algorithm with the per link rate of infection and a recovery rate of . Simulations which died out before an outbreak occurred were removed. The simulations were seeded with a single infectious individual and an outbreak was said to occur if 5 % infectious prevalence was achieved. In all plots simulation results and the solution of ODEs are plotted in solid lines and discrete points, respectively.

To start, we test the performance of our model against existing or state of the art models. To do this, in Fig. 3, we show results for two degree distributions that are homogeneous in the classical sense. Their PGFs are given bywhere the variables correspond to subgraphs given in Fig. 1. Figure 3 shows results from a pairwise model with closures at the level of quadruples (House et al. 2009; House 2010). While the classical clustering is easy to compute, the order-four clustering/transitivity ratios were measured following a recently developed subgraph counting algorithm (Ritchie et al. 2014). These are defined as the ratio of a given subgraph count to all open and closed paths of length four, both counted uniquely. Currently, this model operates using an average or homogenous degree and stores no information about the degree distribution, but does assume random mixing of subgraphs.

Fig. 3

Performance of other models. Lines, circles and squares correspond to simulation average, ODE solution and pairwise ODE solution, respectively. All networks are homogeneous with . The lower peaks correspond to networks generated with each node allocated one of each corner type of a with clustering . Data with higher peak correspond to networks generated with a single and two subgraphs yielding

All models perform well in capturing the epidemic dynamics on networks generated using the PGF given by , see Fig. 3 with higher epidemic peak. However, when networks are created using the PGF given by , see lower peak in Fig. 3, the pairwise model struggles to accurately capture the dynamics, both anticipating and compressing the epidemic’s time scale or duration and underestimating the final size (data not shown). The pairwise model does not encode any information relating to degree or subgraph distribution and hence a homogeneous random set-up, as used here, would be an appropriate choice.

The key advantage of our algorithm over existing ones is that it can handle non-fully connected subgraphs. To test this, in Fig. 4, we utilise networks models C1–C4 as specified in Sect. 2.2. Figure 4 shows plots of simulation average compared to the ODE’s solution for the four network types. We observe that the epidemic behaviour of networks composed of increasingly large cycles quickly converge to that of the random null case. It has previously been observed that for networks with the same degree distribution, an increasing level of clustering slows the epidemic transmission and requires a higher transmission rate in order to observe a successfully spreading epidemic (Keeling 1999b; Green and Kiss 2010). This occurs for two reasons: (1) subgraphs that are densely connected share fewer connections to the rest of the network so an initial seed will be restricted to one part of the network and (2) this same effect leads to infectious nodes competing for susceptible nodes. While this may make transmission more efficient locally, it does limit further seeding in fully susceptible parts of the network. Figure 4 shows that the effect of is similar to that of the clustered network, but less pronounced; both the time and size of peak infectious prevalence is delayed and reduced when compared to the null case. For cycles larger than four nodes this behaviour is less pronounced and the epidemics for larger cycles converge to the null case, as observed with .

Fig. 4

Clustering and cycles. Solid lines and markers correspond to simulation average and ODE solution, respectively. From darkest to lightest, the solid lines correspond to: , , and , i.e., each network used has an identical degree distribution given by . Clustering is and for the and other networks, respectively. For clarity, ODE solutions for only the two extreme cases, the null and triangle network, have been included. Note that the output from the network composed of is close to that of the null case. Epidemics corresponding to cycles of length six have been computed but omitted due to their similarity to the null case. Only two ODE solutions have been included for upper and lower cases

To highlight the flexibility of our model and its wide-ranging applicability to systematically investigating the impact of higher-order network structure, in Fig. 5, we consider four networks with the same first and second moments, and the same classical clustering but generated using different families of subgraphs, see models 1–4 Sect. 2.2. Figure 5 shows simulation average for all four networks and the solution of ODEs for the upper and lower cases, models 1 and 4 respectively.

Fig. 5

Clustering via differing subgraphs. Solid lines and markers correspond to simulation average and ODE solution, respectively. From darkest to lightest the solid lines correspond to: ; , ; , and , , where and denote complete pentagon and hexagon subgraphs, respectively. The networks were generated so that , and . The downward trend of peak prevalence corresponds to networks composed of complete subgraphs of increasing size. The larger subgraphs lead to more connections within the group rather than to the rest of the network

Figure 5 shows a clear trend whereby larger subgraphs lead to epidemics with smaller peak prevalence. A second more subtle trend shows a delay in time until peak prevalence. Subgraphs of larger size lead to a significant difference in the behaviour of epidemics and echo what was observed for increasing levels in clustering. This could be explained by considering a subgraph with average degree . When the network will exhibit extreme clustering, where isolated structures are increasingly densely connected at the cost of becoming isolated. This effect is more subtle than clustering but it can be significant. This suggests that the accuracy of future models would improve if they can correctly account for networks’ subgraph composition, particularly subgraphs beyond that of triangles.

Finally, the data in Fig. 5 has been produced using networks that do not have the same degree distribution but do have equal first and second moments, and clustering. To better understand how the non-equal higher moments may have affected the results, we have simulated epidemics on the corresponding random networks, Model and Model , see Appendix 6.7. This plot shows that the differences observed in Fig. 5 cannot be explained by the difference in the degree distribution alone. Thus, generating identical clustering but using different subgraphs can lead to non-negligible differences in epidemic dynamics. This behaviour has been previously observed, see (Green and Kiss 2010; Ritchie et al. 2014), and warrants a more detailed investigation that will be the subject of a follow-up publication.

Discussion

Higher-order structures, captured for example as different subgraph compositions and arrangements in a network, have been identified as features of real networks. Examples include households, social interactions and biological networks. These building blocks of networks have been shown to play a key role in defining a network’s topology and can have significant impact on the functions of the network or on the dynamical processes unfolding on the network. Despite this, the modelling toolset in this direction is underdeveloped. Here, we provided an approach that considerably extends the scope of the current modelling framework by enabling us to consider arbitrary sets of exotic subgraphs as building blocks for the network. Our approach also offers control over the arrangements of subgraphs and, more importantly and uniquely, an automated way of generating a system of ODEs that accurately capture the prevalence profile for a wide range of subgraph sets, as shown in the results section.

The previous section has shown how higher-order structures may be investigated using this model. Moreover, we provided the first example of generating classes of networks constructed using different subgraph sets while keeping degree, variance and clustering constant. For example, we showed that epidemics on networks with no clustering, but exhibiting open loops, display features which are significantly different to those observed in classical random networks with effectively no clustering. Equally, we have shown that different subgraph combinations or arrangements can create higher-order structure that may significantly affect the epidemic dynamics. Our work opens the possibility to carry out a wide-ranging and systematic investigation of the impact of subgraphs and higher-order structure on dynamics on networks. When presented with real world network data whose structure can be explained by a set of subgraphs, all that will be needed in order to apply our framework is to extract the subgraphs and their distribution around nodes. A possible limitation to the widest applicability is the number of nodes in the largest subgraph. However, as shown by our results when going from squares to pentagons, it is likely that the effect of higher-order structures will decay, or be less marked, as their size increases.

There are two key ways in which this work may be extended: (a) generalisation to dynamics. Due to the definition of it is currently not possible to apply this model to dynamics. However, all the framework relating to network structure is independent from this variable and may therefore still be appropriate. (b) The subgraph approach is highly suitable for adaptation to household models. Household models typically specify a distribution of household sizes overlaid on a contact network to produce a well-connected network (House and Keeling 2009; Ball and Sirl 2012). A successful incorporation of such network in our framework could lead to a highly relevant set of household models.

Supplementary material 1 (zip 20 KB)