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. Electronic supplementary material The online version of this article (doi:10.1007/s00285-015-0884-1) contains supplementary material, which is available to authorized users.


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).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 wellperforming 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 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 nonfully connected subgraphs, such as overlapping triangles.Using such networks, Volz et al. (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 & Newman (Karrer and Newman, 2010) have provided final epidemic size results for networks built by using different mixtures of subgraphs.Furthermore, House and colleagues (House et al., 2009;House, 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 & Newman's approach, which combines a wider variety of subgraphs, can only characterise large-time limits.Finally, to our knowledge, House et al.'s approach has not been compared to stochastic simulations and it will perform poorly for heterogeneous networks (House et al., 2009).
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-fullyconnected 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-infectedrecovered (SIR) epidemic model on random graphs following Volz's original PGF-based derivation (Volz, 2008).This is followed by a rigorous formalisation of the hyperstub configuration model that was first presented by Karrer & Newman (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.Sec.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 PGFbased extension from configuration and unclustered to clustered networks (Volz et al., 2011), but generalised to incorporate arbitrary subgraphs.Finally, sec. 3 provides an algorithm that automatically generates and solves ODEs presented in sec.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 p(k) is the probability that a randomly chosen node has k links.Before we can demonstrate the usefulness of storing the network in this compact way, we need to define the survivor function, θ(t).First, we define infectious contact to be the event whereby an infected node v transmits to its neighbour u, 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 u and v at its ends, and define a direction from node v to node u.Let θ(t) be the probability that there has never been infectious contact from node v to node u by time t.Since an infectious contact does not depend on the state of the receiving node, Volz proposed the following simplifying assumption, "we disallow infectious contact from node u to node v".Otherwise, u may be infected by some other source, and in turn, infect v, thus increasing the probability of infectious contact from v to u.This definition effectively implies that θ(t) is independent across all edges.For example the probability that a degree two node is susceptible at time t is given by θ(t) 2 , or more generally where S(t) is the fraction of susceptibles at time t.To analytically describe θ(t), we need to consider the rate at which a node with degree one becomes infected.This yields where M S (t) and M SI (t) denote the expected degree of a susceptible node and the expected number of SI edges per node at time t.Hence, M SI (t)/M S (t) denotes the probability that a susceptible and infected node are connected at time t.In other words, a node which up to time t is susceptible will, on average, become infected at rate τ M SI (t)/M S (t).It turns out that M S (t) can be computed using the PGF and is given by which can be interpreted as the expected degree conditional on nodes being susceptible.To compute M SI (t) 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 by As before it is possible to condition this on susceptible nodes and thus to compute the expected excess degree of susceptible nodes By 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 by where M SI (t)(τ + γ), 2δ S (t)M SS (t) and δ S (t)M SI (t) denote the I infecting the S or the I recovering, M SI being created by a node in a SS edge being infected by an external source to that SS edge and, finally, the susceptible in a SI 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. extended this methodology to clustered networks by defining a joint probability distribution which describes the typical number of lines and triangles allocated to nodes (Volz et al., 2011).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 App.7.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 G contains two distinct orbits {x 14 , x 17 } and {x 15 , x 16 }.
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 f (x, y) = p x,y may denote the probability of a node having x × G 0 and y × G .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 G must be divisible by three.Otherwise, the sequence needs to be re-generated.For asymmetric subgraphs, e.g., 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.G 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 k i appears k i 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, sec.2.3).For network generation the PGF takes the general form, where ẑ = (z 1 , z 2 , . . ., z m ) is a placeholder and ĥ = (h 1 , h 2 , . . ., h m ) denotes the number of h i hyperstubs assigned to a node.The symbolic form of the PGF provides more flexibility for computation.Let us consider subgraphs distributed as follows: G 0 ∼ P ois(λ 1 ), G ∼ P ois(λ 2 ) and G ∼ P ois(λ 3 ) (both hyperstubs of G are Poisson distributed with parameter λ 3 ).The PGF of such a network is From this PGF, the average number of subgraphs a node belongs to may be computed By replacing z i with z a , where a is the number of stubs contained within the hyperstub h i , the PGF of the classical degree distribution can be recovered The z 5/2 term accounts for the fact that G 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 V ar(aX) = a 2 X.As well as recovering the degree distribution, it is possible to determine the expected number of triangles per node: = λ 2 + 3/2λ 3 , since on average each node in G is incident to 3/2 triangles.To summarise, we have By 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 k = 4, V ar(k) = 8 and = 2, we can form the underdetermined system, where the columns of the LHS matrix correspond to contributions to k , V ar(k) and respectively and G ic denotes a complete subgraph of i nodes.From this system it is possible to obtain two valid solutions: (1) G ∼ P ois( 2) and (2) G 0 ∼ P ois(9/2), G 6c ∼ P ois(3/10).Moreover, by replacing G 6c 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 G then G and so on.If the hyperstubs are distributed such that h i ∼ P ois(2) then the classical degree distribution for each network will be such that only even degrees are possible, i.e., P (degree = 2k) = P (degree = k|P ois(2)) denoted G 0 ∼ 2P ois(2) for convenience.It is also possible to include a null, random, model for comparison, i.e., a network with degree distribution given by G 0 ∼ 2P ois(2) but connected at random.In our investigation we shall be using the following cycle-based networks: where G and G 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 SIR 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. (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 {G 1 , G 2 , . . ., G M }.For example, Fig. 1 shows M = 5 different subgraphs, which result in m = 17 distinct node positions, where m 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 α = (α 1 , α 2 , . . ., α m ) is a placeholder and ŷ = (y 1 , y 2 , . . ., y m ) is such that y i is the number of times a node appears in position x i , i = 1, . . ., m.
For each subgraph its state at time t is denoted by G x (S, I, . . ., R).This not only describes a subgraph and its state but also the expected number of the given subgraph in the given state at time t, i.e., when appended with a state this notation has numerical meaning.Since G x (S, I, . . ., R) accounts for the state of node, it will always explicitly depends on t.To describe the flux between different subgraph states, infectious events within and between subgraphs need to be considered.This requires a generalisation of θ(t) which was first given in sec.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 u, regardless of its state, becomes infected by one of its adjacent nodes within that subgraph.θ(t) now needs to reflect a node's position in the subgraph.Hence, we define define θ i (t) to be the probability that the group of edges connecting a node u in position x i to the parent subgraph have not allowed for infectious contact from any infectious node in the subgraph to u by time t.Again, we impose that u cannot transmit infection to the subgraph in question.Under these assumptions, the infectious contact through hyperstubs to position x i is now independent.A node that appears only k times in position x i remains susceptible with probability θ k i (t).By geometrically compounding all θ i (t) into a PGF, it is possible compute the fraction of the susceptible population.This is given by This probability is equal to the fraction of susceptible nodes in the population at time t (Volz, 2008).θ(t) 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, T i , at which infection is transmitted to a node in position x i through the corresponding subgraph.We note that while T 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 x i has a T i variable associated with it.
The following examples show these rates for positions x 1 , x 2 and x 3 , see Fig. 1: (3) (4) To generate the above identities, we consider a susceptible node in position x i and list all possible corresponding subgraph states that allow this node to be exposed to infection.T = (T 1 , T 2 , . . ., T m ) 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 T i τ −1 by the number of states that involve a susceptible at position x i : A, B, C, D) .
The expected degree of a susceptible node at position x i is given by where θ = (θ 1 , θ 2 , ..., θ m ).To compute the expected degree for every position of every subgraph, one can take the Jacobian of ψ evaluated at x = θ, The i th entry of this vector evaluated at α = θ shall be denoted J i .A susceptible node in position x i will have remained susceptible up to time t, with probability θ i after which infection may be transmitted at rate T i /J i .This information may be used to form the following equation: θi (t) decays at the rate at which a subgraph transmits infection to its node in position x i , 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's excess degree prior to the infection.For the full derivation of susceptibles' excess degree refer to App. 7.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., k i , i = 1, 2, . . ., m.The expected excess degree for susceptible nodes is given by where H(ψ) is the Hessian of the PGF.∆ i,j denotes the expected number of x j positions associated with a node that has been selected at random, but proportionally to the number of x i 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 (SI) links in a network of size N is given by [SI] = N G 0 (SI).To form the ODE for the subgraph state G 0 (SI), we consider all possible ways in which this state may be created or destroyed, namely where (T ∆) 1 denotes the first entry of the vector that is the product of the matrix ∆ multiplied from the left by vector T .Conceptually (T ∆) i denotes the expected number of nodes in position x i 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 I infecting the S or the I recovering.The second term stands for this state being destroyed by the S being infected by an outside source.Finally, the last term corresponds to this state being created by the second node of G 0 (SS) being infected by a source external to the subgraph.To further illustrate this, the equations for G 0 (SS) and G 0 (IS) 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 (T ∆) terms.
To compute network-level prevalences, we recall that S(t) can be computed at any time by Eq. (2).İ(t) is computed directly by differentiating S(t).Namely, since susceptibles become infected and and infected nodes recover at rate γ, we have The total number of equations is given by 2 + m + M i=1 3 |G i | , where | • | denotes the number of nodes in a subgraph.In App.7.2 we give or more example ODEs and in App.7.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, = I 0 /N , where I 0 is the number of initially infected nodes and N is the network size.Initial conditions for the I and R populations are given by At time t = 0 no hyperstub has transmitted infection, therefore, θ i (t = 0) = 1.For a subgraph that contains a single infected node, G(t = 0) = k i where k i is the expected hyperstub degree.For the subgraph with every node susceptible we set G(t = 0) = (1 − ) k i .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 ?? gives a brief summary of the variables that need to be generated, listed in the order they are generated in this section.
Let G denote the vector of states of a subgraph G with G i denoting a specific state of G.For the SIR model, G has 3 |G| elements.To generate T i from G, the following steps are needed: (1) cycle through G, (2) for each infectious contact to node i in state G j , update T i to T i = T i + G. Using T 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, Z.This matrix encodes all information relating to the given subgraph, namely the excess degrees, rates of infection over subgraphs T , epidemiological parameters τ and γ, and implicitly encodes the subgraph's adjacency matrix g.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 Z as a square matrix with all entries set to zero.The i th column and row of Z correspond to state G i .Once populated, the entry Z i,j contains the rate at which state i transitions to state j.
To illustrate how to generate Z, we consider the G 0 subgraph, see Fig. 1, with states G = (SS, SI, SR, IS, II, IR, RS, RI, RR).We associate the state G 0 (SS) with the first row and column of Z. Moving along the top row, when a column index is reached that corresponds to a state that G 0 (SS) may transition to, we update the entry with the appropriate rate.The first row of Z is all zero except for Z 1,2 = (T ∆) 2 and Z 1,4 = (T ∆) 1 .The second row, corresponding to state G 0 (SI), has entries Z 2,3 = γ and Z 2,5 = τ + (T ∆) 1 , see Eq. ( 7).Fill every row of the matrix Z in this way, refer to the SI text for the full matrix corresponding to G 0 .The algorithm for this process is given for an arbitrary subgraph in App.7.6, and the corresponding Matlab code is provided as supplemental material.
Using the rate matrix, the ODE for the subgraph state The final step to generating the full system is to set the initial conditions.Only the initial conditions for subgraph states need computing as I(0), R(0) and θ i (0) are fixed as per the previous section.This can be done by cycling through each element of G.If (a) G i is a purely susceptible state then we set G i 0 = J i (1 − ), and if (b) G i contains a single infectious individual and is otherwise susceptible, we set G i 0 = J i .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 App. 7.5.Typically we generated 500 networks of size N = 15000 and computed a single realisation of the epidemic, according to the Gillespie algorithm with the per link rate of infection τ = 1 and a recovery rate of γ = 1.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 by where the variables α i 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.
All models perform well in capturing the epidemic dynamics on networks generated using the PGF given by ψ 1 , see Fig. (3) with higher epidemic peak.However, when networks are created using the PGF given by ψ 2 , 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, but predicting correctly the peak prevalence.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 sec.2.2.Fig. 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.Fig. 4 shows that the effect of G 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 G .
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 sec.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.
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 k s .When k < k s 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 1 : G 0 ∼ 2P ois(2) and Model 4 : G 0 ∼ P ois(3) + 5P ois(1/5), see App. 7.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.

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, all in the classic sense, fixed.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 SIS dynamics.Due to the definition of θ(t) it is currently not possible to apply this model to SIS 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 et al., 2012).A successful incorporation of such network in our framework could lead to a highly relevant set of household models.
6 Figures and tables x 9 x 6 x 7 x 8 x 10 x 13 x 11 x 12 x 14 x 16 x 15 x 17 Figure 1: Subgraph notation and position labelling.Subgraphs are denoted by G followed by a symbolic subscript for ease of reference.
Figure 2: Graphical representation of (T ∆) i .∆ and T 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 (T ∆) i is being used to model the expected number of x i 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 red node.The right hand side of the figure graphically represents the expected excess degree of G subgraphs for the red node.G ∼ P ois(2); G 0 ∼ P ois(2), G ∼ P ois(2/3); G 0 ∼ P ois(8/3), G cp , ∼ P ois(1/3) and G 0 ∼ P ois(3), G ch , ∼ P ois(1/5), where cp and ch denote complete pentagon and hexagon subgraphs, respectively.The networks were generated so that k = 4, var(k) = 8 and φ = 0.2.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.
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).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).
Table 1: Summary of the key system variables and their generation.

