# Mixing patterns in graphs with higher-order structure

Peter Mann<EMAIL_ADDRESS>Lei Fang  Simon Dobson  School of Computer
Science, University of St Andrews, St Andrews, Fife KY16 9SX, United Kingdom

###### Abstract

In this paper we examine the percolation properties of higher-order networks
that have non-trivial clustering and subgraph-based assortative mixing (the
tendency of vertices to connect to other vertices based on subgraph joint
degree). Our analytical method is based on generating functions. We also
propose a Monte Carlo graph generation algorithm to draw random networks from
the ensemble of graphs with fixed statistics. We use our model to understand
the effect that network microstructure has, through the arrangement of
clustering, on the global properties. Finally, we use an edge disjoint clique
cover to represent empirical networks using our formulation, finding the
resultant model offers a significant improvement over edge-based theory.

###### pacs:

Valid PACS appear here

††preprint: Draft

##  I Introduction

Traditional network theory is limited to pairwise interactions. Higher-order,
non-dyadic interactions lead to fresh insight into collective dynamics that
cannot be described by edge-based theory PhysRevLett.124.218301 ;
iacopini_petri_barrat_latora_2019 ; giusti_ghrist_bassett_2016 ; XXXXXXXJJJJJJ
. Theoretical models that include higher-order interactions become
increasingly complicated compared to edge-based models and generalisations of
the definitions of traditional concepts, such as assortativity or clustering
are often required.

The ability to generate a random graph from an ensemble of random graphs that
have statistically equivalent properties allows us to model, and subsequently
to understand, the topological features that drive the behaviour of dynamical
processes that occur over networks. Topological features can be extrinsic
(finite-size effects), whilst others are intrinsic, such as the clustering
coefficient, $\mathcal{C}$ PhysRevE.68.026121 ; PhysRevLett.103.058701 , or
the degree assortativity PhysRevE.67.026126 ; PhysRevLett.89.208701 ;
balogh_palla_kryven_2020 . Some properties are locally defined, and as such
are subject to local anisotropy; whilst others characterise global attributes
of the network that are averaged over all vertices. Finally, ensemble
properties manifest for a collection of graphs that are statistically
equivalent in some manner. Often, the intrinsic properties of networks are not
independent of one another. For instance, two networks with equivalent
clustering coefficients could exhibit different responses to a given dynamical
process due to the organisation of their triangles. Clustering is known to
increase the critical point of the bond percolation process and reduce the
size of the expected giant connected component (GCC) PhysRevE.80.020901 ;
gleeson_2009 ; PhysRevE.83.056107 ; PhysRevE.81.066114 ; PhysRevE.104.024303 ;
PhysRevE.103.012313 ; PhysRevE.101.062310 ; PhysRevE.103.012309 ;
PhysRevE.104.024304 . However, vertices involved in many triangles that tend
to connect to other vertices with a high triangle count may induce assortative
correlations among the overall degrees. In turn, the critical point of
positively assorted networks is reduced compared to the neutral model
PhysRevE.67.026126  and other properties can also be affected
PhysRevE.67.027101 . The observed behaviour of a percolation process depends
on the interplay between the clustering and the assortativity
PhysRevLett.103.058701 . Therefore, it is important to model both the role of
clustering as well as the role of degree assortativity.

One modelling framework is the configuration model (CM) that enables the
construction of a particular locally tree-like random graph from a
distribution of degrees newman_strogatz_watts_2001 . CM networks are regarded
to be infinitely large, such that fluctuations from ensemble averages do not
effect the statistics of the model. This model has received significant
attention in the literature, and has been extended to the so-called
generalised configuration model (GCM) to incorporate non-trivial clustering
and higher-order structure among vertices PhysRevE.82.066118 . In the GCM, a
set $\vec{\tau}$ of substrate motifs is defined (such as cliques or chordless
cycles of different sizes) and the number of edge-disjoint motifs of each
topology that a given vertex belongs to is described by a tuple, known as its
joint degree that takes values in $\mathbb{Z}_{\geq 0}^{n}$ where $n$ is the
number of distinct motif topologies in the model. For example, for a model
containing only 2- and 3-cliques, a vertex that belongs to 3 ordinary edges, 2
triangles has joint degree $(s,t)=(3,2)$, see Fig 1. Vertex joint degrees are
drawn from a multivariate distribution of joint degrees, $p_{s,t,h,\dots}$,
that defines the probability that a vertex belongs to $s$ ordinary edges, $t$
triangles, $h$ 4-cycles etc. When $\vec{\tau}$ is chosen to contain only
ordinary edges (2-cliques), the GCM collapses to the standard CM.

Figure 1: A snapshot of a 2- and 3-clique clustered graph with vertices
labelled with their joint degree tuples $(s,t)$.

The CM is known to create networks with no assortativity among the degrees of
neighbouring vertices due to the stochastic nature of the construction
algorithm PhysRevE.82.066118 ; vertex degrees are independent random
variables. Assortativity can be inserted, in a controlled manner, into the CM
via a degree-preserving Markov Chain Monte Carlo (MCMC) rewiring algorithm
where the correlations converge to a target value PhysRevLett.89.208701 ;
PhysRevE.67.026126 .

Hasegawa et al PhysRevE.101.062310  examined the role of clustering and
assortativity in the percolating component of tree-triangle networks.
Recently, Mann et al investigated the inter-subgraph correlation properties
that arise naturally in the percolating component of GCM networks that are
composed of arbitrary sized cliques that have neutral mixing patterns in the
substrate graph PhysRevE.105.044314 .

In this paper, we extend that work to examine the inter-subgraph properties of
GCM networks that have arbitrary clique clustering and subgraph correlations
in the substrate network. To do this, we develop an analytical method based on
generating functions and introduce an MCMC rewiring algorithm to simulate such
networks. These steps allow the synthesis of random graphs drawn from an
ensemble of graphs with a given clustering coefficient $\mathcal{C}$ and
correlation structure among the joint degree tuples. This model can be used to
probe the role that both clustering and assortativity have on the properties
of the bond percolation process and will likely lead to the detailed study of
higher-order structure for other kinds of dynamical processes on networks
beyond the percolation model. Additionally, our model allows the clustering
and correlation structure from an empirical network to be gathered through
means of an edge-disjoint clique cover burgio_arenas_gomez_matamalas_2021 ;
PhysRevE.105.044314 . Through our MCMC algorithm, synthetic networks, possibly
containing additional vertices to the empirical network, can then be created.
In this way, our model allows ensemble representations with fixed statistics
of empirical networks to be created and studied. We remark, that the
information in our model can readily be compressed from inter-subgraph
correlations to correlations between overall degrees; however, there is a
corresponding information loss following this process.

##  II Theoretical

Correlations can arise via three distinct processes in GCM networks. Firstly,
correlation between the joint degrees of neighbouring vertices that belong to
the same clique type. For instance, three vertices in a triangle must all have
$t\neq 0$. In this manner, vertices with membership of a given 3-clique are
somewhat assorted together; they will never connect to vertices that don’t
belong to triangles. However, this correlation does not extend beyond the
scope of the motif to which the edge belongs and so, we refer to it as the
trivial correlation. In other words, the excess joint degrees (the remaining
joint degree tuples excluding the motif to which the edge belongs) of two
vertices at the ends of a randomly selected edge are not correlated in the
default GCM construction. And, excluding the trivial assortativity due to
belonging to the same motif, there are no specific correlations between
vertices at play.

We can also discuss assortativity between the overall degrees of vertices in
the GCM. The overall degree of a vertex that belongs to $s$ 2-cliques, $t$
3-cliques, $w$ 4-cliques (and so on) is given by $k=s+2t+3c+\dots$. This
assortativity often arises in an uncontrolled manner. For example, a GCM
network composed of 2-cliques and 10-cliques may be absent of any correlation
structure among the distribution in the numbers of cliques a vertex belongs
to; however, by virtue of the increase in overall degree, vertices that belong
to 10-cliques will likely be of higher-degree than vertices that belong to
just 2-cliques, and as a result, will exhibit positive correlations among the
overall degrees if this is the case PhysRevE.80.020901 ; PhysRevE.101.062310 .

Finally, we can inject correlation into GCM networks deliberately, in a
detailed and controlled fashion, by specifying how likely it is that a vertex
with a given joint degree tuple $(s_{0},t_{0},\dots,c_{0})$ has a neighbour
with a certain joint degree tuple $(s^{\prime},t^{\prime},\dots,c^{\prime})$.
In general, this probability is different depending on the topology of the
edge (to which clique it belongs to). This kind of correlation is the most
general of the three types that arise in the GCM; since the other correlation
types are subsets of this; and so, it is the focus of this paper.

We restrict our attention to GCM networks that are composed of 2- and
3-cliques, symbolised by $\vec{\tau}=\\{\bot,\Delta\\}$, respectively. Let
$e_{\tau,s_{0}t_{0},s^{\prime}t^{\prime}}$ be the probability that a randomly
chosen $\tau$-edge leads to vertices with excess (remaining) joint degrees of
$(s_{0},t_{0})$ and $(s^{\prime},t^{\prime})$ where $s$ and $t$ are the
numbers of 2- and 3-cliques, respectively, that each vertex belongs to, see
Fig 2. The joint excess joint degree distribution is the data input to the
model; other familiar quantities such as excess joint degree distribution,
$q_{\tau,s,t}$ and the joint degree distribution, $p_{s,t}$ can be derived
from these information-rich function.

Figure 2: The joint excess joint degree distribution of topology $\tau$ is
created by selecting all $\tau$-edges and measuring the joint excess degrees
of the terminal vertices, excluding the motif the selected edge belongs to. In
the example, the selected edge (yellow dashed) has vertices of excess joint
degrees $(2,1)$ and $(1,1)$.

There is an $e_{\tau,s_{0}t_{0},s^{\prime}t^{\prime}}$ function for each
$\tau\in\vec{\tau}$; since, we can follow $\text{card}(\vec{\tau})$ distinct
edge-topologies. We assume that

