MCMC Sampling of Directed Flag Complexes with Fixed Undirected Graphs

Constructing null models to test the significance of extracted information is a crucial step in data analysis. In this work, we provide a uniformly sampleable null model of directed graphs with the same (or similar) number of simplices in the flag complex, with the restriction of retaining the underlying undirected graph. We describe an MCMC-based algorithm to sample from this null model and statistically investigate the mixing behaviour. This is paired with a high-performance, Rust-based, publicly available implementation. The motivation comes from topological data analysis of connectomes in neuroscience. In particular, we answer the fundamental question: are the high Betti numbers observed in the investigated graphs evidence of an interesting topology, or are they merely a byproduct of the high numbers of simplices? Indeed, by applying our new tool on the connectome of C. Elegans and parts of the statistical reconstructions of the Blue Brain Project, we find that the Betti numbers observed are considerable statistical outliers with respect to this new null model. We thus, for the first time, statistically confirm that topological data analysis in microscale connectome research is extracting statistically meaningful information.


Motivation
Inspecting spike correlation data, i.e. groups of neurons preferably spiking together, is a natural approach to investigate information processing in the brain.As the data comes in the form of a simplicial complex, the choice of using topological data analysis (TDA) for in-depth inspection is natural.This approach has seen recent success in [7,18], but is severely limited by the data available: acquiring large-scale, yet detailed, in-vivo neuronal activity records remains infeasible, which restricts analysis to experiments of small scale or to in-silico simulations.
Why investigate connectomes?Analyzing connectomes, i.e. the directed adjacency graph on the level of single neurons and synaptic connections between them instead, might mitigate these shortcomings: not only do they offer additional information due to the directionality added, it also makes massive amounts of data readily available.One kind would be recently developed largescale sparse connectome reconstructions [2,12].These reconstructions stochastically generate connectomes from a large -but in comparison to biology still limited -set of biological facts (e.g. by using the position and orientation of different neuron types and their respective probabilities to form synaptic links), resulting in networks which match these biological statistics.Even more excitingly, there has been rapid progress in the field of dense reconstructions [9,13,17,22].In this approach the connectome of small volumes of actual brain tissue is determined.Electron microscopes together with advanced data processing methods achieve a resolution which allows one not only to map neurons (as with traditional light-based microscopy, which is unable to pick up synapses) but even allows to locate and count synaptic vesicles which are several magnitudes smaller.This scientific and technical progress enables an accurate and detailed reconstruction of networks of neurons and synapses far beyond researchers wildest dreams just a mere decade ago.With this kind of indisputable biological ground truth data becoming increasingly available, the accompanying development of analytical methods is most desirable.As it has seen success in the related analysis of spike correlation data, one of the potential methods trying to establish itself is TDA.
Connectomic analysis with TDA.When thinking about graphs through the lens of TDA, it is not immediately clear how to go about extracting complex information: in contrast to simplicial complexes gained by activity correlation, connectomic graphs themselves are just a simplicial complex of maximum dimension one (as vertices are of dimension 0 and edges are of dimension 1) and are therefore topologically not very interesting.
Thus one approach to enable extraction of information is to first construct the so-called flag complex.This represents a very specific case of a simplicial complex which contains all simplices compatible with the given graph.The space generated by this approach has a geometrical realisation and could be interpreted as the space describing potential signal information flow in the brain, where higher-dimensional simplices represent robust, dense, yet directed paths.
Once the flag complex has been computed by greedily searching for all simplices present, further topological analysis is possible, e.g. the calculation of Betti numbers in different dimensions in order to get an idea of the topology of the entire structure [15]; -analysis, i.e. looking for connected simplicial pathways [16]; or the simplex completion probability [19], which endeavours to shed light on simplex generation mechanisms.
The above approach was pioneered in [15], where, amongst others, the dense reconstruction of the neural network of the roundworm Caenorhabditis Elegans (C.Elegans) [4,20] was investigated (cited in Figure 1, left).Here, two findings stand out: compared to an Erdős-Rényi graph (ER-graph) of equal size and edge density, the connectome of C. Elegans seems to have, in all higher dimensions, • more simplices and • higher Betti numbers.But are these really two independent findings?Nontrivial homology expressed by nonzero or high Betti numbers in a given dimension requires, by definition, simplices in this dimension, which the ER-graph used for comparison does not provide.Undoubtedly aware of the problem, the authors of [15] never claimed to have found statistically interesting topology.
In spite of the traction gained by aboves TDA-based approach of connectomic analysis, the same fundamental question remains: What kind of Betti numbers are to be expected in such simplex-rich graphs?

The Need for Sampleable Null Models
As a purely theoretical analysis currently seems to be intractable, we seek at least experimental, statistical insight.For that approach, one needs a suitable null model -in this case the set of all graphs with the identical simplex count seems appropriate.Whilst a parameterized null-model might be envisioned as well, we kept to a nonparametric version as in line with other research in network science.In particular, the condition of keeping the identical simplex counts enforces the same number of vertices and edges already, making this null model a more restricted version of the equally nonparametric ER-graphs so popular in network science.
Our desired insight into the expected topology may then be acquired straightforwardly by uniformly sampling a large number of graphs from that null model and analyzing their topology with the established toolchain.
Sampling from simplex-rich graphs is difficult.The problem lies in quick, efficient and uniform sampling of graphs from this set.Whilst constructing graphs with a specified number of simplices is straightforward, the resulting graphs are not uniformly distributed (or it would be very difficult to construct and prove).
Rejection sampling, i.e. the approach to first sample many ER-graphs and reject all those which don't have the right amount of simplices would work in principle, but is computationally infeasible.Connectomes are so simplex rich that they lie hundreds of standard deviations away from the expected mean and thus, in practice, every single graph would be rejected.This can be easily observed in Figure 1, right, where the distribution in a very similar scenario is depicted.
Another approach takes inspiration from stub matching algorithms, which maintain the degree sequence as required for the configuration model, yet another popular approach in network science.In this approach all edges are cut, leading to a stub per edge at each participating vertex.These stubs are then randomly paired, which results in retaining the number of edges per vertex, the so called degree sequence, of the original graph.Adapted for simplicial complexes (which are temporarily interpreted as hypergraphs), this approach has the potential to rapidly sample directed simplicial complexes with a comparable number of maximal simplices.Yet, for us it is unsuitable (see next paragraph), as the resulting complexes are very, very far from being flag1 .Thus rejection sampling the results becomes, once again, computationally infeasible.
With no easy solution in sight, it is prudent to first investigate preceding research.

Related Work.
Due to the novelty of the construction of directed flag complexes, there are no competing null models for this specific scenario.The closest established null model to the one we envision here is the so called Simplicial Configuration Model [23].It is an extension of the configuration model, which maintains the degree sequence of an undirected graph (i.e. each vertex maintains the same number of edges).As such, it maintains the number of maximal simplices and additionally for each node the number of maximal simplices attached to it.While useful in many other scenarios (e.g. the initially mentioned [18]), this model does not satisfy our needs: besides being undirected and focusing on maximal simplices, it would not respect the property that our original simplicial complex is flag.The last difference especially is crucial, as it is already known that flag complexes are really specific edge cases of simplicial complexes.Thus any statistical difference between the connectomes and this null model might be attributed due the samples not being flag.
Nonetheless, their work demonstrated how a restrained random walk with small and potentially erroneous steps could be interpreted as a Markov Chain Monte Carlo approach which guarantees, under some reasonable assumptions, uniform sampling.