Appendix
In this Appendix we (a) give a more detailed explanation of the excess degree, (b) provide ODEs for an example network, (c) show how our generalised model reduces to a previous model under certain conditions, (d) provide an example state transition matrix, (e) give pseudocode for both the subgraph-based configuration model and the algorithm used to generate the state transition matrix and, finally, (f) compare epidemic dynamics on two configuration model networks with their degree distributions being different but with the same mean and variance.

Excess degree
Recall the probability generating function (PGF) of a network's hyperstub degree distribution with m nodal positions: where α = (α 1 , α 2 , . . ., α m ) is a placeholder, and ŷ = (y 1 , y 2 , . . ., y m ) is such that y i denotes the number of times a node appears in position x i , i = 1, 2, . . ., m.The PGF of the excess degree distribution is a critical component in our derivation and it is illustrative to see how it is computed.To see this, we first compute the expected excess degree, select a node at random but proportional to its number of x i hyperstubs, y i p ŷ. Next, to obtain the expected x i degree, sum this product over all nodes The above sum considers each and every node from which an x i hyperstub originates.However, in hyperstub configuration model networks there is usually more than one type of hyperstub and this adds an additional level of detail to the excess degree.The excess degree now may incorporate two different hyperstubs into its calculations.It is now possible to describe a nodes x i degree but conditional on it being selected through ones of its x j hyperstubs.More formally we can compute the expected excess degree using conditional expectation, E(x j |x i = y i ), which yields where δ x j ,x i denotes the expected x i hyperstub degree observed from a node selected proportionally to its x j hyperstub degree.The denominator is given by Eq. ( 11), and the numerator is specified by