$\sum_{s_{0}}\sum_{t_{0}}\sum_{s^{\prime}}\sum_{t^{\prime}}e_{\tau,s_{0}t_{0},s^{\prime}t^{\prime}}=1,$
(1)

with all indices running from $0,\dots,\infty$ unless otherwise stated. We can
recover the excess joint degree distribution, which is the probability that an
edge of topology $\tau$ is chosen at random and leads to a vertex of remaining
degree $(s_{0},t_{0})$, as

$q_{\tau,s_{0}t_{0}}=\sum_{s^{\prime}}\sum_{t^{\prime}}e_{\tau,s_{0}t_{0},s^{\prime}t^{\prime}}.$
(2)

We can invert $q_{\tau,s,t}$ to obtain $p_{st}$ for each
$\tau\in\\{\bot,\Delta\\}$. For the tree-triangle model this is given by

$\displaystyle p^{\bot}_{s,t}=$
$\displaystyle\frac{\bigg{[}q_{\bot,s-1,t}\big{/}s\bigg{]}}{\bigg{[}\sum\limits_{s^{\prime}t^{\prime}}q_{\bot,s^{\prime}-1,t}\big{/}s^{\prime}\bigg{]}},$
(3a) $\displaystyle p^{\Delta}_{s,t}=$
$\displaystyle\frac{\bigg{[}q_{\Delta,s,t-1}\big{/}t\bigg{]}}{\bigg{[}\sum\limits_{s^{\prime}t^{\prime}}q_{\Delta,s^{\prime},t-1}\big{/}t^{\prime}\bigg{]}},$
(3b)

which can be verified by inserting $q_{\bot,s,t}={(s+1)p_{s+1,t}}/{\langle
s\rangle}$ into Eq 3a or $q_{\Delta,s,t}={(t+1)p_{s,t+1}}/{\langle t\rangle}$
into Eq 3b where $\langle s\rangle$ is the average number of 2-cliques a
vertex belongs to

$\langle s\rangle=\sum_{s}\sum_{t}sp_{s,t},$ (4)

and similarly for $\langle t\rangle$. We remark that, in general, the inverted
distributions, $p_{s,t}^{\tau}$, are not numerically equivalent to one
another, nor do they necessarily contain all $(s,t)$ tuples. This is because,
the excess joint degree distributions are not equivalent to one another for
different topologies. To see this, $q_{\bot,s,t}$, (the joint excess
distribution obtained from randomly selected 2-clique edges) can yield no
information on the marginal distribution of how vertices with zero 2-cliques,
but non-zero 3-clique counts interact; and vice versa for $q_{\Delta,s,t}$;
since, we cannot follow a 2-clique edge to a vertex with no 2-cliques.
However, the actual joint degree distribution, $p_{s,t}$, can be recovered by
multiplying each partially observed joint degree distribution by the fraction
of the network (by vertex count) that was used to create the observation. For
instance, $p_{s,t}=n_{\tau}p^{\tau}_{s,t}$ where $n_{\tau}$ are the fraction
of vertices that belong to at least one $\tau$ motif. For a coherent model, we
have $n_{\tau}p^{\tau}_{s,t}=n_{\nu}p^{\nu}_{s,t}$ for a given $(s,t)$ joint
degree and for $\tau,\nu\in\vec{\tau}$. Each $n_{\tau}$ can readily be found
as one minus the fraction of vertices that do not belong to any $\tau$
cliques; for instance $n_{\bot}=1-\sum_{t}p_{0,t}$. When all vertices belong
to both 2- and 3-cliques, $n_{\tau}=n_{\nu}$, and the inverted joint
distributions coincide. Finally, it should be noted that we can obtain the
joint excess degree distribution between overall degrees $e_{j,k}$ from
$e_{\tau,s_{0}t_{0},s^{\prime}t^{\prime}}$ by application of a Kronecker delta

$\displaystyle e_{j,k}=$
$\displaystyle\sum_{s_{0},t_{0}}\sum_{s^{\prime}t^{\prime}}w_{\bot}e_{\bot,s_{0}t_{0},s^{\prime}t^{\prime}}\delta_{j,s_{0}+2t_{0}}\delta_{k,s^{\prime}+2t^{\prime}}$
$\displaystyle+\sum_{s_{0},t_{0}}\sum_{s^{\prime}t^{\prime}}w_{\Delta}e_{\Delta,s_{0}t_{0},s^{\prime}t^{\prime}}\delta_{j,s_{0}+2t_{0}+1}\delta_{k,s^{\prime}+2t^{\prime}+1}.$
(5)

The $+1$ terms in the final two delta functions account for the contribution
of the other edge within the 3-clique towards the overall degree. The
weighting factors $w_{\tau}$ account for the required renormalisation from the
independently normalised distributions to a single distribution. These are
given by

$w_{\bot}=\frac{\sum\limits_{st}sp_{st}}{\sum\limits_{st}sp_{st}+2\sum\limits_{st}tp_{st}},$
(6)

and

$w_{\Delta}=\frac{2\sum\limits_{st}tp_{st}}{\sum\limits_{st}sp_{st}+2\sum\limits_{st}tp_{st}}.$
(7)

This expression compresses the information available from the model and
collapses to the formulation derived in PhysRevLett.89.208701 . Additionally,
the overall degree distribution $p_{k}$ can readily be found from a joint
degree distribution PhysRevLett.103.058701

$p_{k}=\sum_{st}p_{st}\delta_{k,s+2t}.$ (8)

###  II.1 $G_{0}$ generating function

We can use $e_{\tau,s_{0}t_{0},s^{\prime}t^{\prime}}$ to find ensemble
properties by deriving a generating function formulation for the model.
Suppose that a vertex of joint degree $(s_{0},t_{0})$ is chosen from the
network. Now suppose that we traverse the $s_{0}$ 2-clique edges and record
the excess joint degrees $(s^{\prime},t^{\prime})$ of those neighbouring
vertices. Let
$\\{r_{\bot}\\}=\\{r_{\bot,s_{1},t_{1}},\dots,r_{\bot,s_{n},t_{n}}\\}$ be the
configuration of those neighbour joint degrees we might record; such that
there are $r_{\bot,s^{\prime},t^{\prime}}$ neighbours of excess joint degree
$(s^{\prime},t^{\prime})$. The probability of a given configuration of excess
joint degrees surrounding the focal vertex following $s_{0}$ 2-clique edges is

$P(\\{r_{\bot}\\}\mid
s_{0})=s_{0}!\prod_{s,t}\frac{1}{r_{\bot,st}!}\left(\frac{e_{\bot,s_{0}-1,t_{0},s,t}}{\sum\limits_{s,t}e_{\bot,s_{0}-1,t_{0},s,t}}\right)^{r_{\bot,st}}.$
(9)

Similarly, let us traverse each of the $2t_{0}$ 3-clique edges from the focal
vertex to the neighbours along the corners of each triangle and record the
excess joint degree of those vertices. We record the number of each unique
excess joint degree tuple as $r_{\Delta,s^{\prime}t^{\prime}}$ to find the
configuration $\\{r_{\Delta}\\}$ surrounding the focal vertex along all of its
3-clique edges. The probability of a given configuration is the multinomial

$P(\\{r_{\Delta}\\}\mid
t_{0})=2t_{0}!\prod_{s,t}\frac{1}{r_{\Delta,st}!}\left(\frac{e_{\Delta,s_{0},t_{0}-1,s,t}}{\sum\limits_{s,t}e_{\Delta,s_{0},t_{0}-1,s,t}}\right)^{r_{\Delta,st}}.$
(10)

The probability that a vertex with joint degree $(s_{0},t_{0})$ has a
neighbour configuration $\\{r_{\bot}\\}$ and $\\{r_{\Delta}\\}$ is simply the
product of these terms

$P(\\{r_{\bot}\\},\\{r_{\Delta}\\}\mid s_{0},t_{0})=P(\\{r_{\bot}\\}\mid
s_{0})P(\\{r_{\Delta}\\}\mid
t_{0})\delta\left(s_{0},\sum\limits_{s,t}r_{\bot,st}\right)\delta\left(2t_{0},\sum\limits_{s,t}r_{\Delta,st}\right),$
(11)

where $\delta(i,j)$ is the Kronecker delta that ensures only configurations
that sum to $s_{0}$ and $2t_{0}$ 2- and 3-clique edges, respectively are
retained. We can then generate this expression by summing over all of the
possible configurations along each edge topology

$\mathcal{P}_{s_{0},t_{0}}(z_{\bot},z_{\Delta})=\sum_{\\{r_{\bot}\\}}\sum_{\\{r_{\Delta}\\}}P\left(\\{r_{\bot}\\},\\{r_{\Delta}\\}\mid
s_{0},t_{0}\right)\left(\prod_{s,t}z_{\bot,st}^{r_{\bot,st}}\right)\left(\prod_{s,t}z_{\Delta,st}^{r_{\Delta,st}}\right),$
(12)

where

$\sum_{\\{r_{\tau}\\}}:=\sum_{r_{\tau,s_{1}t_{1}}}\sum_{r_{\tau,s_{1}t_{2}}}\cdots\sum_{r_{\tau,s_{n}t_{m}}},$
(13)

and where $s_{n}$ $(t_{m})$ is the largest 2-clique (3-clique) degree. The
arguments of $P_{s_{0},t_{0}}$ are vectors containing an element for each
$(s^{\prime},t^{\prime})$ pair in $e_{\tau,s_{0},t_{0},s^{\prime}t^{\prime}}$
for each clique topology $\\{\bot,\Delta\\}$ such that
$z_{\tau}=z_{\bot,s_{1}t_{1}},\dots,z_{\bot,s_{n}t_{m}}$. We then apply the
multinomial theorem