The Proposed Approach
Restricting the null model further.Straightforwardly adapting [23] turned out to be unviable: the difficulty lies not only with the additional directionality, but primarily with the additional complexity introduced with the upwardclosure of flag-complexes (see Remark 2.7).
In order to start somewhere, we decided to conceptually factorize modifications of directed graphs into two distinct ways: connectivity and directionality.This allows us to treat the two resulting subproblems separately, resulting in the original problem becoming more approachable.It further allows for more detailed insight of the origin of the topological anomaly, should one be present.On the other hand, should one isolated aspect exhibit significant nonrandom behaviour already, one would expect that to still be visible in the general model.
The connectivity is an undirected property.The notion of simplices reduces to the notion of cliques and the flag complex becomes the clique complex.Our previously envisioned null model would then require graphs with the same number of cliques.We believe that the approach used in [23] could be adapted to this model, although that has yet to be investigated in detail.
The directionality without the connectivity then respects a given undirected graph (i.e. the information whether  and  are connected at all), and is flexible only with respect to the direction of edges and the position of reciprocated edges (later referred to as double edges).In this work, the focus lies exclusively on directionality.

Proposed Null Model
We thus propose a null model which enables uniform sampling of graphs that, with respect to a given graph, have: • the same number of vertices and directed edges, • a similar amount of simplices in the flag complex, • the same underlying undirected graph.
Here 'similar' means an allowed, user-adjustable deviation from the simplex count, e.g."within 1% of the original graph".The deviation is not meant as a parameter, more like a small error deemed acceptable, as maintaining the precise simplex count is extremely difficult.
Overview of the utilized approach.Even with the additional restriction of retaining the undirected graph, rejection sampling as explained above remains infeasible (see Figure 1, right).So far our research indicates that approaches on the level of simplicial complexes do not maintain the flag property.Thus, as the flag complex is completely determined by the graph, we have to modify the graph itself, then compute the flag complex and observe the resulting changes.For this approach, it is crucial that local changes in the graph result only in local changes in the flag complex.Exploiting this locality property is what keeps our approach computationally feasible.
Inspired by [23], we utilize, superficially, the same technique: a restricted random walk with resampling on the set of desired graphs.After a sufficient number of random steps, the point we end up at is uniformly distributed on the set of desired graphs.This is under the assumption that each point is as likely to enter as to exit (double stochastic) and all the desired points are connected (irreducible).
Our small steps conceptually consist of the so-called simple moves: flipping edges and moving double edges (depicted in Figure 3).The edges to be modified are chosen uniformly at random, and independent of previous steps, making the process a Markov Chain.As every change is reversible with the same probability, we automatically get double stochasticity.These moves do have the potential to change the simplex count though, which over time might accumulate, and subsequently lead to a loss of the desired property of a similar simplex count.
We thus reject changes which would lead the random walk too far astray and simply resample the current point instead (see Figure 2, red and blue arrows).
A common problem with restricted random walks is restricting too much, and thus splitting up the desired set into connected components, which are not interconnected by small steps anymore.In terms of Markov processes, this is expressed as the Markov Chain no longer being irreducible.To tackle this problem we utilize two tricks which extend the approach used in [23].The first is rather general and could be used in similar settings: we allow the random walk to stray as far away from our originally desired set as necessary for it to become connected again (the yellow zone in Figure 2).This implies an additional, subsequent rejection sampling step, which is necessary to filter out undesired graphs, but incurs significant additional computational cost.Second, we allow for larger modifications, the so called complex moves (depicted in Figure 6 and  7).Though independently discovered, they bear quite some resemblance to the stub-shuffling approach mentioned above (although restricted to small parts of the graph).Besides helping with irreducibility, they have the additional benefit of speeding up the random walk by combining dozens of small, local steps into one big move.This not only speeds up computation, but also guarantees some stability of the simplex count as well.
The choice of a rather basic MCMC approach over more refined versions is based on the following observations: Hamilton-MCMC or the Wang-Landau algorithm and all derivations thereof are not applicable, as our setting is very discrete in nature.Metropolis-Hastings and Reversible Jump MCMC both coincide with our approach in this particular setting: the acceptance ratio of Metropolis-Hastings binarizes and thus boils down to our chosen strategy of resampling.Should the dimension of the space searched change at any point then reversible jumps are the method of choice (see e.g.[14]).But in the case investigated here the number of vertices and edges and thus the dimension is fixed, which negates any benefit of reversible jumps.

Structure of the Paper
In Section 2, we briefly recall basic definitions and results on graphs, homology, and MCMC-sampling.Furthermore, we describe general methods to adapt MCMC-sampling to sets only very indirectly describable, as is the case here.
In Section 3 we apply this general framework from the viewpoint of Markov Chains to our problem, first describing and analyzing very simple transition steps useful for a theoretical approach to the problem and then more complex ones, which have, amongst others, computational advantages.We further discuss boundary relaxations necessary to achieve irreducibility of the transition matrix.
In Section 4 we adopt the algorithmic standpoint and briefly discuss our high-performance, publicly available, Rust-based implementation.We then empirically investigate the mixing speed of this process by analyzing the Hamming distance over time and conducting a  2 test for sample independence.
In Section 5 we apply this null-model to the connectome of C. Elegans, showing for the first time that this connectome seems to have an interesting, statistically abnormal topology compared to the null model designed and acquired in the previous sections.We further apply it to parts of the connectome of the Blue Brain Project, demonstrating that even large graphs with millions of edges are feasible with our implementation.
In Section 6 we discuss potential spin-offs, side problems, generalisations and the path to a null-model which is no longer connectivity-restrained.

Graphs
The following definitions of graphs and subgraphs are primarily stated for selfcontainedness and to establish notation.Most readers may skip them if aware of the subject.Details of the definition of cliques in directed graphs (as induced complete subgraphs) and simplices (as not-necessarily induced complete, acyclic, oriented subgraphs) are necessary to follow Section 3 and Section 4 in all technical detail, though.The same holds for lemmata on the locality property of flag-complexes.
• A simple undirected graph is a pair  = ( , ) of a finite set of vertices  and a symmetric relation  ⊆ ( ×  )\Δ  , where Δ  = {,  |  ∈  }.We refer to the set of simple undirected graphs by Graph.By excluding Δ  from  we exclude self-loops from  to itself.
• A simple directed graph is a pair  = ( , ) of a finite set of vertices  and a relation  ⊆ ( ×  )\Δ  , where Δ  = {(, ) |  ∈  }.We refer to the set of simple directed graphs by DiGraph.Here an edge (, ) ∈  is interpreted as an edge from  to .We call an edge a single edge, when (, ) ∈ , but (, ) ∉ .It is not forbidden that both (, ) and (, ) lie in , in that case we call it a double edge.
• We call pr() of a directed graph  the underlying undirected graph.Cliques are usually defined over undirected graphs, where their definition is analogous: pr() is a clique in an undirected graph.Simplices could also be defined as complete, oriented, acyclic subgraphs of , respectively oriented, acyclic subgraphs of a clique in .Remark 2.4.While Example 2.3 might, on first glance, imply that -simplices in (+1)-cliques are quite common, they are actually very rare in higher dimensions.Indeed, let  =  2 + 2 be the number of edges in a -simplex, leading us to 2  potential (single-edge only)-configurations of the edges in a (+1)-clique.At the same time there are "only" ( + 1)! potential total orders on the set of vertices, leading to a ratio of of -simplices per (+1)-clique.That function shrinks quite quickly (values rounded): This effect is somewhat countered by occurrences of double edges and the rich number of lower-dimensional cliques in a high-dimensional clique.Still, there is a good chance that flipping a single edge destroys a high-dimensional simplex.
Remark 2.5.As Example 2.3 demonstrates, one clique may exhibit multiple simplices if the graph contains double edges.Indeed, an -clique with  = (  2 ) double edges (i.e.every edge is a double edge) contains exactly ! simplices of dimension -1, as every order on the vertices could be expressed with the available edges.With -1 double edges, the number of (-1)-simplices is exactly ! 2 .With two missing double edges the minimum amount of simplices drops to ! 4 , whereas the maximum amount of simplices seems to follow the integer sequence A058298 2 .Once there are 3 double edges missing, one can enforce a 3-cycle, disabling all (-1)-simplices within this -clique.This leads to the following table:

Directed Flag Complexes and their Locality Property
Definition 2.6 (Directed Flag Complexes and the Number of Simplices).Let  = ( , ) be a simple directed graph and  the dimension of the highestdimensional simplex in .The directed flag complex ℱ() is a tuple of length +1 containing in the respective position the set of all -dimensional simplices of .
We then define where   yields the number of simplices of dimension  in ℱ().
Remark 2.7.Flag complexes are simplicial complexes, but the reverse is not necessarily true.A flag complex could be understood (or even defined) as the maximal simplicial complex over a given graph.This has a significant, and difficult to predict, impact on small modifications done to the complex: adding a high-dimensional simplex to a standard simplicial complex merely implies adding lower dimensional simplices as well, as every subset of a simplex must be contained (this property is referred to as downward closure).This is in contrast to flag complexes, where adding a (high-dimensional) simplex might not only require to add new lower-dimensional simplices, but then, as additional consequence of edges added due to the downward closure, new higher-dimensional simplices might appear in the graph and are added to the flag complex as well.One could refer to this property as upward-closure.This additional upward closure makes it difficult to adapt samplers designed for simplicial complexes to flag complexes.
The construction of flag-complexes is quite local: if there is a small local change in the directed graph, then the flag complex of that graph will remain the same for anything not in the direct vicinity of the changes.This definition straightforwardly extends to multiple edges, with  ⊆ : 2 https://oeis.org/history?seq=A058298 Lemma 2.9.Let  = ( , ) and  ′ = ( ,  ′ ) be two directed graphs with the same underlying undirected graph.Let  ∶= ( ′ ⧵ ) ∪ ( ⧵  ′ ) be the (symmetric) difference between these two graphs.Then Proof.To keep the equations more readable, we shortcut  ∶= [nbhd()] and  ∶=  ′ [nbhd()].We start with the following complement-defining equations: Our desired statement holds as soon as ℱ()  = ℱ()  as we can then substitute ℱ()  with ℱ() ⧵ ℱ() in the second equation.Let w.l.o.g. ∈ ℱ()  .Then  is not a simplex in ℱ() and thus pr() not a (undirected) clique in pr().As the underlying undirected graphs of  and  are the same, pr() is also not a (undirected) clique in pr().As a simplex is by definition also a clique in the underlying undirected graph,  cannot be a simplex in ℱ().
Still,  ∈ ℱ( ′ ) as all the edges of  are disjunct to  and are thus both in  and  ′ .Thus  ∈ ℱ()  .The same holds true for the other direction.This result might look technical, but is crucial for the implementation, as it allows for highly efficient calculations of updated simplex counts after local modifications of the graph.See Section 4.1.3for the application.

Unrestricted MCMC-Sampling
We recap the basics of Markov Chain Monte Carlo (MCMC) sampling, as it is the foundation of everything onward.For a concise introduction already in the context of graph sampling we refer to [5].For a more exhaustive, general work with focus on mixing behaviour we refer to [8].Definition 2.10 (Time-Independent Markov Chains).Let  be a  × stochastic matrix, i.e. it is non-negative and the sum of every row is 1.Let  be the starting distribution over  ∶= {1, … ,  }.This forms a stochastic process, the time-independent Markov Chain ℳ(,  ,  ).It yields for time  ∈ ℕ a random variables   which is distributed according to    .
One could describe (,  ) as a weighted directed graph by interpreting  as an adjacency matrix, but to avoid confusion with the graphs of our null model (which are points in ) we stick to the notion and notation of stochastic matrices.

Definition 2.11 (Aperiodic, Doubly Stochastic, Irreducible). Transition matrices may have the following properties:
aperiodic If there is a point  ∈  from which we start and return to at times  1 ,  2 , … then the greatest common divisor of  1 ,  2 , … shall be 1 for aperiodicity of .On the level of graphs: the greatest common divisor of the length of all cycles shall be 1.

doubly stochastic
The transpose of  is stochastic as well, i.e. the sum over each row and each column is one.The corresponding graph property would be regularity.
irreducible For each ,  ′ ∈  there exists a  ∈ ℕ such that  ′ is reachable from  in  steps.In graph terms it would be strongly connected.
A matrix which is aperiodic and irreducible is called ergodic (as its Markov Chain is ergodic).
Ergodic Markov Chains have an equilibrium distribution, and if the transition matrix is doubly stochastic, this equilibrium distribution is uniform: Theorem 2.12 (Uniform Equilibrium Distribution).Let  be a finite set of states and  a transition matrix which is ergodic and doubly stochastic.Then the random variable   defined by the Markov Chain ℳ(,  ,  ) converges, no matter the starting distribution , to the uniform distribution, i.e. with  = 1 || and   being the distribution of   : This theorem has an important practical application: let us assume a random walk obeys the rules of an ergodic, doubly stochastic Markov Chain and starts somewhere in .Then the theorem guarantees that the distribution it would end up with is approximately uniform.This approximation becomes more accurate the longer the random walk.
This process can be used to uniformly sample from  (assuming an initial starting point is supplied).The approach is called Markov Chain Monte Carlo Sampling.

Mixing Time and Sampling Distance
Waiting till infinity is terribly boring.For practical applications concerned with uniformly sampling from  the Markov Chain is usually run and subsampled every  steps, where  is called the sampling distance.Describing the (asymptotic) relationship between time and distance to the uniform distribution is referred to as studying the mixing behaviour, where rigorous theoretical results are often non-trivial.A more modest and applied approach -and thus the one pursued in this work -tries to estimate autocorrelation empirically.

Restricted MCMC-Sampling
Sometimes one strives to sample from a desired subset  ⊆ .If  can only be described in a very implicit way, the straightforward adaption of MCMCsampling is to first try a step, with the same probabilities as before, and to check if one is still in the desired subset .One has to be careful how to treat a misstep into   : simply rejecting a forbidden move would result in boundary effects on  -boundary regions would be underrepresented as they are harder to reach (only from one side).The resulting distribution of the MCMC-sampling process would no longer be uniform, as the transition matrix would no longer be doubly stochastic.
A simple solution is the resampling-strategy: instead of merely rejecting a forbidden move and treating it as it has never happened, one resamples the current state instead.This results in border regions which area hard to enter also being harder to leave, ensuring -at least up to the connected component -the desired uniform sampling quality.Definition 2.13 (Restricted Transition Matrix).Let  be a transition matrix over a finite set of states  and  the subset of allowed states.Then with ,  ′ ∈  the restricted transition matrix   is defined by This trick is known in the context of degree-maintaining Markov Chains as "swap-and-hold" or "trial-swap" (see [1,5]) and ensures that the resulting restricted transition matrix retains most of its MCMC-properties: Lemma 2.14.Let  be a transition matrix which is aperiodic and doubly stochastic.Then the restricted transition matrix   is aperiodic and doubly stochastic as well.
Note that irreducibility is not so easily inherited -connected components of , even if they were connected by , are not necessarily connected by   anymore (See Figure 2).Indeed, this is the main concern throughout this paper.We developed the following general strategy to deal with it: Definition 2.15 (-connecting subset).We call  ⊆  -connecting iff  is a subset of one irreducible component of (,   ).
Instead of running the Markov Chain ℳ(,   , ) and potentially being stuck in one irreducible component of , we run the Markov Chain on an connecting subset , i.e.ℳ(,   , ) and subsequently rejection-sample the resulting samples for membership in  (i.e. a sample  is accepted if  ∈  and otherwise rejected).
As  is a subset of the connected component of , this results in uniformly sampled samples of .
To optimize this approach one should strive for smaller, ideally even the minimum -connecting subset.This is, in general, a hard problem and no universal solution is in sight.