ODEs for an example network
The following provides ODEs for a simple example network composed of only G 0 and G .When deriving ODEs by hand listing out equations for T i is a good starting point as they include many of the subgraph states, i.e.G 0 (SI), and can be used as the start of a check list when listing state equations It is important to node the above equations will not list every subgraph state and that a subgraph composed of n will have 3 n state equations.For example, the first few state equation for G 0 are given by Ġ0 with equations for the following being omitted { Ġ0 (SR), Ġ0 (II), Ġ0 (IR), Ġ0 (RS), Ġ0 (RI), Ġ0 (RR)}, Similarly, sample ODEs for the G subgraph, taken from a system of 27 ODEs, are: with equations for the following being omitted { Ġ0 (SSR), Ġ0 (SII), Ġ0 (SIR), Ġ0 (SRS), Ġ0 (SRI), Ġ0 (SRR), Ġ0 (ISI), Ġ0 (ISR), Ġ0 (IIS), Ġ0 (III), Ġ0 (IIR), Ġ0 (IRS), Ġ0 (IRI), Ġ0 (IRR), Ġ0 (RSS), Ġ0 (RSI), Ġ0 (RSR), Ġ0 (RIS), Ġ0 (RII), Ġ0 (RIR), Ġ0 (RRS), Ġ0 (RRI), Ġ0 (RRR)}, Each hyperstub will have a survivor function and a corresponding ODE describing its evolution, as follows: The fraction of the population that is susceptible or infected is computed by compounding θ i into the PGF.Symbolically, this is computed by the following where ψ is the generating function that generates the hyperstub degree distribution and θ = (θ 1 , θ 2 , θ 3 , θ 4 , θ 5 ) is the probability that infection via subgraphs of types one to five has not been transmitted.The total system size for this example network is given by with each term in the above corresponding to G 0 , G , survivor functions and epidemic prevalence, respectively.In general, the total number of equations is given by: where G i denotes a subgraph, |G i | is the number of nodes in a subgraph, and m is the total number of subgraphs.