$(z_{1}+z_{2}+\dots+z_{m})^{n}=\sum_{k_{1}}\sum_{k_{2}}\cdots\sum_{k_{m}}\delta\left(\sum\limits_{t}^{m}k_{t},n\right)\frac{n!}{k_{1}!k_{2}!\cdots
k_{m}!}\prod_{t=1}^{m}z_{t}^{k_{t}},$ (14)

to Eq 12 and average over the probability of choosing a vertex of joint degree
$(s_{0},t_{0})$ to arrive at a generating function for a randomly selected
vertex

$\displaystyle\mathcal{G}_{0}(z_{\bot},z_{\Delta}^{2})=$
$\displaystyle\sum_{s_{0}}\sum_{t_{0}}p_{s_{0}t_{0}}\left(\frac{\sum\limits_{s,t}e_{\bot,s_{0}-1,t_{0},s,t}z_{\bot,st}}{\sum\limits_{s,t}e_{\bot,s_{0}-1,t_{0},s,t}}\right)^{s_{0}}\left(\frac{\sum\limits_{s,t}e_{\Delta,s_{0},t_{0}-1,s,t}z_{\Delta,st}}{\sum\limits_{s,t}e_{\Delta,s_{0},t_{0}-1,s,t}}\right)^{2t_{0}},$
(15a) $\displaystyle=$ $\displaystyle
G_{0}\left(\frac{\sum\limits_{s,t}e_{\bot,s_{0}-1,t_{0},s,t}z_{\bot,st}}{\sum\limits_{s,t}e_{\bot,s_{0}-1,t_{0},s,t}},\left[\frac{\sum\limits_{s,t}e_{\Delta,s_{0},t_{0}-1,s,t}z_{\Delta,st}}{\sum\limits_{s,t}e_{\Delta,s_{0},t_{0}-1,s,t}}\right]^{2}\right),$
(15b)

with $G_{0}(x,y)=\sum_{s}\sum_{t}p_{st}x^{s}y^{t}$. This is the fundamental
generating function of the model to obtain the properties of randomly selected
vertices. We can extract the properties of individual joint degrees from this
distribution as well as vertices with overall degree $k$, irrespective of
their joint degree by multiplying by a delta
$G_{0}(z)=G_{0}(z,z^{2})=\sum_{k}p_{k}z^{k}=\sum_{k}\sum_{st}p_{st}z^{s+2t}\delta_{k,s+2t}$.
In this case our expressions reduce to those found by Newman
PhysRevE.67.026126 .

We now calculate the size $S$ of the GCC by introducing a generating function
for each joint excess joint degree distribution. For this we need to generate
the probability of following an edge of a given topology to a vertex and that
vertex failing to belong to the GCC. For this purpose, let

$\displaystyle\mathcal{F}_{\tau,s_{0}t_{0}}(z_{\bot},z_{\Delta})$
$\displaystyle=\frac{\sum\limits_{s^{\prime}}\sum\limits_{t^{\prime}}e_{\tau,s_{0}t_{0},s^{\prime}t^{\prime}}z_{\bot,s^{\prime}t^{\prime}}^{s^{\prime}}z^{t^{\prime}}_{\Delta,s^{\prime}t^{\prime}}}{\sum\limits_{s^{\prime}}\sum\limits_{t^{\prime}}e_{\tau,s_{0}t_{0},s^{\prime}t^{\prime}}}.$
(16)

We define $u_{\bot,s_{0},t_{0}}$ to be the probability that a 2-clique
connected to a vertex of remaining degree $(s_{0},t_{0})$ at one end leads to
another vertex that does not belong to the giant component. Similarly, we
define $u_{\Delta,s_{0},t_{0}}$ to be the probability that a randomly selected
3-clique edge that has one vertex of joint excess degree $(s_{0},t_{0})$ fails
to connect to the GCC via the other end of the edge. We can write self-
consistency expressions for these probabilities by summing over all possible
neighbouring vertex excess joint degrees multiplied by the probability that
their edges fail to connect that vertex as

$\displaystyle u_{\bot,s_{0},t_{0}}=$
$\displaystyle\mathcal{F}_{\bot,s_{0}t_{0}}(g_{\bot},g^{2}_{\Delta}),$ (17a)
$\displaystyle u_{\Delta,s_{0},t_{0}}=$
$\displaystyle\mathcal{F}_{\Delta,s_{0}t_{0}}(g_{\bot},g^{2}_{\Delta}),$ (17b)

where $g_{\bot}$ and $g_{\Delta}$ are vectors with elements
$g_{\bot,s,t}=1-\phi+u_{\bot,s,t}\phi$ and
$g_{\Delta,s,t}^{2}=(1-\phi+u_{\Delta,s,t}\phi)^{2}-2\phi^{2}(1-\phi)u_{\Delta,s,t}(1-u_{\Delta,s,t})$,
which are the probabilities that a 2-clique edge and both edges in a 3-clique
fail to connect to the GCC PhysRevE.80.020901 , respectively for bond
occupancy probability $\phi$. The GCC is then found from

$S=1-G_{0}\left(\frac{\sum\limits_{s,t}e_{\bot,s_{0}-1,t_{0},s,t}g^{s}_{\bot,st}g_{\Delta,s+1,t-1}^{2(t-1)}}{\sum\limits_{s,t}e_{\bot,s_{0}-1,t_{0},s,t}},\left[\frac{\sum\limits_{s,t}e_{\Delta,s_{0},t_{0}-1,s,t}g_{\bot,s-1,t+1}^{s-1}g^{2t}_{\Delta,st}}{\sum\limits_{s,t}e_{\Delta,s_{0},t_{0}-1,s,t}}\right]^{2}\right).$
(18)

Note, some the indices contain $s\pm 1$ and $t\mp 1$; since, given an excess
2-clique joint degree, $(s,t)$, we obtain the excess 3-clique joint degree as
$(s+1,t-1)$; similarly, given an excess 3-clique joint degree of $(s,t)$, we
obtain the 2-clique excess joint degree as $(s-1,t+1)$. We extend this model
to GCM networks comprised of larger clique sizes in Appendix A.

###  II.2 Percolation threshold

Like other percolation models on random graphs, the model exhibits a 2nd order
critical point in $\phi$ with order parameter $S$. This indicates that the
order parameter is continuous at the transition, but its derivative is not.
The value of $\phi$ at the critical point, $\phi_{\text{crit}}$, can be found
by perturbing around the point $u_{\tau}=1$ $\forall\tau\in\vec{\tau}$, which
corresponds to the probability that all edges, regardless of topology, fail to
lead to the GCC. As $\phi\gtrsim\phi_{\text{crit}}$ each $u_{\tau}$ drops
below 1 and we can expand both expressions in Eqs 17 using a Taylor series
PhysRevLett.103.058701 ; 10.1371/journal.pone.0071321  about
$z_{\tau,st}=1-u_{\tau,st}$ to obtain

$\displaystyle z_{\tau,s_{0}t_{0}}\approx\ $ $\displaystyle
1-\mathcal{F}_{\tau,s_{0}t_{0}}(\bm{1},\bm{1})-\sum_{s}\sum_{t}\frac{\partial\mathcal{F}_{\tau,s_{0}t_{0}}}{\partial
g_{\bot,st}}\frac{\partial g_{\bot,st}}{\partial
z_{\bot,st}}\bigg{|}_{z_{\bot,st}=0}-2\sum_{s}\sum_{t}\frac{\partial\mathcal{F}_{\tau,s_{0}t_{0}}}{\partial
g_{\Delta,st}}\frac{\partial g_{\Delta,st}}{\partial
z_{\Delta,st}}\bigg{|}_{z_{\Delta,st}=0},$ $\displaystyle\approx\ $
$\displaystyle\phi\sum_{s}\sum_{t}\mathcal{F}^{(\bot,st)}_{\tau,s_{0}t_{0}}(\bm{1},\bm{1})z_{\bot,st}+2\phi(1+\phi-\phi^{2})\sum_{s}\sum_{t}\mathcal{F}_{\tau,s_{0}t_{0}}^{({{\Delta,st}})}(\bm{1},\bm{1})z_{\Delta,st}+(\mathcal{O}^{2}),$
(19)

where $\mathcal{F}_{\tau,st}^{{\nu,st}}$ is the partial derivative of
$\mathcal{F}_{\tau,st}$ with respect to $g_{\nu,st}$ and where we have treated
$g_{\tau,st}(z_{\tau,st})$ with $g_{\tau,st}(0)=1$. We have assumed that
$\mathcal{F}_{\tau,st}$ is well-behaved at $u=1$, with finite derivatives,
which does not hold for power law degree distributions. The system can be
expressed as a matrix equation $\bm{z}=\bm{\mathcal{F}}\bm{\phi}\cdot\bm{z}$
where $\bm{\mathcal{F}}$ is the Jacobian of $\mathcal{F}_{\tau,st}$
$\forall\tau\in\\{\bot,\Delta\\}$, $\bm{z}$ is the vector
$(z_{\bot,s_{1}t_{1}},\dots,z_{\Delta,s_{n}t_{m}})$ and

$\bm{\phi}=\sum_{i=1}\mathbf{c}_{i}^{T}\mathbf{g}\mathbf{c}_{i}\mathbf{c}_{i}^{T},$
(20)

where
$\mathbf{g}=(g^{\prime}_{\bot,s_{1}t_{1}},\dots,g^{\prime}_{\Delta,s_{n}t_{m}})$
is the vector of derivatives of $g_{\tau,st}(z_{\tau,st})$ at $z_{\tau,st}=0$,
$\mathbf{c}_{i}$ is the i-th basis vector of
$\mathbb{R}^{n_{\bot}m_{\bot}+n_{\Delta}m_{\Delta}}$ and $T$ denotes the
transpose. For a given $\phi$ the system exhibits non-trivial solution for
$\bm{z}$ when $\det(\bm{\mathcal{F}}\bm{\phi}-\bm{I})=0$ where $\bm{I}$ is the
identity matrix. The scalar value $\phi_{\text{crit}}$ (that uniquely
determines matrix $\bm{\phi}$) that satisfies this expression can then be
found from the characteristic equation.