Constructing the Transition Matrix
As motivated in the introduction, we pursue to sample uniformly from the following null model: let  be a simple, directed graph with simplex count  = ().We desire to sample from all graphs with the same underlying undirected graph, simplex counts close to , i.e. between  − ≤  ≤  + and exactly the same amount of vertices ( 0 ) and edges ( 1 ).With the formal definition introduced below, we refer to this set as   +  − ().But how do we sample uniformly?In Section 2.2 we have introduced the theoretical background, restricted MCMC-processes with resampling and extended boundaries.
In this section we will discuss its application to our problem.This is primarily done by designing the required transition matrix.Here we show that it sticks to the desired subset   +  − () and exhibits the desired properties of aperiodicity, double stochasticity and irreducibility necessary for an MCMC-process.Admittedly the last property is the most tricky one, instead of rigorous proofs we present a detailed discussion.This is achieved by first investigating the so-called simple moves, a simpleto-understand and analyzable set of atomic moves which are enough to modify the directionality of a graph.They are complemented by the so-called complex moves, which have both theoretical and computational benefits by making use of the combinatorial nature of simplices.Both do not respect the simplex count by themselves though, so the random walk must be restricted as explained in the previous section.As these restriction might split   +  − () into multiple disconnected subsets, the last part of this chapter discusses potential choices of connecting boundaries.
This section is complemented by Section 4, where everything is viewed from a more practical, algorithmic perspective and the mixing time is determined experimentally.With  ∞ 0 () we refer to the (not-really) boundaries that are 0 respectively ∞ for all dimensions higher than 1.This describes the set of all directed graphs with the same underlying graph as  and the same number of (directed) edges, but no further restrictions on the simplex count.The following properties are straightforward to check and ensure that the image and coimage of SEF and DEM coincide, namely  ∞ 0 (): Lemma 3.3.Let  be transformed to  ′ by either SEF or DEM:

The move retains the number of directed edges:
Definition 3.4 (Unrestricted Simple Move Transition Matrix).Let  = ( , ) be a directed graph.We define the unrestricted transition matrices  u SEF and  u DEM : •  u SEF (,  ′ ) is given by the probability of sampling a single edge (, ) ∈  by a uniform distribution over all single edges such that  ′ = SEF , ().It is 0 if there is no such edge.
•  u DEM (,  ′ ) given by the probability of sampling a single edge (, ) ∈  and a double edge (, ) uniformly (i.e. by a uniform distribution over the single respectively double edges) such that  ′ = DEM , , ().It is 0 if there are no such edges.Furthermore, with  SEF ,  DEM ∈ (0, 1) and  SEF +  DEM = 1, the unrestricted transition matrix of simple moves is the convex combination of each moves transition matrix: As  is finite, these matrices are finite.
Lemma 3.5.Let  be a simple, directed graph.The unrestricted transition matrix of simple moves on ,  u SM , is doubly stochastic and irreducible on  ∞ 0 ().If  further contains at least one double edge,  u SM is aperiodic as well.
Proof.The previous lemma ensured that we do not leave  ∞ 0 ( ) with our moves.We check the remaining properties: • Aperiodicity: The aperiodicity is due to the fact that there are cycles of length 2 and length 3 (see Figure4), whose greatest common divisor is 1.
• Double Stochasticity: The move SEF , has the inverse SEF , and analogously, DEM , , has the inverse DEM , , .The probability of picking the inverse move is the same as picking the original one: the number of single respectively double edges is retained through all moves and thus is the probability of picking one, as they're selected via a uniform distribution.This makes  u SEF and  u DEM symmetric and thus doubly stochastic.The same holds true for their convex combination  u SM .• Irreducibility: We have to check that every graph  ′ ∈  ∞ 0 () is connected to  by our set of moves.This can be done by defining a standard form which is reachable by all graphs in  ∞ 0 (): let us assume some arbitrary total order on the vertex set . Then SEF and DEM could be used to flip all edges such that their ingoing vertex is ordered higher than their outgoing vertex (i.e.flip them to the upper right triangle of the adjacency matrix where possible) and push all double edges such that the tuple of outgoing vertex, ingoing vertex is minimal (i.e.push them as far top left as possible in the adjacency graph).This is a potential standard form reachable by all graphs in  ∞ 0 () using only SEF and DEM.Thus  u SM is irreducible.

Complex Moves
As we have witnessed e.g. in Figure 5 the simple moves SEF and DEM easily create and/or destroy simplices.As the formation of (especially) high-dimensional simplices by random chance is quite unlikely (Remark 2.4) we present here more complex moves to modify graphs.These combine many SEF and DEM at once while trying to retain () at the same time.
This is achieved by applying permutations of the nodes onto the edges, i.e. (, ) ↦ ((), ()).That way the number of total orders on the vertices in a clique (which correspond to simplices) is retained.However, operations of that character could change the underlying undirected graph if applied to anything else than cliques.Whereas Clique Permute applies within one clique, Clique Swap goes one step further: two equal-sized cliques swap all their edges.This not only permutes edges in two cliques, but has the potential to transfer many double edges at once, keeping the overall number of simplices highly depended on the number of double edges comparatively stable.Definition 3.7 (Clique Swap).Let  = ( , ) be a directed graph.Let  1 ,  2 ⊆  be equal-sized sets of vertices such that they form, together with their induced edge sets  1 and  2 , maximal cliques in .Let further  ∶  1 ∪  2 →  1 ∪  2 be a permutation with the additional property that ( 1 ) = ( 2 ) (and thus also ( 2 ) =  1 ).With that we define the Clique Swap where we define ((, )) = ((), ()) for ,  ∈  1 or ,  ∈  2 .
The intended purpose and easiest case of the Clique Swap is to swap faraway cliques as sketched in Figure 7.However, these two cliques might overlap, requiring the additional technicalities above.On the other hand, if  1 and  2 overlap completely, we've just got the Clique Permute.
The requirement of the cliques to be maximal is not strictly necessary.In fact, by dropping it, CS is flexible enough to define SEF, DEM and CP at the same time.In our implementation, however, we stick to maximal cliques for a variety of reasons (see Section 4) and thus define them like above.