Equivalence to previous model for complete subgraphs
The PGF formulation originally proposed by Volz et al. (Volz et al., 2011) is equivalent to our proposed model in the case of complete subgraphs.Consider an arbitrary complete subgraph composed of l nodes and a network that is composed only of this subgraph.If positions within the subgraph are labelled distinctly, {x 1 , x 2 , . . .x l }, as we have done in our approach, then the PGF of such a network is given by where ŷ = (y 1 , y 2 , . . ., y l ).Volz et al.'s framework treats all topologically equivalent positions as one single position.Thus, in this case, the subgraph has a single label, x, that corresponds to a single count, y, and the PGF takes the following form We now show how one may obtain Eq. ( 15) from Eq. ( 14).Since both PGFs describe the same network, the rate at which our formulation allocates position x i must be 1/l the rate at which Volz et al.'s formulation allocates x.If we replace the unique position labels of Eq. ( 14) with a single position marker (such as in Volz et al.'s model), the following expression is obtained where the following substitutions, y i = y/l and α i = α, were made so that α y is the result of the above product.Now, every time an x i is allocated, we allocate an x instead.Finally, since p ŷ is a joint distribution of l identically distributed independent random variables, i.e., ŷ = (y/l, y/l, . . ., y/l), we get: It is also possible to translate between the two models elsewhere in the derivation.As an example, in our approach, infection over lines is given by T 1 and T 2 , as per Eq.(3).By summing these values, the equivalent values used in Volz et al.'s formulation may be recovered.Following our derivation, first let G 0 (SI) ≡ G 0 (IS) and: Since each G 0 is generated from a PGF that allocates positions at rate 1/2 that of Volz et al.'s PGF, the 2 will cancel yielding τ G 0 (SI).However, it is only necessary to show equivalence between the two PGFs since all other variables follow from this.

State transition matrix
The state transition matrix for G 0 (lines) is given by:

Algorithm 1 -Hyperstub CM algorithm
Algorithm 1: The hyperstub configuration model.In this implementation, multiple-edges are over written (line 35) but self-edges are permitted.To prevent this, if nodes already share an edge or a node has been selected twice (self-edge) lines 25-28 are repeated until a valid selection is made.This reselection step has been omitted below for readability.for For each hyperstub of g i do

23
% Select unfiformly at random and without

24
% replacement a node incident to each desired hyperstub.

Figure 3 :Figure 4 :Figure 5 :
Figure3: Performance of other models.Lines, circles and squares correspond to simulation average, ODE solution and pairwise ODE solution, respectively.All networks are homogeneous with k = 5.The lower peaks correspond to networks generated with each node allocated one of each corner type of a G with clustering φ = 0.3.Data with higher peak correspond to networks generated with a single G 0 and two G subgraphs yielding φ ≈ 0.

% 7 H:
Each row of K corresponds to single node's hyperstub sequence.6 K: the hyperstub degree sequence, a non-square matrix K ∈ N N ×H 0 , the number of hyperstub types, 8 A: the adjacency matrix of the network, A ∈ {0, 1}

25 n 1 %
= rand-sample(h i 1 ) 26 n 2 = rand-sample(h i 2 )The following compares pairs of the selected nodes 31 % to determine their connectivity in A.

Figure 6 :
Figure6: The effect of higher moments.The solid and discrete plots correspond to the null networks G 0 ∼ 2P ois(2) and G 0 ∼ P ois(3) + 5P ois(1/5) respectively, i.e. the null cases for the triangle and hexagon networks.Both plots have equal first and second moments and clustering equal to that of a random network.The difference observed is a result of non-equal higher moments and is not enough to explain the difference observed in Fig.5.
: a dynamic list of nodes that are incident to H i hyperstubs, 12 g i : the adjacency matrix of a subgraph, g ∈ {0, 1} n i ×n i , N ×N , 9 M : the number of subgraphs, 10 H i : the degree of a specific hyperstub, 11 h i 13 n i : the number of nodes in g i .14 Procedure 15 % The following creates dynamic lists the, 'hyperstub bins'.16 for every node i do 17 for each H j do 18 append K i,j multiples of node(i) to the hyperstub bin(h j ) 19 end 20 end 21 for For each subgraph g i do 22