###  II.3 Joint correlation functions

Traditionally, the information of the overall degree correlation properties of
a network is summarised by an assortativity coefficient,
$r^{\text{Pearson}}_{k}$. This can be obtained from the family of
$e_{\tau,st,s^{\prime}t^{\prime}}$ matrices by first finding $e_{j,k}$ via Eq
5 before evaluating

$r^{\text{Pearson}}_{k}=\frac{1}{\sigma^{2}_{q_{k}}}\sum_{jk}jk(e_{jk}-q_{j}q_{k}),$
(21)

where

$\sigma^{2}_{q_{k}}=\sum\limits_{k}k^{2}q_{k}-\left(\sum\limits_{k}kq_{k}\right)^{2},$
(22)

is the variance of the excess degree distribution of overall degrees $q_{k}$,
which in turn, is given by the summation of the rows $q_{k}=\sum_{j}e_{jk}$
PhysRevE.67.026126 . In general we expect a well-behaved assortativity
coefficient to take its values in $r\in[-1,1]$, with positive (negative)
coefficient indicating assortative (dissasortative) mixing and $r=0$
indicating random mixing. Finally we note that the expected statistical error
in an assortativity coefficient can be calculated by the jackknife method
efron_1979  or the bootstrap method 10.1214/ss/1177013815 .

In this section, we will examine possible extensions of the assortativity
coefficient to encapsulate the mixing properties of vertices based on joint
degree tuples. A natural definition is to simply extend the Pearson
correlation coefficient along each edge topology. For instance, if a network
comprises of 2- and 3-cliques, then we obtain two assortativtiy coefficients;
one for each $\tau\in\vec{\tau}=\\{\bot,\Delta\\}$ given by

$r^{\text{Pearson}}_{\tau}=\frac{1}{\sigma_{q_{\tau}}^{2}}\sum_{s,t,s^{\prime},t^{\prime}}sts^{\prime}t^{\prime}(e_{\tau,st,s^{\prime}t^{\prime}}-q_{\tau,st}q_{\tau,s^{\prime}t^{\prime}}),$
(23)

where $\sigma_{q_{\tau}}$ is the standard deviation of the excess joint degree
distribution $q_{\tau,st}$ in topology $\tau$, normalising the function by its
maximum value

$\sigma^{2}_{q_{\tau}}=\sum\limits_{s,t}({st})^{2}q_{\tau,st}-\left(\sum\limits_{st}stq_{\tau,st}\right)^{2}.$
(24)

The square of the standard deviation (the variance) is strictly positive
unless all vertices have the same joint degree, in which case the coefficient
is undefined.

Let us now introduce an additional assortativity coefficient that classifies
joint degree tuples independently and considers the inter-class mixing
properties. In this case the assortativity coefficient along each edge
topology is

$r_{\tau}^{\text{class}}=\frac{\sum\limits_{st}e_{\tau,st,st}-\sum\limits_{st}\left(\sum\limits_{s^{\prime}t^{\prime}}e_{\tau,st,s^{\prime}t^{\prime}}\right)\left(\sum\limits_{s^{\prime}t^{\prime}}e_{\tau,s^{\prime}t^{\prime},st}\right)}{1-\sum\limits_{st}\left(\sum\limits_{s^{\prime}t^{\prime}}e_{\tau,st,s^{\prime}t^{\prime}}\right)\left(\sum\limits_{s^{\prime}t^{\prime}}e_{\tau,s^{\prime}t^{\prime},st}\right)}.$
(25)

This assortativity measure examines the extent to which each class of joint
degree tuple mixes with itself in the network. More compactly, and with
generalisation to models with additional clique sizes, we write
PhysRevE.67.026126 ; PhysRevE.99.042306

$r_{\tau}^{\text{class}}=\frac{\text{Tr}[e_{\tau}]-||{e_{\tau}^{2}}||}{1-||{e_{\tau}^{2}}||},$
(26)

where $||x||$ is the sum of all elements of matrix $x$ and $\text{Tr}[x]$ is
the trace. This information is collected into a tuple of coefficients
$r^{\text{class}}=(r^{\text{class}}_{\bot},r^{\text{class}}_{\Delta},\dots,r^{\text{class}}_{\gamma})$
that describe the mixing patterns along each motif topology. Defining the
coefficient in this way is equivalent to mapping each distinct joint degree
tuple to a colour and then comparing the mixing patterns of vertex colours for
the different edge types balogh_palla_kryven_2020 . The drawback of this
coefficient is the loss of a metric between pairs of non-identical joint
degrees. For example, considering joint degree tuples that take the form
$(s,t)$, this definition would not find the joint degree $(1,2)$ any more
similar to $(1,3)$ than to $(10,20)$, despite the obvious difference in
connectivity. This classification-based coefficient is therefore stricter than
the degree-based one; since, it can exhibit no correlation, even when overall
degrees are correlated or when either of the marginal dimensions are
correlated.

There are other ways to define a assortativity coefficients PhysRevE.87.022801
; especially when there are vector attributes placed on vertices
pelechrinis_wei_2016 ; rabbany_eswaran_dubrawski_faloutsos_2017 . And, whilst
there is not a universal approach in the network science literature, most
methods embed the vector quantity in a metric space and compare distances
along each component. We leave further examination of assortativity
coefficients for future work.

###  II.4 Choosing $e_{\tau,st,s^{\prime}t^{\prime}}$ matrices

In order to use the formulation, we must find suitable
$e_{\tau,st,s^{\prime}t^{\prime}}$ matrices for each edge type
$\tau\in\\{\bot,\Delta\\}$ and the choice will determine the mixing properties
of the model. By far the simplest approach is to simply extract these
relationships from empirical networks, and we discuss this in section IV. In
this section, we overview manual processes instead.

A matrix $e_{\tau,st,s^{\prime}t^{\prime}}$ is a valid if it satisfies the
following conditions

$\displaystyle e_{\tau,st,s^{\prime}t^{\prime}}$ $\displaystyle\geq 0$ (27a)
$\displaystyle e_{\tau,st,s^{\prime}t^{\prime}}$
$\displaystyle=e_{\tau,s^{\prime}t^{\prime},st}$ (27b)
$\displaystyle\sum_{s^{\prime}t^{\prime}}e_{\tau,st,s^{\prime}t^{\prime}}$
$\displaystyle=q_{\tau,st},$ (27c) $\displaystyle\sum_{st}q_{\tau,st}$
$\displaystyle=1,$ (27d)

with $q_{\tau,st}\in[0,1]$. In addition, we require the further constraint
that the inversion of matrices $e_{\tau,st,s^{\prime}t^{\prime}}$ and
$e_{\nu,st,s^{\prime}t^{\prime}}$ for $p_{st}$ through Eqs 3a and 3b coincide,
thus constraining the excess joint degree distributions $q_{\tau}$ and
$q_{\nu}$.

For neutral mixing, in the absence of any inter-subgraph assortativity, the
joint excess joint degree distribution is given by the product of the excess
joint degrees

$e^{(\text{n})}_{\tau,st,s^{\prime}t^{\prime}}=q_{\tau,st}q_{\tau,s^{\prime}t^{\prime}}.$
(28)

For a perfectly assorted network, in which edges only connect vertices with
the same excess joint degree we have

$e^{(\text{a})}_{\tau,st,s^{\prime}t^{\prime}}=q_{\tau,st}\delta_{s,s^{\prime}}\delta_{t,t^{\prime}}.$
(29)

It is convenient to parameterise the extent of the assortativity with an
infinitesimal value $\epsilon$. We can readily construct these matrices from
the neutral case by re-distributing the row summations (excess joint degree
distribution) such that the diagonal is given by

$e^{(\text{a})}_{\tau,st,st}=\sum_{s^{\prime}t^{\prime}}e_{\tau,st,s^{\prime}t^{\prime}}-\epsilon\sum_{s^{\prime}t^{\prime}}q_{\tau,st}q_{\tau,s^{\prime}t^{\prime}}(1-\delta_{s,s^{\prime}}\delta_{t,t^{\prime}}),$
(30)

whilst off-diagonal entries, $s\neq s^{\prime}$ and $t\neq t^{\prime}$ are
given by

$e^{(\text{a})}_{\tau,st,s^{\prime}t^{\prime}}=\epsilon\sum_{s^{\prime}t^{\prime}}q_{\tau,st}q_{\tau,s^{\prime}t^{\prime}}(1-\delta_{s,s^{\prime}}\delta_{t,t^{\prime}}).$
(31)

For a perfectly dissasortative network, in which a vertex with a given joint
degree will never connect to a neighbour with the same joint degree we have a
zero diagonal element

$e^{(\text{d})}_{\tau,st,st}=0,\qquad\forall(s,t)$ (32)

The off-diagonal elements are given by a linear combination of the row sums
whose coefficients are $\pm 1$. The row sum $q_{i}$ and column sum $q_{j}$ for
element $e_{ij}$ always appears with a positive sign; whilst, the coefficients
of the other row and column sums are such that they sum to zero for any given
row or column. For instance, consider the following 3x3 example

$e^{(\text{d})}_{\tau}=\begin{pmatrix}0&\frac{q_{1}+q_{2}-q_{3}}{2}&\frac{q_{1}-q_{2}+q_{3}}{2}\\\
\frac{q_{1}+q_{2}-q_{3}}{2}&0&\frac{-q_{1}+q_{2}+q_{3}}{2}\\\
\frac{q_{1}-q_{2}+q_{3}}{2}&\frac{-q_{1}+q_{2}+q_{3}}{2}&0\\\ \end{pmatrix},$
(33)

where $q_{i}$ corresponds to the sum of elements on row $i$. The elements are
divided by 2 since there are 2 non-zero values to distribute the row sum over
per row. We present more details in Appendix B.

We can relax the strict zero diagonal property of these matrices by
introducing infinitesimal values $\epsilon_{i}$ along the diagonal such as