All of the above statements hold also true for 𝐺
Proof.We show for each 1.This is proven by showing that The first straightforwardly follows from  being a bijection: if  is a bijection on , then  ×  is also a bijection on  × .
For the second we note that  ⧵ ( We define the unrestricted transition matrix which comprises all the moves introduced so far: Definition 3.9 (Unrestricted Transition Matrix).Let  = ( , ) be a directed graph.Let   be a distribution over all clique sizes occurring in .We define the unrestricted transition matrices  u CP and  u CS the following way: •  u  (,  ′ ) is defined by the product of the probabilities of picking a size  (by   ), then uniformly sampling a maximal clique  and a permutation  ∈   such that  ′ = CP   ().
•  u  (,  ′ ) is defined by the product of the probabilities of picking a size  (by   ), then uniformly sampling two maximal cliques  1 ,  2 and a bijection as described in Definition 3.7 such that  ′ = CS   1 , 2 ().
Furthermore, with ( SEF ,  DEM ,  CP ,  CS ) forming a positive probability vector, the unrestricted transition matrix  u is the convex combination of each moves transition matrix: Lemma 3.10. u is symmetric and double stochastic.
Proof.Maximal cliques are an undirected property of the graph, and as the underlying graph does not change, so neither    .Thus sampling  1 and  2 has the same chance at every time step (assuming    does not change over time).
As we picked our permutation  AB and  BA by a uniform distribution, picking the inverse permutation is just as likely.Thus the inverse move  ′ →  has the same probability as  →  ′ and the transition matrix   CS is symmetric.By definition   CS is stochastic.The same holds true for   CP , and as a convex combination of symmetric stochastic matrices also for   .A symmetric stochastic matrix is double stochastic.

Relaxing Boundaries
The simple moves introduced do not strictly retain the number of simplices (see e.g. Figure 5).Even the more complex moves only retain the number of simplices within their cliques and not in their respective neighbourhood.To avoid straying to far from our desired simplex we must restrict the transition matrix   as described in Section 2.2.3.This leads to an uncertainty if   +  − is still connected.In this section we look for relaxations necessary to connect   +  − .We define the simplex-count equivalent of Definition 2.15: Definition 3.12 (  +  − -connecting relaxed bounds).We define the simplex count bounds  −− and  ++ as   +  − -connecting with respect to an unrestricted transition matrix  u , if for all ,  ′ ∈   +  − there exists a  ∈ ℕ such that (  ++  −− )  , ′ > 0. Corollary 3.13.Let  = ( , ) be a directed graph with at least one double edge.Let  − ,  + be simplex bounds and let  −− and  ++ be   +  − -connecting bounds.Then ℳ(  ++  −− (),  u , ) forms an ergodic Markov Chain with uniform equilibrium distribution.
As we need to subsequently filter the output for membership in   +  − , we are interested in finding tight boundaries, meaning that the complement of   +  − with respect to the irreducible component of   ++  −− that contains   +  − should not be unnecessarily large.

Single-Edge-Only Case
We first look at a special case: graphs that do not have any double edges.This is a property retained throughout all graphs with the same underlying undirected graph and the same number of edges and thus the moves introduced respect this property as well.
In this case we propose the following bounds: Conjecture 3.14.Let  be a directed graph without double edges and  − ,  + simplex count boundaries.Then Assuming that we are looking for graphs with a simplex count far above the average, this is an easy and tight relaxed boundary.To convince ourselves of the tightness, we remind ourselves of the distribution of simplex counts as in Figure 1 (right) and Remark 2.4: While  ++ = ∞ is probably not optimal (could be lower), it does not really matter, as with rising simplex count the number of graphs quickly decay, making it unlikely to sample them with a uniform sampler.
We now outlay our attempts to prove Conjecture 3.14: For easier analysis, we only take SEF into account.
The core idea is to connect ,  ′ through star points which are directed acyclic graphs (DAGs).Let us assume there are DAGs  and  ′ connected to  (respectively to  ′ ).Now  and  ′ do not have any cycles and thus each correspond to a (global!) order on .Thus both are connected by a series of SEF, as each SEF corresponds to a neighbouring transposition which are known to be generators of the symmetric group of permutations.
This leaves the question if and how  is connected to its corresponding DAG .This is equivalent to the question if any  can be transformed to a DAG by a series of SEF, where each individual step may not decrease  2 (every obstruction to a simplex in a clique is perceivable as a 3-cycle so eradicating them is enough).
In the case of  being fully connected (i.e. a clique), we give, with Theorem A.3, a proof to this statement.The general case of non-complete graphs turns out to be surprisingly intricate to prove3 : despite significant effort, we have to state this as the very elementary Conjecture A.4 for now.Yet, both Conjecture 3.14 and A.4 seem to be correct, as a computational search for counterexamples was unsuccessful so far, despite having checked billions of graphs.
Both Theorem A.3 and Conjecture A.4 also make a statement about the number of SEFs required: with them, we can bound the graph diameter of (  ++  −− ,  ) to || ≤ diam((  ++  −− ,  )) ≤ 2||, a potentially useful statement for estimating mixing times.

Bounds for Graphs with Double Edges
In the previous subsection we have discussed suitable relaxations in the case of single-edge-only graphs.There the number of -cliques was a natural and rather strict upper bound to connect everything.This analysis is not enough to tackle the graphs motivating this work, biological connectomes.Indeed, C. Elegans contains 233 ( 11%) double edges and the Layers 1-4 of the BBP reconstruction contain 15258 ( 1%) double edges.Simultaneously, only a few double edges may result in a combinatorial explosion of the number of high dimensional directed simplices: an 8-clique in C. Elegans requires only ( 82 ) = 28 double edges to generate 8! = 40320 7-simplices.This questions the choice of the open upper boundary used in the single-edge-case: while it is still ( − ,  + )-connecting, it is no longer a strict and thus efficient choice.An experiment confirmed that sampling from C. Elegans with an open upper boundary is not practical.
Thus we approach the double-edge-case differently.A rigorous treatment similar to the single-edge-only-case is currently out of scope, so we describe assumptions used in the implementation.They are based on the observation that the graphs investigated are quite large, yet somewhat regular and random.
• The graphs are large enough to have a considerable buffering-capacity: small changes in the simplex count caused by modifications in some parts of the graph can be (approximately) compensated in other, distant parts of the graph.It is then enough to ensure that boundaries are relaxed enough to allow for a single DEM as the graph may use its buffering capacity afterwards to return to the former number of overall simplices.This buffering capacity enables more complex local changes which would otherwise result in a (temporary) boundary violation if performed in direct succession.
• We can treat the simplex counts separately for each dimension: there are, for every dimension, enough maximal cliques available such that the lower dimensional effects of high dimensional changes can be compensated using solely the maximal cliques of these lower dimensions.But as these cliques are maximal the change will not propagate up again.
• In particular, there are many maximal 1-cliques which allow removing or adding double edges from the remaining graph by "hiding" them in these 1-cliques.This allows us to (locally) view DEM as removing or adding an edge.
With these assumptions we have to look for the worst-case scenario of connecting two graphs with similar simplex counts.Our approximation of a worstcase scenario is sketched in Figure 8: two -sized cliques ( being maximal clique size), one (upper left) full with double edges and the other (lower left) void of any double edges.The upper clique thus exhibits ! many simplices and the lower at most one.In the transformed graph (right) the double edges have switched their position.This is exactly the scenario CS would resolve in a single transition and was the original motivation to include this transition in the first place.CS entails modifying multiple edges simultaneously though, a violating critical to the buffering effect.Thus, for the theoretical connectivity analysis, we only consider DEM for the rest of this chapter.
Here we identify the following high-impact DEM.Removing the first edge of a clique full of double edges halves the number simplices in this clique (and analogously, adding the last one missing doubles the number of simplices).This effect is bounded primarily by the size of the clique, leading to Δ (1)  = 1 2 ( + 1)! as a first approximation of a suitable boundary.
Usually, especially in the higher dimensions, ( + 1)! is (much) bigger than  +  , meaning that cliques full of double edges do not appear in the graphs desired.In these cases Δ is then an impractically generous boundary.A tighter one requires the integer sequence A058298 which describes in the th entry the maximal number of -simplices in a large enough clique with  < ( + 1)! double edges (see Remark 2.5).In the adapted situation of Figure 8 the maximal unavoidable simplex count loss or gain in a series of DEM is then (generously) bounded by Δ It is a rather generous estimate as there are many cases where smaller changes in the clique are possible, e.g. by internally moving double edges instead of moving them out of the clique.This handles the scenario depicted in Figure 8 and similar ones.This scenario is a bit optimistic though, as one edge may participate in many cliques and thus simplices at the same time.This becomes a problem if a series of transition requires modifying an edge with a strong impact on the simplex count in multiple cliques and there is no alternative series of transitions where this situation could be avoided or defused first.In the highest dimensions, where in practical applications double edges are rather rare, this should not happen.In the lower dimensions this is quite likely however, thus we approximate this network effect by assuming that the lower dimensional simplices are embedded a fully-double-edge highest-dimensional -clique.With this assumption, there are ( −2 −2 ) cliques of size  which all share the same (undirected) edge, a multiplier applied to our previous boundary relaxations.The dimension-wise minimum of Δ (1) and Δ (2) adjusted by this clustering combinatorial factor should then be wide enough to connect   +  − .We add it the upper end of target boundaries: , Δ ) .