$e^{({\text{r}})}_{\tau}=\epsilon+e^{(\text{d})}_{\tau}$ (34)

where $\epsilon$ is the matrix

$\epsilon=\begin{pmatrix}\epsilon_{1}&\frac{-\epsilon_{1}-\epsilon_{2}+\epsilon_{3}}{2}&\frac{-\epsilon_{1}+\epsilon_{2}-\epsilon_{3}}{2}\\\
\frac{-\epsilon_{1}-\epsilon_{2}+\epsilon_{3}}{2}&\epsilon_{2}&\frac{\epsilon_{1}-\epsilon_{2}-\epsilon_{3}}{2}\\\
\frac{-\epsilon_{1}+\epsilon_{2}-\epsilon_{3}}{2}&\frac{\epsilon_{1}-\epsilon_{2}-\epsilon_{3}}{2}&\epsilon_{3}\\\
\end{pmatrix}.$ (35)

##  III Simulation

With a theoretical framework in place to study the percolation properties of
clustered random graphs that exhibit a given mixing pattern, we would now like
to discuss how to simulate such networks. The data input to the simulation
framework is a coherent set of $e_{\tau}$ matrices for a set of permissible
subgraph topologies $\vec{\tau}$. We have two options: i) extract $e_{\tau}$
from an empirical network; or ii) write $e_{\tau}$ analytically. We have
discussed how to generate matrices analytically in section II.4 and so, we now
overview a workflow to examine empirical networks. The workflow is summarised
in Fig 3 and defined as follows. For a given empirical graph, an edge-disjoint
(non-overlapping) clique cover can be used to cover the network by a
predetermined basis set of clique subgraphs. Once the cover has been applied
to the network, the joint excess joint degree distribution,
$e_{\tau,{s,t,\dots},{s^{\prime},t^{\prime},\dots}}$ can be gathered; from
which, the excess joint degree distributions, $q_{\tau,{s,t,\dots}}$ and the
joint degree distribution $p_{s,t,\dots}$ can be found. Random graphs can be
created according to the GCM prescription and these graphs can be post-
processed using a stochastic MCMC rewiring algorithm, which we introduce in
Section III, to yield the correct mixing patterns.

Figure 3: An overview, from left to right, of the workflow to create a random
graph of arbitrary size with statistically equivalent clustering and inter-
clique mixing properties from a substrate network. A substrate network is
covered by an edge-disjoint clique cover before the mixing matrices are
extracted along each edge type. The excess joint degree distributions can then
be calculated and inverted to find the joint degree distribution $p_{st}$.
Random networks can then be created using the configuration model and the
resulting synthetic network can be rewired by our MCMC stochastic algorithm.

Given the centrality of covering empirical networks with an edge-disjoint set
of cliques to the application of this work, it is appropriate to comment on
existing methods in the literature. To the author’s knowledge, there are two
distinct edge-disjoint clique covers: i) the EECC
burgio_arenas_gomez_matamalas_2021  developed by Burgio et al and ii) the MPCC
PhysRevE.105.044314  developed by Mann et al. Both of these covers are
stochastic heuristics that cover a network in cliques up to some tunable
maximum size such that a given edge in the substrate network belongs to a
single clique. The EECC optimises on the number of cliques that are used to
cover the graph, aiming for as few as possible, such a cover is a maximal
cover. Conversely, the MPCC gives precedence to larger cliques, aiming to
retain as much local structure as possible for high-degree vertices and
neighbourhoods. Details on these algorithms can be found in
burgio_arenas_gomez_matamalas_2021 ; PhysRevE.105.044314 .

Once we have covered the empirical network in motifs, we can then extract the
joint excess joint degree distributions,
$e_{\tau,{s,t,\dots},{s^{\prime},t^{\prime},\dots}}$ from it; subsequently,
the excess joint degree distributions, $q_{\tau,{s,t,\dots}}$ and the joint
degree distribution $p_{s,t,\dots}$ can be obtained. Random graphs $G(V,E)$
can now be created according to the GCM algorithm by drawing vertex joint
degrees from $p_{s,t,\dots}$ to create a sequence of joint degrees of length
$|V|$ PhysRevE.82.066118 ; PhysRevE.103.012313 ; PhysRevE.103.012309 . We note
that the number of vertices $|V|$ is independent of the size of the empirical
network, and so, this method can be used to enlarge network datasets for
subsequent analysis whilst holding the clustering coefficient and the motif
distribution fixed. At this point, the inter-motif correlations are not fixed.
We then create a list $l_{\tau}$ for each distinct clique size
$\tau\in\vec{\tau}$ and enter the identity of each vertex once into those
lists for each distinct motif that it belongs to of size $\tau\in\vec{\tau}$.
The lists of vertex identities are then randomised. We now construct the
network. For each clique we draw without replacement the appropriate number of
vertices from each list; equivalently, we partition the list into segments
with length equal to the number of vertices for a given size clique. We repeat
this for all motifs until the lists are exhausted; all edges have been placed
in the network and the algorithm is terminated.

In general, it is highly unlikely that the same joint degree sequence will be
drawn from $p_{s,t,\dots}$ or that the identity of the vertices that belong to
a motif will be identical upon repeat of this algorithm. This means that we
will obtain a different graph $G^{\prime}(V,E)$ that belongs to the same
ensemble of random graphs $\mathcal{G}$ that share equivalent statistical
properties. We remark that the accidental selection of the same vertex from a
clique list, or the accidental combining of one or more edges between motifs
(such that they are not edge-disjoint) is possible; however, both events are
increasingly uncommon as $|V|$ increases for a fixed joint degree distribution
and so will not effect the statistics of the model.

The graphs created via the GCM are absent of correlations between motifs;
although they may display overall degree correlations due to the sizes of each
clique. To inject a correlation structure between the motifs, we introduce an
MCMC rewiring algorithm whose target distribution yields the correct
$e_{\tau}$ joint excess joint degree distribution for all $\tau\in\vec{\tau}$
that appear in the model. At each stage of the algorithm a clique size is
selected from $\vec{\tau}$. Next, two vertices that belong to cliques of that
size are selected at random from different motifs in the network. The edges
that connect the vertices to the clique are broken. The two vertices are
exchanged and new edges are created to rewire the incident vertices to their
new locations. For example, in Fig 4, a network comprising 2- and 3-cliques is
rewired according to our MCMC algorithm. In this case, vertices $u_{3}$ and
$v_{3}$ are swapped to yield an alternative configuration that preserves the
joint degree sequence. To achieve this, edges
$(u_{1},u_{3}),(u_{2},u_{3}),(v_{1},v_{3})$, and $(v_{2},v_{3})$ are broken
and edges $(u_{1},v_{3}),(u_{2},v_{3}),(v_{1},u_{3})$ and $(v_{2},u_{3})$ are
created. The proposed swap is then accepted with probability
$P=\text{min}(1,\pi)$ where

$\pi=\frac{e_{\Delta,k_{u_{1}},k_{v_{3}}}e_{\Delta,k_{u_{2}},k_{v_{3}}}e_{\Delta,k_{v_{1}},k_{u_{3}}}e_{\Delta,k_{v_{2}},k_{u_{3}}}}{e_{\Delta,k_{u_{1}},k_{u_{3}}}e_{\Delta,k_{u_{2}},k_{u_{3}}}e_{\Delta,k_{v_{1}},k_{v_{3}}}e_{\Delta,k_{v_{2}},k_{v_{3}}}}$
(36)

where $k_{i}=(s_{i},t_{i})$ in $e_{\tau,k_{i},k_{j}}$ is the excess joint
degree of vertex $i$ following topology $\tau$. This Metropolis-Hastings
condition is the quotient of the product of the proposed joint excess joint
degree distributions following the swap and their non-permuted counterparts.

Figure 4: The rewiring process for two triangle corners in the network. Edges
$(u_{1},u_{3}),(u_{2},u_{3}),(v_{1},v_{3})$, and $(v_{2},v_{3})$ are broken
and new edges $(u_{1},v_{3}),(u_{2},v_{3}),(v_{1},u_{3})$ and $(v_{2},u_{3})$
are created.

For larger clique sizes, when vertex $u_{1}$ is swapped with $v_{1}$, the
Metropolis-Hastings acceptance condition is written

$\pi=\prod_{i}\frac{e_{\tau,k_{u_{1}},k_{v_{i}}}e_{\tau,k_{v_{1}},k_{v_{i}}}}{e_{\tau,k_{u_{1}},k_{u_{i}}}e_{\tau,k_{v_{1}},k_{v_{i}}}}$
(37)

where index $i$ runs over the number of edges being broken in each clique
(i.e. its size minus 1).

It is clear to see that this dynamics conserves the joint degree sequence and
is ergodic on the set of graphs having that joint degree sequence. By this, we
mean that any graph with the joint degree sequence can be obtained in a finite
number of Monte Carlo updates. The algorithm also satisfies detailed balance
for any graph configuration $\lambda$ that occurs with probability

$p_{\lambda}=\prod_{\tau}\prod_{i=1}^{M_{\tau}}e_{\tau_{i}}$ (38)

where the product is over all distinct clique sizes $\tau$ in the network and
all edges that belong to cliques of that size $i=1,\dots,M_{\tau}$ in the
graph. Detailed balance can be shown by considering the transition probability
between any two states $\lambda\rightarrow\lambda^{\prime}$ by evaluating
$p_{\lambda}P(\lambda\rightarrow\lambda^{\prime})=p_{\lambda^{\prime}}P(\lambda^{\prime}\rightarrow\lambda)$.
With ergodicity and detailed balance satisfied, the MCMC algorithm yields a
graph that has the required correlation structure as its fixed point.

##  IV Discussion

###  IV.1 $k$-regular graphs

Consider a random graph ensemble that consists of networks that have both 2-
and 3-clique subgraphs. Let the joint degree distribution $p_{st}$ be given by
$(5,1)=1/3,\ (3,2)=1/3$ and $(1,3)=1/3$ such that the global clustering
coefficient and overall vertex degree are fixed for all networks
($k$-regular). In this example, by following 2-clique edges, we can reach
vertices with joint excess degrees $(0,3),\ (4,1)$ and $(2,2)$; i.e. these are
the excess joint degrees in $q_{\bot,st}$. Similarly, by following triangle
edges to a corner, we can reach vertices with joint excess degrees $(1,2),\
(0,5)$ and $(3,1)$; which constitute the excess joint degrees of
$q_{\Delta,st}$.

The mixing matrices $e_{\tau}^{\text{neutral}}$ for the uncorrelated model are
given by

$\NiceMatrixOptions{code-for-first-row=\scriptstyle,code-for-first-
col=\scriptstyle,}e_{\bot}^{(\text{n})}=\pNiceArray[first-row,first-
col]{ccc}&\rotatebox[origin={c}]{60.0}{$\scriptstyle(0,3)$}\rotatebox[origin={c}]{60.0}{$\scriptstyle(4,1)$}\rotatebox[origin={c}]{60.0}{$\scriptstyle(2,2)$}\\\
(0,3)\frac{1}{81}\frac{5}{81}\frac{3}{81}\\\
(4,1)\frac{5}{81}\frac{25}{81}\frac{15}{81}\\\
(2,2)\frac{3}{81}\frac{15}{81}\frac{9}{81}\\\ $ (39)

and

$\NiceMatrixOptions{code-for-first-row=\scriptstyle,code-for-first-
col=\scriptstyle,}e_{\Delta}^{(\text{n})}=\pNiceArray[first-row,first-
col]{ccc}&\rotatebox[origin={c}]{60.0}{$\scriptstyle(3,1)$}\rotatebox[origin={c}]{60.0}{$\scriptstyle(1,2)$}\rotatebox[origin={c}]{60.0}{$\scriptstyle(5,0)$}\\\
(3,1)\frac{16}{144}\frac{24}{144}\frac{8}{144}\\\
(1,2)\frac{24}{144}\frac{36}{144}\frac{12}{144}\\\
(5,0)\frac{8}{144}\frac{12}{144}\frac{4}{144}\\\ $ (40)

For example, element $e^{(\text{n})}_{\bot,(0,3),(0,3)}=1/9\times 1/9$ since
$1/9$ of all tree edges lead to a vertex with joint degree $(1,3)$. Similarly,
element $e^{(\text{n})}_{\Delta,(3,1),(3,1)}=4/12\times 4/12$ since four of
the 12 3-clique edges lead to vertices with joint excess degree $(3,1)$.

We can make assorted matrices $e_{\tau}^{(\text{a})}$ by pushing the
probability mass into the diagonal entries and setting the off-diagonal terms
to a infinitesimal value, according to Eqs 30 and 31.

$\NiceMatrixOptions{code-for-first-row=\scriptstyle,code-for-first-
col=\scriptstyle,}e_{\bot}^{(\text{a})}=\pNiceArray[first-row,first-
col]{ccc}&\rotatebox[origin={c}]{60.0}{$\scriptstyle(0,3)$}\rotatebox[origin={c}]{60.0}{$\scriptstyle(4,1)$}\rotatebox[origin={c}]{60.0}{$\scriptstyle(2,2)$}\\\
(0,3)\frac{9}{81}-\epsilon\frac{5}{81}-\epsilon\frac{3}{81}\epsilon\frac{5}{81}\epsilon\frac{3}{81}\\\
(4,1)\epsilon\frac{5}{81}\frac{45}{81}-\epsilon\frac{5}{81}-\epsilon\frac{15}{81}\epsilon\frac{15}{81}\\\
(2,2)\epsilon\frac{3}{81}\epsilon\frac{15}{81}\frac{27}{81}-\epsilon\frac{3}{81}-\epsilon\frac{15}{81}\\\
$ (41)

$\NiceMatrixOptions{code-for-first-row=\scriptstyle,code-for-first-
col=\scriptstyle,}e_{\Delta}^{(\text{a})}=\pNiceArray[first-row,first-
col]{ccc}&\rotatebox[origin={c}]{60.0}{$\scriptstyle(3,1)$}\rotatebox[origin={c}]{60.0}{$\scriptstyle(1,2)$}\rotatebox[origin={c}]{60.0}{$\scriptstyle(5,0)$}\\\
(3,1)\frac{48}{144}-\epsilon\frac{24}{144}-\epsilon\frac{8}{144}\epsilon\frac{24}{144}\epsilon\frac{8}{144}\\\
(1,2)\epsilon\frac{24}{144}\frac{72}{144}-\epsilon\frac{24}{144}-\epsilon\frac{12}{144}\epsilon\frac{12}{144}\\\
(5,0)\epsilon\frac{8}{144}\epsilon\frac{12}{144}\frac{24}{144}-\epsilon\frac{8}{144}-\epsilon\frac{12}{144}\\\
$ (42)

where $0\lesssim\epsilon$. This procedure establishes an equivalence class
between the excess joint degree distributions of the assorted and neutral
correlation matrices; and consequently, between the inverted joint degree
distribution. In other words, the joint degree distribution and the excess
joint degree distributions, the clustering coefficients and the overall degree
of the vertices of the networks that can be created will all be equivalent;
however, their inter-subgraph mixing patterns exhibit different structure. The
expected size of the GCC for the neutral and assorted mixing matrices is shown
in the top plot of Fig 5. Despite being $k$-regular, assortativity between the
joint degrees changes the percolation properties of the network, reducing the
critical point.

Next, we examine the effect of disassortative mixing in a $k$-regular network
with joint degree distribution $p_{st}$ given by $(5,1)\approx 0.282$,
$(3,2)\approx 0.374$ and $(1,3)\approx 0.344$. The neutral mixing matrices are
given by Eq 28 whilst the dissasorted matrices are given by Eq 33. In Fig 5,
we plot the results of Eq 18 for the two mixing patterns, observing that
inter-subgraph mixing leads to an (albeit small) increase in the critical
point of the model, in line with traditional theory of dissasortative mixing
between overall degrees PhysRevLett.89.208701 . However, due to the vertices
having the same overall degrees, other theoretical models could not exhibit
this behaviour. Importantly, our model provides evidence that assortativity
plays an important role in the observed percolation dynamics, even when the
clustering coefficient is fixed (in addition to the overall degree); settling
a longstanding debate on the subject PhysRevE.80.020901 ; PhysRevE.98.062314 ;
PhysRevE.81.066114 .

Figure 5: The expected size of the GCC for neutral and correlated random
networks in $k$-regular graphs with 2- and 3-clique clustering among vertices
with overall degree $s+2t=7$ and fixed clustering coefficients. Curves are the
theoretical results of Eq 18.

###  IV.2 Empirical graphs

For our next analysis, we examine the correlation properties of a series of
empirical networks. The datasets are the network of hyperlinks between weblogs
on US politics 10.1145/1134271.1134277 ; a network of coauthorships between
scientists posting preprints on the Condensed Matter E-Print Archive between
1995 and 1999 doi:10.1073/pnas.98.2.404 ; the Astrophysics collaboration
network konect:2017:ca-AstroPh ; konect:leskovec107 ; Twitter
NIPS2012_7a614fd0 ; friendship data of Facebook users konect:2017:facebook-
wosn-links ; viswanath09 ; the social network of Douban konect:2017:douban ;
konect:socialcomputing ; user–user friendship relations from Brightkite
konect:2017:loc-brightkite_edges ; konect:cho2011 ; the collaboration graph of
authors of scientific papers from the arXiv’s Astrophysics (astro-ph) section
konect:2017:ca-AstroPh ; konect:leskovec107 ; the network of connections
between autonomous systems of the Internet konect:2017:topology ;
konect:zhang05 ; the network of autonomous systems of the Internet connected
with each other from the CAIDA project konect:2017:as-caida20071105 ;
konect:2017:as-caida20071105 ; coappearance network of characters in the novel
Les Miserables knuth_knuth_knuth_2009 ; the C. Elegans network
watts_strogatz_1998 .

In Fig 6, we plot the size of the GCC for a series of empirical networks under
different clique covers. Scatter points are the result of Monte Carlo
simulations of bond percolation on the original graph; whilst the (red) dashed
lines are the theoretical results of traditional edge-based theory (Eq 9 of
PhysRevLett.89.208701 ) for 2-clique covered networks; finally, the (green)
dashed and dotted lines are the results of 2- and 3-clique covered networks in
combination with Eq 18. Across the series, it is immediately evident that the
higher-order mixing properties are crucial in capturing the percolation
response of the empirical network. In each case, our 2-3 clique model is in
strong agreement with the dynamics of the empirical graph. It is interesting
that in some case, such as C. Elegans and the LesMis network, that the
2-clique theory is in close agreement with the simulation; whilst, for social
networks, such as the Twitter and Douban networks, the edge-based theory
offers poor agreement with simulation.

Figure 6: The size of the GCC for a series of empirical network datasets. The
scatter points are the average of Monte Carlo bond percolation on the original
network. The (red) dashed lines are the expected theoretical curves when using
a 2-clique approximation PhysRevLett.89.208701  ; whilst the (green) dashed
and dotted curves are the results for the 2- and 3-clique theory of Eq 18. In
each case, including higher-order interactions from the empirical network
yields a better theoretical representation.  Figure 7: Scatter plots of the
assortativity coefficients $r_{\tau}^{\text{Pearson}}$ and
$r_{\tau}^{\text{Class}}$ given by Eqs 23 and 25, respectively of 6 empirical
networks with a linear regression model fit. The shaded regions are the
bootstrap confidence interval for the regression estimate whilst the markers
are the data points. To collect the data, the networks were covered using an
MPCC before the $e_{\tau}$ mixing matrices were extracted and used to
calculate the assortativity coefficient along each edge topology.