Implementation and Mixing Times
In this section we adopt the perspective of random walks instead of distributions.We start with the rough sketch of the algorithm used for sampling and then discuss the experimental setup used to establish reasonable mixing time.
line supplied integer serving as a seed for the pseudo random number generator is set.

Calculating a Transition
At every step, we first choose the kind of transition to apply.The distribution can be selected by a command line option and defaults to 10% SEF, 10% DEM, 60% CP, and 20% CS, values that seemed reasonable but where never optimized.
A separate function for each kind of transition returns a list of directed edges to be deleted or inserted.
• Single Edge Flip: A SEF is determined by a uniformly sampled single edge (, ).We return (, ) marked for deletion and (, ) marked for insertion.
Here we can safely assume that not all edges are double edges, as there would be only one such graph.
• Double Edge Move: The directed single edge and the undirected double edge are sampled uniformly together with a bit indicating the direction to remove from the double edge.When there is no double edge, we resort to a transition that does nothing.
• Clique Permute: We only consider maximal cliques of the underlying undirected graph.Instead of sampling them uniformly, we first decide the size, weighted by the fifth root of the number of maximal cliques with this size in the graph to increase the chance of larger cliques.Then, for the selected size, a maximal clique is drawn uniformly.The permutation is then chosen uniformly.
• Clique Swap: Two equally sized cliques (with vertex sets  and ) are chosen as above, and the bijection required by Definition 3.7 is constructed from three uniformly chosen permutations: One for the intersection of  and , the second for the bijection from  ⧵  to  ⧵ , and the third for the bijection from  ⧵  to  ⧵ .It is possible that the same cliques are chosen, in which case this forms a Clique Permute.

Applying the Transition
With the list of changes generated we then try the transition and check if the simplex count is still inside the given bounds.If not, the transition is reversed.
Mutating the graph in-place avoids copying it, as the graph is typically much larger than the transition changeset.Instead of counting (and therefore enumerating, a computationally costly process) all simplices after each transition we calculate the effect on the simplex count in the neighbourhood nbhd of the affected edges (see Definition 2.8).We then calculate the change of the simplex count by restricting the graph to nbhd.It follows from Lemma 2.9 that ( ′ ) = () + ( ′ | nbhd ) − (| nbhd ).
Exploiting this locality property is key to making the transition step fast and thus MCMC-sampling computationally feasible.

Initialization and Data Structures
To make the operations stated above efficient we precompute and cache properties only dependent on the underlying undirected graph, as it is invariant with respect to the transitions: we search for the maximal undirected cliques via the Bron-Kerbosch algorithm [21], and save them as arrays of vertex indices.
We also build a hash map from the undirected edges (sorted pairs of vertex indices) to their neighbourhoods saved as an array of vertices according to Definition 2.8.This speeds up calculating the neighbourhood of the transitions.
The current graph is represented both as a dense and a sparse adjacency matrix.The sparse adjacency matrix is an indexable hash set of edges (integer pairs).It provides (1) edge insertion, deletion, and uniform sampling.Just like flagser, we maintain a dense matrix (| | 2 bits) as it is faster for read only lookups in practice.To uniformly sample double edges quickly, another indexable hash set is maintained.This combined graph representation is implemented in the EdgeMapGraph structure in the crate flag-complex.

Connecting Boundaries and Sampling Distance
The values for  −− ,  ++ are calculated as described in Section 3.3 but can easily be hardcoded as well.The sampling distance is command-line and defaults to 2 1 log 2  1 , where  1 describes the number of directed edges.This is a reasonable lower bound on sampling distances (compare for mixing times for random walks on hypercubes, see [8]).We dedicate the next section to validating these sampling distances in computational experiments.

Sampling Distance
A rigorous theoretical, yet practically useful result on the mixing time is out of scope for this paper, as it would require extensive theoretical insight into the geometry of the space of directed graphs and simplex counts.Still, the question of what sampling distance to choose is crucial in practice.In this section we describe and conduct two experiments to estimate reasonable defaults used in our implementation and for the experiments described in Section 5: analyzing the distance between two graphs (here the Hamming distance on the adjacency matrix) over time and a test to check for independence of samples.

Hamming Distance over Time
Definition 4.1 (Hamming distance).Let  1 = ( ,  1 ),  2 = ( ,  2 ) be two graphs over the same vertex set.The Hamming distance is defined by the number of edges unique to one of the graphs This definition is equal to the Hamming distance of the bit strings given by the graphs adjacency matrices.
We measure the development of the distance with respect to the number of Markov Chain transitions  ∈ {0, ⋯ ,  }, where  ∈ ℕ is the maximal temporal distance considered.To avoid measuring only the distance to the original starting graph, we sample one very long chain (893) and investigate not only from the starting graph, but also from the offsets  0 ∈ {0,  , 2 , ⋯}.The behaviour of the distance over time is then described by  ↦   (  0 ,   0 + ).
For Figure 9 we started the sampler with the connectome of C. Elegans, target bounds of ±1% and default relaxed bounds and seed.We plot the range of observed values (blue, shaded) and mean (orange line) over all offsets  0 ∈ {0, ⋯ , 892}.The behaviour with  0 = 0 (starting graph is C. Elegans) is plotted explicitly (blue line).The red horizontal line is the (experimental) limit behaviour of the Hamming distance.A quick visual analysis seems to suggest that the mean is sufficiently close to the limit after roughly 50000 ≃ 2 1 log  1 steps.
We further tested the sampler restricted to simple moves as defined in Definition 3.2 (see green line in Figure 9).Initially simple moves move away from the initial state much faster, as they exhibit roughly half the rejection rate compared to complex moves and thus resample less often.But ultimately the green line can be observed to struggle to converge to the limit.We interpret that as a justification of the introduction of the complex moves, even though they take on average 8 times longer to compute per step.
We ran a similar analysis on the BBP layer 1-3 connectome with similar results (not shown).

𝜒 2 -Test for Independence
Assume the Markov Chain used by the sampler yields the sequence of graphs   with  ∈ ℕ 0 .We want to find the number of steps  required between two graphs   and  + such that we cannot reject the hypothesis "  and  + are stochastically independent".
We check for multiple sampling distances  which are then fixed for each hypothesis test, using a binarizing functions  ∶ DiGraph → {0, 1} to derive bit sequences   = (  ) with  ∈ ℕ 0 .We now count the pairs ( 2 ,  2+1 ) ∈ {00, 01, 10, 11} of the bit sequences and use Pearson's  2 test for independence [3]: if they are independent, the frequency of the pairs should match the frequencies expected by counting single bits.Thus the final test hypothesis is " 2 and  2+1 are stochastically independent for all applicable ".
As binary binning functions we split   (  ) by the median for each dimension .
The test results are shown in Figure 10.We show the values of the test statistic while varying the sampling distance  and the simplex dimension  used for binning.For C. Elegans the data consists of 704964 graphs sampled with distance 100, and we must reject the independence hypothesis for sampling distances smaller than 40000.For BBP layers 1-3 the data consists of 108846 graphs sampled with distance 10000, and we must reject the independence hypothesis for sampling distances smaller than 6 ⋅ 10 6 .This test confirms, both for C. Elegans and BBP layer 1-3, the choice of sampling distances  ≃ 2 1 log 2 ( 1 ) presented in Section 4.1.5.