In Fig 7 we analyse the inter-subgraph mixing properties for a range of
empirical networks. For each dataset we first find a full MPCC cover of the
network before extracting the mixing matrices $e_{\tau}$. These matrices are
then used to calculate each $r_{\tau}^{\text{Pearson}}$ and
$r_{\tau}^{\text{Class}}$ through Eqs 23 and 25, respectively. Scatter points
are the values of the empirically obtained correlation coefficients for each
clique size, whilst lines are fitted regression plots with shaded confidence
intervals. We observe a range of mixing patterns in the networks, in
particular we notice that vertices that belong to cliques of size $2<c<7$ tend
to exhibit assortative mixing. For example, the increase in the Pearson
coefficient for vertices in the Twitter network. However, larger clique sizes
tend to be dissasortatively mixed, for instance, cliques of size larger than
$10$ in the condensed matter coauthorship network. This is almost certainly
due to the smaller number of larger cliques in these finite sized networks.

##  V Conclusion

In this paper we have introduced a theoretical framework, including a
generating function formulation and rewiring algorithm, that accounts for
higher-order mixing patterns in networks with clustering. We have applied our
framework to random graph ensembles, conceptual models and to empirical graphs
by use of edge-disjoint clique covers from the literature. Members (graphs)
within the ensembles have the following fixed statistics: joint degree
distribution, clustering coefficient and inter-motif correlation structure. We
found that our model offers excellent agreement with Monte Carlo simulation
and yields a significant improvement beyond the results of ordinary edge-based
theory. We found that non-trivial correlation structure of empirical networks
as a function of cliques size, indicating the presence of structure at
different scales.

We have provided an analytical description of how to write the input data to
this model, $e_{\tau}$ matrices as well as how to mine them from empirical
networks. To simulate such networks we have introduced an MCMC rewiring
algorithm that, in conjunction with the well-known GCM algorithm, yields
synthetic networks with the statistically equivalent mixing properties as the
empirical models.

This method can be used to draw samples of synthetic network datasets with
fixed statistical properties that are representative of the original data in
terms of clustering and mixing patterns. Our formulation is focused around
clique subgraphs; however, we anticipate the generalisation of this work to
hypergraphs in future work, and we believe this would offer an increasingly
realistic model of empirical networks and greater understanding of their
properties. We anticipate the application of our model to a variety of
dynamical processes and networked datasets, which we leave for future work.

For subgraphs larger than 2-cliques, our method inherently captures beyond
nearest-neighbour connections within the inter-subgraph correlation structure.
However, it falls short of being a sufficient model of long range inter-
subgraph correlations; a feature that would allow higher-order equivalence
classes among tensor quantities that capture the $n$-hop correlation
structures to be written. By increasing the rigidity of the local environment
around a vertex, the degrees of freedom of the GCM model and consequently the
number of possible graphs in the random graph ensemble is reduced. Properties
such as connectedness, core-periphery sub-structure, beyond nearest-neighbour
inter-subgraph correlations and modularity are crucial mesoscale aspects of
network topology whose control would reduce the number of configurations
within an ensemble of random graphs, and therefore yield models that are
increasingly representative of empirical data.

The role of correlations of this nature has implications for the spread of
epidemic diseases, and their non-trivial thresholds, and should be studied in
this context PhysRevE.103.062308 ; PhysRevE.103.042307 . We leave this for
future work.

##  VI ACKNOWLEDGMENTS

This work was partially supported by the UK Engineering and Physical Sciences
Research Council under grant number EP/N007565/1 (Science of Sensor Systems
Software).

## References

  * (1) A. P. Millán, J. J. Torres, and G. Bianconi, “Explosive higher-order kuramoto dynamics on simplicial complexes,” Phys. Rev. Lett., vol. 124, p. 218301, May 2020.
  * (2) I. Iacopini, G. Petri, A. Barrat, and V. Latora, “Simplicial models of social contagion,” Nature Communications, vol. 10, no. 1, 2019.
  * (3) C. Giusti, R. Ghrist, and D. S. Bassett, “Two’s company, three (or more) is a simplex,” Journal of Computational Neuroscience, vol. 41, no. 1, p. 1–14, 2016.
  * (4) U. Alvarez-Rodriguez, F. Battiston, G. F. de Arruda, Y. Moreno, M. Perc, and V. Latora, “Evolutionary dynamics of higher-order interactions in social networks,” Nature Human Behaviour, vol. 5, no. 5, p. 586–595, 2021.
  * (5) M. E. J. Newman, “Properties of highly clustered networks,” Phys. Rev. E, vol. 68, p. 026121, Aug 2003.
  * (6) M. E. J. Newman, “Random graphs with clustering,” Phys. Rev. Lett., vol. 103, p. 058701, Jul 2009.
  * (7) M. E. J. Newman, “Mixing patterns in networks,” Phys. Rev. E, vol. 67, p. 026126, Feb 2003.
  * (8) M. E. J. Newman, “Assortative mixing in networks,” Phys. Rev. Lett., vol. 89, p. 208701, Oct 2002.
  * (9) S. G. Balogh, G. Palla, and I. Kryven, “Networks with degree–degree correlations are special cases of the edge-coloured random graph,” Journal of Complex Networks, vol. 8, no. 4, 2020.
  * (10) J. C. Miller, “Percolation and epidemics in random clustered networks,” Phys. Rev. E, vol. 80, p. 020901, Aug 2009.
  * (11) J. P. Gleeson, “Bond percolation on a class of clustered random networks,” Physical Review E, vol. 80, Oct 2009.
  * (12) A. Hackett, S. Melnik, and J. P. Gleeson, “Cascades on a class of clustered random networks,” Phys. Rev. E, vol. 83, p. 056107, May 2011.
  * (13) J. P. Gleeson, S. Melnik, and A. Hackett, “How clustering affects the bond percolation threshold in complex networks,” Phys. Rev. E, vol. 81, p. 066114, Jun 2010.
  * (14) P. Mann, V. A. Smith, J. B. O. Mitchell, and S. Dobson, “Symbiotic and antagonistic disease dynamics on networks using bond percolation,” Phys. Rev. E, vol. 104, p. 024303, Aug 2021.
  * (15) P. Mann, V. A. Smith, J. B. O. Mitchell, and S. Dobson, “Percolation in random graphs with higher-order clustering,” Phys. Rev. E, vol. 103, p. 012313, Jan 2021.
  * (16) T. Hasegawa and S. Mizutaka, “Structure of percolating clusters in random clustered networks,” Phys. Rev. E, vol. 101, p. 062310, Jun 2020.
  * (17) P. Mann, V. A. Smith, J. B. O. Mitchell, and S. Dobson, “Random graphs with arbitrary clustering and their applications,” Phys. Rev. E, vol. 103, p. 012309, Jan 2021.
  * (18) P. Mann, V. A. Smith, J. B. O. Mitchell, C. Jefferson, and S. Dobson, “Exact formula for bond percolation on cliques,” Phys. Rev. E, vol. 104, p. 024304, Aug 2021.
  * (19) A. Vázquez and M. Weigt, “Computational complexity arising from degree correlations in networks,” Phys. Rev. E, vol. 67, p. 027101, Feb 2003.
  * (20) M. E. J. Newman, S. H. Strogatz, and D. J. Watts, “Random graphs with arbitrary degree distributions and their applications,” Physical Review E, vol. 64, no. 2, 2001.
  * (21) B. Karrer and M. E. J. Newman, “Random graphs containing arbitrary distributions of subgraphs,” Phys. Rev. E, vol. 82, p. 066118, Dec 2010\.
  * (22) P. Mann, V. A. Smith, J. B. O. Mitchell, and S. Dobson, “Degree correlations in graphs with clique clustering,” Phys. Rev. E, vol. 105, p. 044314, Apr 2022.
  * (23) G. Burgio, A. Arenas, S. Gómez, and J. T. Matamalas, “Network clique cover approximation to analyze complex contagions through group interactions,” Communications Physics, vol. 4, no. 1, 2021.
  * (24) M. E. J. Newman and C. R. Ferrario, “Interacting epidemics and coinfection on contact networks,” PLOS ONE, vol. 8, pp. 1–8, 08 2013.
  * (25) B. Efron, “Computers and the theory of statistics: Thinking the unthinkable,” SIAM Review, vol. 21, no. 4, p. 460–480, 1979.
  * (26) B. Efron and R. Tibshirani, “Bootstrap Methods for Standard Errors, Confidence Intervals, and Other Measures of Statistical Accuracy,” Statistical Science, vol. 1, no. 1, pp. 54 – 75, 1986.
  * (27) G. T. Cantwell and M. E. J. Newman, “Mixing patterns and individual differences in networks,” Phys. Rev. E, vol. 99, p. 042306, Apr 2019.
  * (28) N. Litvak and R. van der Hofstad, “Uncovering disassortativity in large scale-free networks,” Phys. Rev. E, vol. 87, p. 022801, Feb 2013.
  * (29) K. Pelechrinis and D. Wei, “Va-index: Quantifying assortativity patterns in networks with multidimensional nodal attributes,” PLOS ONE, vol. 11, no. 1, 2016.
  * (30) R. Rabbany, D. Eswaran, A. W. Dubrawski, and C. Faloutsos, “Beyond assortativity: Proclivity index for attributed networks (prone),” Advances in Knowledge Discovery and Data Mining, p. 225–237, 2017.
  * (31) S. Mizutaka and T. Hasegawa, “Disassortativity of percolating clusters in random networks,” Phys. Rev. E, vol. 98, p. 062314, Dec 2018.
  * (32) L. A. Adamic and N. Glance, “The political blogosphere and the 2004 u.s. election: Divided they blog,” in Proceedings of the 3rd International Workshop on Link Discovery, LinkKDD ’05, (New York, NY, USA), p. 36–43, Association for Computing Machinery, 2005.
  * (33) M. E. J. Newman, “The structure of scientific collaboration networks,” Proceedings of the National Academy of Sciences, vol. 98, no. 2, pp. 404–409, 2001.
  * (34) “arxiv astro-ph network dataset – KONECT,” Oct. 2017.
  * (35) J. Leskovec, J. Kleinberg, and C. Faloutsos, “Graph evolution: Densification and shrinking diameters,” ACM Trans. Knowl. Discov. from Data, vol. 1, no. 1, pp. 1–40, 2007.
  * (36) J. Leskovec and J. Mcauley, “Learning to discover social circles in ego networks,” in Advances in Neural Information Processing Systems (F. Pereira, C. Burges, L. Bottou, and K. Weinberger, eds.), vol. 25, Curran Associates, Inc., 2012.
  * (37) “Facebook friendships network dataset – KONECT,” Oct. 2017.
  * (38) B. Viswanath, A. Mislove, M. Cha, and K. P. Gummadi, “On the evolution of user interaction in Facebook,” in Proc. Workshop on Online Soc. Netw., pp. 37–42, 2009.
  * (39) “Douban network dataset – KONECT,” Oct. 2017.
  * (40) R. Zafarani and H. Liu, “Social computing data repository at ASU,” 2009.
  * (41) “Brightkite network dataset – KONECT,” Oct. 2017.
  * (42) E. Cho, S. A. Myers, and J. Leskovec, “Friendship and mobility: User movement in location-based social networks,” in Proc. Int. Conf. on Knowl. Discov. and Data Min., pp. 1082–1090, 2011.
  * (43) “Internet topology network dataset – KONECT,” Oct. 2017.
  * (44) B. Zhang, R. Liu, D. Massey, and L. Zhang, “Collecting the Internet AS-level topology,” SIGCOMM Comput. Communication Review, vol. 35, no. 1, pp. 53–61, 2005.
  * (45) “Caida network dataset – KONECT,” Oct. 2017.
  * (46) D. E. Knuth, D. E. Knuth, and D. E. Knuth, The Stanford graphbase: A platform for combinatorial computing.  Addison-Wesley, 2009.
  * (47) D. J. Watts and S. H. Strogatz, “Collective dynamics of ‘small-world’ networks,” Nature, vol. 393, no. 6684, p. 440–442, 1998.
  * (48) P. Mann, V. A. Smith, J. B. O. Mitchell, and S. Dobson, “Two-pathogen model with competition on clustered networks,” Phys. Rev. E, vol. 103, p. 062308, Jun 2021.
  * (49) P. Mann, V. A. Smith, J. B. O. Mitchell, and S. Dobson, “Cooperative coinfection dynamics on clustered networks,” Phys. Rev. E, vol. 103, p. 042307, Apr 2021.

##  Appendix A larger cliques

We now generalise the model developed in section II.1 to account for GCM
networks that contain larger clique subgraphs beyond 2- and 3-cliques. Suppose
that a vertex of joint degree $(s_{0},t_{0},\dots,c_{0})$ is chosen from a
random graph that is comprised of edge-disjoint cliques up to some maximum
size $C$. Now suppose that we traverse each of the edges of a given clique
size $\tau$ and record the excess joint degrees
$(s^{\prime},t^{\prime},\dots,c^{\prime})$ of the neighbouring vertices and
then repeat this for all $\tau\in\vec{\tau}$. Let
$\\{r_{\tau}\\}=\\{r_{\tau,s_{1},t_{1},\dots,c_{1}},\dots,r_{\tau,s_{n},t_{n},\dots,c_{n}}\\}$
be the configuration of those neighbour joint degrees we might record; such
that there are $r_{\tau,s^{\prime},t^{\prime},\dots,c^{\prime}}$ neighbours of
excess joint degree $(s^{\prime},t^{\prime},\dots,c^{\prime})$ and so on. The
probability of a given configuration of excess joint degrees surrounding the
focal vertex is

$\displaystyle\mathcal{P}_{s_{0},\dots,c_{0}}=$ $\displaystyle\ p_{s_{0}\dots
c_{0}}\sum_{\\{r_{\bot}\\}}\cdots\sum_{\\{r_{C}\\}}\prod_{\tau}(\tau-1)k_{0}!\prod_{s,\dots,c}\frac{1}{r_{\tau,s\dots
c}!}$
$\displaystyle\times\left(\frac{e_{\tau,s_{0},\dots,k_{0}-1,\dots,c_{0},s,\dots,k,\dots,c}}{\sum\limits_{s,\dots,c}e_{\tau,s_{0},\dots,k_{0}-1,\dots,c_{0},s,\dots,k,\dots,c}}z_{\tau,s\dots
c}\right)^{r_{\tau,s\dots
c}}\delta\left((\tau-1)k_{0},\sum\limits_{s,\dots,c}r_{\tau,s\dots c}\right)$
(43)

and after application of Eq 14 we have

$\displaystyle\mathcal{P}_{s_{0},\dots,c_{0}}=p_{s_{0}\dots
c_{0}}\prod_{\tau}\left(\frac{\sum\limits_{s,\dots,c}e_{\tau,s_{0},\dots,k_{0}-1,\dots,c_{0},s,\dots,k,\dots,c}z_{\tau,s,\dots,c}}{\sum\limits_{s,\dots,c}e_{\tau,s_{0},\dots,k_{0}-1,\dots,c_{0},s,\dots,c}}\right)^{(\tau-1)k_{0}}$
(44)

which yields the average properties of the first order neighbours of a vertex
with joint degree $(s_{0},\dots,c_{0})$ that has been randomly selected from
the network. This is an information-rich distribution function from which we
can extract the average properties $\mathcal{P}_{k}$ of vertices with overall
degree $k$ by filtering out those joint degrees that sum to $k$ as

$\mathcal{P}_{k}=\sum_{s}\dots\sum_{c}\mathcal{P}_{s\dots
c}\delta\left(k,\sum_{\tau}(\tau-1)\cdot k_{\tau}\right)$ (45)

where $k_{\tau}$ is the number of $\tau$ cliques the focal vertex belongs to.

##  Appendix B Perfectly dissasortative matrices

Our definition of perfectly dissasortative symmetric matrices includes a zero
diagonal $e_{ii}=0$ whose off diagonal elements $e_{ij}$ sum to a given value
$\sum_{ij}e_{ij}=q_{i}$ such that all entries are positive. Once the structure
of a given matrix has been found, the excess joint degree labels can be
arranged to fit the required mixing patterns. Extending the example from Eq
33, consider a 5x5 matrix $M_{\tau}$ with zero diagonal entries and with off-
diagonal elements $e_{i,j}$ given by

$\displaystyle e_{1,2}=$ $\displaystyle\ {q_{1}+q_{2}+q_{3}+q_{4}-q_{5}}$
(46a) $\displaystyle e_{1,3}=$ $\displaystyle\
{q_{1}+q_{2}+q_{3}-q_{4}-q_{5}}$ (46b) $\displaystyle e_{1,4}=$
$\displaystyle\ {q_{1}-q_{2}-q_{3}+q_{4}+q_{5}}$ (46c) $\displaystyle
e_{1,5}=$ $\displaystyle\ {q_{1}-q_{2}-q_{3}-q_{4}+q_{5}}$ (46d)
$\displaystyle e_{2,1}=$ $\displaystyle\ {q_{1}+q_{2}+q_{3}+q_{4}-q_{5}}$
(47a) $\displaystyle e_{2,3}=$ $\displaystyle\
{q_{1}+q_{2}+q_{3}-q_{4}+q_{5}}$ (47b) $\displaystyle e_{2,4}=$
$\displaystyle\ {-q_{1}+q_{2}-q_{3}+q_{4}-q_{5}}$ (47c) $\displaystyle
e_{2,5}=$ $\displaystyle\ {-q_{1}+q_{2}-q_{3}-q_{4}+q_{5}}$ (47d)
$\displaystyle e_{3,1}=$ $\displaystyle\ {q_{1}+q_{2}+q_{3}-q_{4}-q_{5}}$
(48a) $\displaystyle e_{3,2}=$ $\displaystyle\
{q_{1}+q_{2}+q_{3}-q_{4}+q_{5}}$ (48b) $\displaystyle e_{3,4}=$
$\displaystyle\ {-q_{1}-q_{2}+q_{3}+q_{4}-q_{5}}$ (48c) $\displaystyle
e_{3,5}=$ $\displaystyle\ {-q_{1}-q_{2}+q_{3}+q_{4}+q_{5}}$ (48d)
$\displaystyle e_{4,1}=$ $\displaystyle\ {q_{1}-q_{2}-q_{3}+q_{4}+q_{5}}$
(49a) $\displaystyle e_{4,2}=$ $\displaystyle\
{-q_{1}+q_{2}-q_{3}+q_{4}-q_{5}}$ (49b) $\displaystyle e_{4,3}=$
$\displaystyle\ {-q_{1}-q_{2}+q_{3}+q_{4}-q_{5}}$ (49c) $\displaystyle
e_{4,5}=$ $\displaystyle\ {q_{1}+q_{2}+q_{3}+q_{4}+q_{5}}$ (49d)
$\displaystyle e_{5,1}=$ $\displaystyle\ {q_{1}-q_{2}-q_{3}-q_{4}+q_{5}}$
(50a) $\displaystyle e_{5,2}=$ $\displaystyle\
{-q_{1}+q_{2}-q_{3}-q_{4}+q_{5}}$ (50b) $\displaystyle e_{5,3}=$
$\displaystyle\ {-q_{1}-q_{2}+q_{3}+q_{4}+q_{5}}$ (50c) $\displaystyle
e_{5,4}=$ $\displaystyle\ {q_{1}+q_{2}+q_{3}+q_{4}+q_{5}}$ (50d)

The matrix $e_{\tau}=0.25M$ is a perfectly dissasortative mixing matrix that
satisfies the constraints on the system given by Eqs 27. Following this
prescription, only matrices with odd order, $n$, can be found. This is
because, for odd $n$ the number of positive elements per row is $n-1$, which
must be even for pairs of elements whose net sum is zero to be found. It is
unclear if other methods of generating perfectly dissasortative matrices that
satisfy Eqs 27 exist.