Homology Analysis for C. Elegans
C. Elegans is a well studied specimen in neuroscience.Its neural wiring is determined completely by its genome and thus it has been studied extensively, leading to an exhaustive map of all its neurons and synapses.We explore the data in the version described in [20], with 279 neurons and 2194 directed synapses.
For our connectomic analysis we ignore any additional information associated with the neurons and synapses (like type, size or location in 3-space) and abstract everything away which cannot be represented in a directed graph in the following way: A neuron forms a vertex and a directed chemical synapse with  as its presynaptic and  as its postsynaptic neurons becomes an edge (, ) in a graph henceforth referred to as .
Using flagser [6,11] we calculate the simplex count and Betti numbers, with same results as in Figure 1, left.As already discussed in that study, the simplex count and Betti numbers are very distinct to the ones found in comparable Erdős-Rényi graphs.
We pursue a target-relaxation of ±1% and use boundaries as described in Section 3.3. − 2408 samples remained.This low yield is due to the comparatively small number of 7-simplices in C.Elegans, where ±1% leaves a very small margin for error.
On the remaining 2408 graphs we analyzed the homology using flagser.The result is pretty clear:  C.Elegans exhibits significantly lower Betti numbers in low dimensions (1 to 3), but higher Betti numbers in high dimensions (6 and 7) than expected on the average of the null model.Note to the Reviewer: We will make the sampled, filtered graphs available on the long-term scientific data hosting service of the TU Graz.We still need to figure out some details, a footnote with the URL will be added in time.

Homology Analysis for Layers 1-4 of a Rodent Somatosensory Cortex
Computationally, C. Elegans is nowadays seen more as a toy example.To demonstrate that the MCMC-approach is feasible for significantly bigger connectomes as well, we perform the same analysis on the connectome of layers 1-4 of the Blue Brain Project statistical reconstruction of the somatosensory cortex of a rat.With 12519 neurons and 1431930 directed synapses it about half (vertices) to a quarter (edges) of the full reconstruction.We invested the computational power of 32 AMD EPYC 7543 cores over roughly 24 days to sample around 300 samples per core, each sample differs by 114554400 transitions, twice the sampling distance recommended in Section 4.2 to make up for the low transition acceptance rate of around 25%.After rejection-filtering these 9600 samples, 407 remained.We present here the Betti numbers of the original graph as well as mean and standard deviation of the graphs of the null model: Interestingly, the BBP graph exhibits significantly less holes in every dimension except  = 2.We can therefore conclude that the geometric spaces described by the flag complex are remarkably flat for graphs with a comparable number of simplices.We hypothesize that this is due to the strong overall flow of information: The majority of synapses go from lower levels to the higher ones.Note to the Reviewer: Links to samples missing, see above.

Summary and Outlook
Summary The preceding sections began with establishing the need for more refined null models than ER-graphs when investigating directed flag complexes.We decided to first investigate a very specialized null model -it consists of directed graphs with a specified range of simplices in the flag complex and keeps the number of edges and the underlying undirected graph.We then designed, implemented and tested an MCMC-based sampler and made it publicly available.In the special case of oriented graphs (i.e.directed graphs without double edges) we state an elementary graph-theoretical conjecture which, if/once proven, implies correctness and efficiency of our sampler.An investigation of the statistical properties of the homology (i.e.Betti numbers) of two connectomic graphs well-known in Neuroscience then followed and we discovered that they, with some significance, indeed behave extraordinary.

Implications
We hope that this work leads to a more accurate and quicker assessment of the worth of TDA methods -not only Betti-numbers, but connectivity and simplicial completions probability -and by that the application of TDA to connectomes and thus Neuroscience in general.In particular, now that we know that a variety of connectomes exhibits a topology significantly nonrandom, we should be encouraged to invest more resources into investigating the how and why of the emergence of holes (or the suspicious lack thereof).We furthermore hope that the insights, approaches and code are of help in developing other models in the area of directed and undirected flag complexes as well as directed (non-flag) simplicial complexes.
Limitations From a more applied perspective interested in the graph analysis the null model presented here might be too narrow.By restricting the underlying undirected graph and the number of edges we only allow freedom in the direction of edges and, in the case of non-oriented graphs, the placement of double edges.This corresponds more to concepts of "flow" or "direction".In contrast, the actual "connectivity" is not changed and thus no statements concerning it may be derived from this model.Simultaneously, the model presented here is too loose for a more rigorous analysis: unlike more traditional approaches like the graph/simplicial configuration model we do not maintain any notion of directed vertex or edge participation (the undirected one is trivially maintained due to the fixed undirected backbone).Thus our model does not respect simplex clustering (or a lack thereof) as well as the existing model for undirected simplicial complexes or the related graph configuration model.It could very well be that a null model taking vertex participation into account would not yield comparable results to the one investigated here.
From a computational standpoint there is room for improvement as well: not only very small graphs, where  − ,  + becomes really narrow, and large graphs, due to the estimated time of  ⋅  1 log( 1 ) steps, but also (locally) dense graphs where single transitions become very expensive due to the large neighbourhood require substantial computational effort to sample.Even with state-of-the-art CPUs a useful number of samples requires weeks to months of computation time.The most promising and reachable lever is decreasing the rejection rate, i.e. finding better  ++ ,  −− -bounds.
Both the quest for better boundaries while maintaining irreducibility as well as restricting the model further to respect vertex participation require a substantial better insight into the geometry of   ++  −− in the double-edge case.Incidentally, this is also the biggest theoretical weakness right now: so far we do not have rigorous mathematical proofs of irreducibility of the transition matrices.Note that this does not rob this approach of all value -even if the transition matrices used were not irreducible, the statement that the connectomes investigated have unusual topology remains true, although only w.r.t. the connected component.Nonetheless, it would be satisfying and useful to be able to formalize, and thus understand, the double-edge case (see Section 3.3.2) as well or even better as the single-edge-case is already understood.

Potential Research Directions
As discussed in Section 3.3.1 and Appendix A, the question of connectivity of  ∞  −− for oriented graphs can be stated in a innocuous-looking graph-theoretical way.Proving or refuting Conjecture A.4 is thus an obvious open problem.More generally, the geometry of that space is probably rather tame and worth investigating.Thorough insights would probably enable formal descriptions of mixing time and rejection rate.
The double-edge-case seems considerably more challenging.Simultaneously it is more rewarding, as all practical applications found so for fall into that category.Formalizing, sharpening and proving the hand-wavy assumptions and rules given in Section 3.3.2 is out of scope for the current paper, but should be tackled as soon as possible.
In practical applications though, most of the rejections are due to violations in the highest dimension.From that perspective refining the rules of thumb would be useful enough already.In particular, utilizing knowledge about the undirected graph (i.e. the clique complex, which is fixed!) might be enough to increase the samplers yield by several orders of magnitude.
Once the yield is increased and the necessary computational power is acquired the analysis performed in Section 5 can be applied to bigger and denser connectomes.In particular, the full BBP connectome and the larvae of Drosophila [22] should be addressed.
Similarly we would like to test if -connectedness and the simplex completion probability [16] are as able as Betti-number-analysis in picking up the particularities of actual connectomes w.r.t. the null model here.
Furthermore we hope to combine this work with an (upcoming) one on sampling undirected flag complexes while maintaining clique count (and potentially vertex participation) to construct a sampler which samples directed flag complexes and does not respect the underlying undirected graph.This would form a null model less narrow and artificially restricted.
Proof.The proof works by two loops: The outer loop is an induction over the number of vertices , as with  = 1 the statement is trivially true.An inner loop reduces the indegree of a vertex : starting with a (potentially non-unique) vertex  of minimal indegree, we flip incoming edges such that they become outgoing edges to further reduce indeg() until it becomes 0. The vertex  is then a source and the subproblem of ⧵{} is considered, forming the inductions step  ↦  −1 of the outer loop.The step where  is removed constructs a total order on : the elements removed earlier (which only have outgoing connections to the remaining graph) are considered to be smaller.This results in a DAG.
It remains to verify that this indegree-reducing SEF is legal, i.e. it strictly reduces the number of 3-cycles in the graph.Let  ∈  with indeg() minimal, i.e. indeg() ≤ indeg() for all  ∈ .By assumption indeg() > 0, otherwise  would already be a source and we could consider the strictly smaller subproblem of  ⧵ {}.Let  = In() be the set of all incoming vertices and  = Out() the set of all outgoing vertices of .Now consider  ∈ : Flipping the edge (, ) creates as many 3-cycles as outdeg  (), the outdegree of  with respect to , and destroys as many 3-cycles as indeg  (), the indegree of  with respect to .Thus a modification of  via a flip of (, ) strictly decreases the number of 3-cycles and thus adheres to our restriction iff indeg  () > outdeg  ().We observe (justifications below): As  is complete,  = {} ⊍  ⊍  and thus indeg() = indeg  () + indeg  (), leading to the first line.Another implication of completeness is the relation of indegree and outdegree: as, between to different vertices, there must be exactly one edge in one direction or another, indeg() = | | − 1 − outdeg().This statement still holds true when restricted to , thus the identity indeg  () = || − 1 − outdeg() leads to the second line.The identity || = indeg() leads, together with the initial assumption that  has minimal indegree, to the lower bound proposed in the third line.This proves conformity of the inner loop of our requirement to reduce the number of 3-cycles each step.
Note that once an edge (, ) has been flipped, it will not be flipped again:  will be the unique element of minimal degree in the next step, and  lies now in Out().Once  has become a source, all connections with  will no longer be considered.Thus each edge must be flipped at most once.Furthermore, an element with minimal indegree may have at most  2 incoming edges, leading to the additional factor of 1  2 in the number of steps required.Does this Theorem also hold in the case of non-complete graphs?No.If one weakens the requirement of a strictly monotonic decrease in the number of 3-cycles to just a monotonic decrease in the number of 3-cycles (i.e each SEF may not increase the number of 3-cycles) then the result still seems to be true: The graph can be transformed into a DAG while flipping each edge at most once.
We tested this hypothesis on ER-graphs configured with ( = 30,  = 0.8), ( = 60,  = 0.4) and ( = 300,  = 0.1) on 150 million graphs each and did not find any counterexample.However, despite ongoing effort, the following conjecture still remains unproven: Conjecture A.4.A finite simple oriented graph  = ( , ) of  vertices can be transformed into a directed acyclic graph by at most || applications of SEF such that the number of 3-cycles does not increase with each individual step.

Figure 1 :
Figure 1: Left: Simplex Count and Betti numbers of C. Elegans compared to an ER-graph.Figure taken from [15].Right: Histogram of the simplex count in dimension 2 of random ER-graphs with the same underlying undirected graph as C. Elegans.

Figure 2 :
Figure 2: A sketch of a restricted random walk.Dots on the grid represent states in .The green area represent the desired states in , the yellow area represents a potential -connecting subset .The red area represents   and is forbidden: Instead of entering the red area (red arrow) we immediately resample (blue arrow).

Figure 4 :
Figure 4: There is a cycle of length 3 (DEM from left to middle, DEM to the right, then SEF from right to left) and a cycle of length 2 (SEF from left to right and then immediately back).

Figure 7 :
Figure 7: Clique Swap: On the left hand side we observe two 5-cliques with differing number of double edges.Clique Swap now swaps (via the bijection (0, 1, 2, 3, 4) ↦(D,C,B,E,A)) all the edges from the upper clique to the lower clique and vice versa.The remaining graph and their edges to these cliques remain untouched.

Figure 8 :
Figure 8: A potential worst-case scenario to find relaxed simplex bounds for.

Figure 10 :
Figure 10:  2 Test statistic  , using the bit sequences derived from simplex counts in different dimensions.Thin, colored: Values from individual tests.Thick, same color: smoothed.Black, dashed: Acceptance limit  2 1− for  = 1%.If the samples are independent, the probability of   <  2 1− is 99%.The dimensions shown here are exemplary, including the longest rejecting ones.Left: Starting with C. Elegans.Right: Starting with BBP layers 1-3.

Figure 11 :
Figure 11: Histogram on the value of the Betti numbers in the null model in dimension 1 and 3 (left) and 6 and 7 (right).Vertical lines show respective values of C.Elegans.

Figure 12 :
Figure 12: Histogram on the value of the Betti numbers in the null model in dimensions 1 (left), 2 (middle), and 3 (right).Vertical lines show respective values of BBP, layers 1-4.

Definition 2.2 (Cliques
• Let  = (  ,   ) and  = ( , ) be (directed or undirected) graphs with   ⊆  and   ⊆ .Then  is called a subgraph of .Furthermore, if   is maximal w.r.t   and , then  is called the induced subgraph and is denoted with [].and Simplices).Let  = ( , ) be a directed graph.• Let   ⊆  be a vertex set such that for all ,  ∈   ,  ≠  there is (, ) ∈  or (, ) ∈ .Then  = [  ] is called a clique of order  = |  |, or an -clique. is a subgraph  = (  ,   ) of a clique together with a total order ≺ on   represented by its edges, i.e. edges point from lower to higher vertices w.r.t.≺:   = {(, ) ∈  | ,  ∈   and  ≺ }.With  + 1 = |  | we refer to  as the dimension of  and call  a simplex.A simplex is often denoted by giving the vertices in rising order:  = [ 0 , …   ] implies that that   ≺   iff  < .
4. This is clear by the two previous statements.4.A simplex is a clique together with directed edges corresponding to a local total order.As the undirected graph is retained, we must check only for the number of orders present.Let  = ( 0 ,  1 , … ,   ) an -simplex in [ 1 ] ∪ [ 2 ]and  be a permutation with the required properties.Then the simplex  is transformed into the simplex () = (( 0 ), ( 1 ), … , (  )), showing that no simplex is lost by CS.A similar argument with  −1 shows that no simplex is gained as well.Thus the number of simplices is retained.5.CP  is actually a special case of CS where  1 =  2 .Above results then follow immediately.
2, with the addition of  ++ 8 = 10 as  has cliques of size 9 but no simplices of dimension 8. Thus the simplex counts and boundaries are:  respects the boundaries  −− ,  ++ by resampling instead (around 40% of transitions were accepted) and chooses its moves SEF, DEM, CP, CS with the probabilities [0.1, 0.1, 0.6, 0.2].For the sampling distance 48704 ≃ 2 1  2 ( 1 ) was chosen as motivated in Section 4.2.Investing 8 days of computation time using a 32 core AMD EPYC 7543 CPU we sampled 10000 samples of   ++  −− for each of the 128 starting seeds.After rejection everything not in   +