Cliques of Neurons Bound into Cavities Provide a Missing Link between Structure and Function

The lack of a formal link between neural network structure and its emergent function has hampered our understanding of how the brain processes information. We have now come closer to describing such a link by taking the direction of synaptic transmission into account, constructing graphs of a network that reflect the direction of information flow, and analyzing these directed graphs using algebraic topology. Applying this approach to a local network of neurons in the neocortex revealed a remarkably intricate and previously unseen topology of synaptic connectivity. The synaptic network contains an abundance of cliques of neurons bound into cavities that guide the emergence of correlated activity. In response to stimuli, correlated activity binds synaptically connected neurons into functional cliques and cavities that evolve in a stereotypical sequence toward peak complexity. We propose that the brain processes stimuli by forming increasingly complex functional cliques and cavities.

The Blue Brain Project (BBP) has recently generated the first draft digital reconstruction and simulation of a microcircuit of neurons in the neocortex of a two-week-old rat (Figure 1A) [7].This reconstruction is made available through the Neocortical Microcircuit Portal (https://bbpnmc.epfl.ch)[11].Based on sparse anatomical and physiological data for neurons and synapses and on a variety of biologically motivated organizing principles, the complete connectivity between neurons belonging to a neocortical microcircuit was digitally reconstructed -a microconnectome.The structural properties of the reconstruction have been extensively validated against independent data, and simulations of the reconstruction reproduced multiple in vitro and in vivo experiments without adjusting any parameter, further validating its biological accuracy.
In this article we apply methods from topology to the analysis of 42 variants of the digital reconstruction, grouped in six sets of seven microciruits each.The first five sets of microcircuits take into account biological variability in layer heights, proportions of cell types, and cell densities from five individual rats, while the sixth set is based on the average reconstruction across the five individuals.To form each set of microcircuits, seven statistically varying instantiations of the microcircuit were reconstructed [12].The 42 microcircuits are therefore all distinct, though the degree of resemblance within each set is higher than that between sets.The structural connectivity of each reconstructed microcircuit can be represented as a directed graph with approximately 3 × 10 4 vertices and 8 × 10 6 edges, while its functional connectivity can be represented as a time series of subgraphs formed by functionally effective connections.
Our topological analysis of the detailed structural and functional connectivity of these 42 neural microcircuits led to a number of surprising observations.Firstly, we found that the distribution of directed cliques (directed all-to-all connected subsets) of neurons by size is highly significantly different from both that in Erdős-Rényi random graphs with the same number of vertices and the same average connection probability and that in more sophisticated random graphs, constructed either by taking into account distance-dependent probabilities varying within and between cortical layers or morphological types of neurons, or according to Peters' Rule [9], [10] (Figure 1D).In particular, we found that directed cliques of up to eight neurons are highly prominent motifs in the reconstructed microcircuits: the average microcircuit incorporates approximately 10 8 3-cliques and 4-cliques, approximately 10 7 5-cliques, approximately 10 5 6-cliques, and approximately 10 3 7-cliques.Taking the alternating sum of the numbers of directed cliques of various sizes, we computed the Euler characteristic (EC) [5] of the 42 reconstructed microcircuits, obtaining in each case a value on the order of 10 7 , indicating a preponderence of directed cliques consisting of an odd number of neurons (Figure 2).
Another topological metric that we considered in this analysis are the Betti numbers (SI, Supplementary Text, ST1.3) associated to a graph via its directed flag complex (Figure 1C).These are a sequence of natural numbers β 0 , β 1 , β 2 , ... that measure the higher-order organizational complexity of the network, detecting "cyclic" chains of intersecting directed cliques.For each graph considered here we determined its homological dimension, i.e., the maximum n such that β n = 0. We showed that the reconstructed microcircuits have homological dimension 5 (Figure 2D), whereas the random graphs considered have homological dimension at most 4, strongly indicating that the microcircuits possess a higher degree of organizational complexity than the random graphs.
Topological methods also enabled us to distinguish functional responses to different input patterns fed into the microcircuit through thalamo-cortical connections.We ran simulations of neural activity in one of the reconstructed microcircuits during one second, over the course of which a given stimulus was applied every 50 ms (Figure 3).We then binned the output of the simulations by 5 ms timesteps and associated to each timestep a transmission-response graph, the vertices of which are all of the neurons in the microcircuit and the edges of which encode connections in the microcircuit whose activity in that time step leads to firing of the postsynaptic neuron (Figure 4).The size of the time bins and the precise rule for formation of the transmission-response graph for each time bin are biologically motivated, as explained in more detail in the Supplementary Methods section (SI, Supplementary Methods, SM1).
From the time series of transmission-response graphs for each of 20 trials of two different stimuli (called Circle and Point for geometric reasons (Figure 4A), we derived time series of two non-topological metrics (mean firing rate and number of edges in the transmission-response graph) and five topological metrics (the number of 3-cliques, EC, β 0 , β 1 , and β 2 ) and applied a Gaussian Bayes classifier (SI, Supplementary Methods, SM2) to determine how successfully each of the metrics classified the 40 trials in the time bins corresponding to the first two stimulations and in the time bins immediately following those stimulations (Figure 5).In each of those crucial time bins, the metrics that were most successful at classification the number of 3-cliques (denoted 2D in the figure), β 2 , and, in one case, the Euler characteristic (Figure 2A).
We expect the methods applied here will prove useful for studying networks in general.

Structural topology
We computed the binary adjacency matrices of all 42 digitally reconstructed microcircuits and then generated the associated directed flag complexes (SI, Supplementary Text, ST1.2), which are oriented simplicial complexes encoding the connectivity of all orders of the underlying directed graph: to each directed nclique (SI, Supplementary Text, ST1.2) in the underlying graph corresponds to an oriented (n − 1)-simplex in the flag complex, and the faces of a simplex correspond to the directed subcliques of its associated directed clique (Figure 1 B and C).For each neuron in the microcircuit, there is a vertex in the underlying directed graph that is labelled with the unique global identification number (GID) of the neuron.The (j, k)-coefficient of the structural adjacency matrix is 1 if and only if there is a directed connection in the microcircuit from the neuron with GID j to the neuron with GID k.We refer to this adjacency matrix as the structural matrix of the microcircuit and to its associated directed flag complex as a neocortical microcircuit complex or N-complex .
Having computed each of the 42 N-complexes, we counted the simplices in each dimension.For comparison with non-biological matrices, we generated five Erdős-Rényi random graphs [4] of a comparable size (31,000 vertices) and connection probability 0.8%, the same as the average arising from the structural matrices of the microcircuits (SI, Supplementary Methods SM3.1).We refer to the associated directed flag complexes as ER-complexes.
To have a more biological control, we also generated 20 adjacency matrices, given by partly randomizing the structural matrix of one of the average microcircuits, taking into account its biologically meaningful division into six layers in 10 cases and into 55 morphological neuron types (m-types) [7] in 10 cases.The randomization was carried out so that the distance-dependent connection probability for all pairs of layers (respectively, pairs of m-types) was identical to that of the original matrix, i.e., for each pair of layers (respectively, m-types) the number of connections between them was the same as that of the original and for each 25 µm distance bin the number of connections was identical.The matrices are completely random otherwise (SI, Supplementary Methods SM3.2, SM3.3).We call the associated directed flag complexes L-complexes (respectively, M-complexes).Note that since each m-type is restricted to a fixed layer, the M-complex should retain more of the structure of the original N-complex than the L-complex.Our final and most biological control consisted in the generation of 10 connectivity matrices for 31,000 neurons according to Peters' Rule [9], [10] for which the associated directed flag complexes are called P-complexes (SI, Supplementary Methods SM3.4).Having carried out the computations for 10 control matrices out of each randomized set of 20, the very small variance in the results convinced us that no further computations should be needed.
The resulting distribution of simplices displayed highly consistent behavior among the N-complexes, all of which we computed, with a small variation among the samples arising from different rats.Note that Figure 1 represents the analysis only of the seven N-complexes arising from the average reconstruction because the randomizations are based on those microcircuits.The ER-complexes showed almost identical behavior among the different instances, as did the L-complexes, M-complexes, and P-complexes.On the other hand, the N-complexes exhibited remarkably different distributions from the various random complexes (Figure 1 D), with much greater numbers of simplices and simplices of significantly higher dimension.We computed the Euler characteristic of all N-complexes, as well as that of the various random complexes, obtaining large positive values in all cases, due to the predominance of even-dimensional (particularly 2-dimensional) simplices.
The Betti numbers (SI, Supplementary Text, ST1.2) of a simplicial complex provide a much finer and more sophisticated measure of its organizational complexity than the dimension-wise simplex count or the Euler characteristic.The n-th Betti number, β n , counts the number of chains of simplices intersecting along faces to create an "n-dimensional hole" in the complex, which requires a certain degree of organization among the simplices.On the other hand, computation of the Betti numbers is much more expensive than that of the directed flag complex of a directed graph or its Euler characteristic.In fact, the sheer size of the complexes we considered here made it practically impossible to do so on a computer with 256 GB of RAM.We succeeded in computing the highest nonzero Betti numbers of the Ncomplexes, however, by restricting our attention to the 5-th and 6-th coskeleta (SI, Supplementary Text, ST1.2).The top Betti number in all N-complexes appeared in dimension 5, with β 5 varying between 1 and 80 (Figure 2D).By contrast, β n = 0 for all n > 3 for all ER-complexes and P-complexes considered, while β n = 0 for all n > 4 for all L-complexes and M-complexes.Moreover β 4 varies between 0 and 6 for all L-complexes and M-complexes, so that these Betti numbers are almost negligible.

Functional topology
We tested our methods on active microcircuits as well.In an experiment that we call the Point vs. Circle test, we activated in a simulation the incoming thalamocortical fibers of one of the average that the stimulated fibers formed first a point shape, then a circle shape [7].The size of the point shape was chosen such that the average firing rate of the neurons was essentially the same as for the circle shape, and in both cases the fibers were activated regularly and synchronously with a frequency of 20 Hz for one second, similar to the whisker deflection approximation in [7,Figure 17A].We performed 20 trials of each stimulus (Figure 3).The trials of each stimulus exhibit biological trial-to-trial variability in the neural response, due to the stochasticity of the synapse models and of some of the ion channel models.The aim of this experiment was to determine whether our topological methods were able to classify the two different stimuli, the point and the circle better than the firing rate, which is largely overlapping for the first two stimulations (see Figure 3D).
After a systematic analysis to determine the appropriate time bin size and conditions for probable spike transmission from one neuron to another (SI, Supplementary Methods, SM1.4), we divided the activity of the microcircuit into 5 ms time bins for 1 second after the initial stimulation and recorded for each 0 ≤ n < 200 a functional connectivity matrix A(n) for the times between 5n ms and 5(n + 1) ms.The (j, k)-coefficient of the binary matrix A(n) is 1 if and only if the following three conditions are satisfied, where s j i denotes the time of the i-th spike of neuron j.
(1) The (j, k)-coefficient of the structural matrix is 1, i.e., there is a structural connection from the neuron with GID j to the neuron with GID k.
(2) There is some i such that 5n ms ≤ s j i < 5(n + 1) ms, i.e., the neuron with GID j spikes in the n-th time bin.
(3) There is some l such that 0 ms < s k l − s j i < 7.5 ms, i.e., the neuron with GID k spikes after the neuron with GID j, within a 7.5 ms interval.
We call the matrices A(n) transmission-response matrices, as it is reasonable to assume that the spiking of neuron k is influenced by the spiking of neuron j under conditions (1)-(3) above.
The goal of the Point vs. Circle test was to determine whether topological metrics, such as simplex counts, Betti numbers and Euler characteristic, could classify correctly two groups of stimuli of a similar nature and whether these metrics contain more information than the mean firing rate.In Figure 4C we provide plots of the time series of the average zeroth, first, and second Betti numbers, of the average numbers of 1-and 2-simplices, and of the average Euler characteristic for 20 trials of each stimulus.We applied a Gaussian Bayes classifier (SI, Supplementary Methods, SM2) to each metric in each time bin, to determine their success rate at classifying the various trials of the stimuli.To compare, we also classified the stimuli according to the mean firing rates.To allow for a fair comparison, we used three mean firing rates (between t to t + 5, t + 5 to t + 10, and t + 10 to t + 15 ms) for the classification at each time step t, since the transmission-response edges for time step t are based on information from up to t + 12.5 ms.
As illustrated by Figure 5 A, none of the metrics considered, topological or otherwise, succeeded very well at classifying the stimuli for times between 10 ms and 50 ms after the initial stimulation, which is not surprising given the strong similarity between the spatial propagation of activity of the two stimuli during this period (Figure 3).On the other hand, in the very first time bin, immediately after the initial stimulation, the 1-and 2-dimensional simplex counts and β 1 and β 2 all classify very well.In the second time bin the 2-dimensional simplex count and β 2 continue to classify very well, and the Euler characteristic classifies even better.Immediately after the second stimulation, from 50 ms to 55 ms after the initial stimulation, none of the metrics performs very well, but the 2-dimensional simplex count and β 2 still have the highest success rate.In the next time bin, from 55 ms to 60 ms after the initial stimulation, the 2-dimensional simplex count and β 2 again classify very well and are the only metrics to do so.In all of these cases, the topological metrics far outperform the metric based on firing rate.

Discussion
We have introduced topological analysis of directed graphs encoding structural or functional connectivity of digital reconstructions of neural microcircuits.We showed in particular that these directed graphs differed significantly from random graphs of both Erdős-Rényi-type and types taking into account biologically constrained, distance-dependent connection probabilities.The topological analysis revealed not only the existence of high-dimensional simplices representing the most extreme form of circuit "motifs" -all-to-all connectivity within a set of neuronsthat have so far been been reported for brain tissue, but also that there are a surprisingly huge number of these structures.We established moreover that topological methods effectively distinguish functional responses to distinct thalamic stimuli, introducing a new measure of the spatio-temporal activity responses generated by neural tissue.The results of our topological analysis of biologically realistic digital  reconstructions provide a convincing argument for considering topology as a useful mathematical tool for analyzing the structural and functional connectome of neural circuits.
Our results lead naturally to many new questions, most notably concerning the biological significance of the high-dimensional simplices and homology classes we have discovered in the digitally reconstructed neocortical microcircuits.We intend to explore these questions in future studies.In particular we hypothesize that the time series of different topological metrics could reveal an evolving spatio-temporal code that goes beyond either rate or timing information to one that incorporates the structural organization.Such metrics could yield a deeper understanding of how the structural organization constrains emergent functional states.Age-dependent changes in such digital reconstructions may help reveal even more complex topological structures with development, and changes introduced by synaptic plasticity may reveal structures associated with learning and memory.
We expect the topological approach to studying directed graphs that we implement here will also prove useful in applications of network science outside of neuroscience, in the study of networks exhibiting intricate directed connectivity patterns, such as gene and protein networks, VLSI circuits, and electrical grids.The obvious utility of the directed flag complex in these applications may also encourage theorists to establish results analogous to those established by Kahle concerning Betti numbers of undirected flag complexes of random graphs [6].

4.1.
Computation of flag complexes and their Betti numbers.We represent the directed flag complex of a directed graph by a reference-based data structure, using vectors to store the references to the simplices in the simplicial complex.The required storage space grows linearly with the number of vertices and with the number of edges.A publicly available C++ implementation of the code will be available on http://neurotop.gforge.inria.fr/.All homology computations carried out for this paper were made with F 2 coefficients, using the boundary matrix reduced by an algorithm from the PHAT [2] library.For further details, please see (SI, Supplementary Text, ST2).4.2.The Point vs. Circle experiment.The stimulated reconstructed microcircuit is innervated by 310 VPM fibers, whose horizontal centers of innervation are evenly distributed over the microcircuit (one fiber per mini-column).It is therefore possible to activate the microcircuit with topographically different stimuli by selecting only a subset of these 310 fibers.Here we used two different stimuli, a point and a circle, which were calibrated by adjusting the respective number of fibers to evoke an overall similar mean firing rate (i.e., close enough to prevent clearly distinguishing between the two stimuli simply by the mean population firing rate).The microcircuit was stimulated by synchronous spikes, similar to the whisker deflection experiment described by Markram et al. (2015).The point stimulus consisted of synchronous spikes in the 46 neighboring fibers of the center of the microcircuit, whereas the circle stimulus involved 56 fibers near the periphery of the microcircuit.The stimulation was repeated every 50 ms, but only the firing rates after the first two stimulations (at 0 and 50 ms) are overlapping.
We used a Gaussian naïve Bayes classifier [8], where we performed 500 classification trials, randomly choosing 15 trials of each stimulus to be part of the training data, and five trials of each stimulus to be part of the test data.We then obtained the mean ratio of successfully classified test data points using 500 different training and test sets.The classification of the firing rate used the firing rates of three consecutive time bins, to make it a fairer comparison, since the edges may contain firing rate information of more than two time bins, over a range of 12.5 ms.

4.3.
Computation of transmission-response matrices.Transmission-response matrices were calculated according to the specifications mentioned above, using a custom-written program in the Python programming language.It combined the matrix of synaptic connections (structural matrix), constructed as part of the standard reconstruction process of the BBP, with the spiking output of a simulation run and user-defined values for time steps ∆t 1 and ∆t 2 (5 and 7.5 ms in our analyses).For further details, please see (SI, Supplementary Methods ,SM1).

SI Appendix
To

SM1. Optimization of the parameters for the transmission-response matrices
The transmission-response matrices that allow us to analyze activity in an experiment (cf. the section on Functional Topology in the main body of the article) form a sequence depending on two parameters, ∆t 1 and ∆t 2 .The number of matrices in the sequence is the duration of the experiment divided by ∆t 1 .In other words for a given experiment of duration T and fixed ∆t i , we obtain a sequence of matrices , where N is the integer value of T /∆t 1 .For fixed values of ∆t 1 and ∆t 2 , the corresponding sequence {A(n)} N n=1 is obtained as follows.The spiking output of the simulation is first converted into lists of spike times, one for each neuron.Standard histogram methods, binning by ∆t 1 , are applied to each list to determine in which time steps a presynaptic neuron fired.For each time bin in which a particular neuron fired, the exact timing of its first spike in that bin is then compared to the full list of spike times of each neuron it innervates, to ascertain which of them had spiked at most ∆t 2 ms after the presynaptic neuron.(Spiking of a pair of neurons within ∆t 2 ms is ignored if they are not structurally connected.)For all pre-postsynaptic pairs satisfying this constraint on spike timing, the corresponding entry in the transmission-response matrix for that time step is set to 1 and all others to 0. More precisely, the (j, k)-coefficient of the binary transmission-response matrix A(n) corresponding to the n-th time bin is 1 if and only if the following three conditions are satisfied, where s j i denotes the time of the i-th spike of neuron j.
(1) The (j, k)-coefficient of the structural matrix is 1, i.e., there is a structural connection from the neuron with GID j to the neuron with GID k, so that they form a pre-post synaptic pair.
(2) There is some i such that n∆t 1 ms ≤ s j i < (n + 1)∆t 1 ms, i.e., the neuron with GID j spikes in the n-th time bin.
(3) There is some l such that 0 ms < s k l − s j i < ∆t 2 ms, i.e., the neuron with GID k spikes after the neuron with GID j, within a ∆t 2 ms interval.Starting with firing data from spontaneous activity in the reconstructed microcircuit, we generated sequences of 20 transmission-response matrices for ∆t i ∈ {1, 2, 5, 10, 20, 50, 100} ms, thus creating 49 such sequences corresponding to every possible choice of the pair (∆t 1 , ∆t 2 ).We refer to each of these sequences as the true transmission-response matrices corresponding to the pair (∆t 1 , ∆t 2 ).
In the rest of this section, we describe the procedure for optimizing the choice of the time intervals ∆t 1 and ∆t 2 so that the associated true transmission-response matrices best reflect the actual successful transmission of signals between neurons in the microcircuit.SM1.1.Properties of the transmission-response matrix.The nonzero coefficients in a transmission-response matrix are a subset of those in the structural matrix.Due to the partly stochastic behavior of the in silico microcircuit, the subset will vary even for subsequent applications of the same stimulus.In fact, even an exact repetition of the same conditions will lead to different transmission-response matrices, if the random number generator is seeded differently.It follows that the generation of the transmission-response matrices for a given stimulus should be considered as a stochastic process.With the correct choice of the parameters ∆t i , the matrices should reflect how the microcircuit processes a stimulus and thus take into account parameters of neural processing, such as pre-post synaptic interaction.
To find parameters ∆t 1 and ∆t 2 that maximize the degree to which neural processing is captured by the transmission-response matrices, we first develop a stochastic model for synaptic firing that takes into account neural processing and that depends on ∆t 1 and ∆t 2 .For the purpose of this analysis, we assume that the true transmission-response matrices are compatible with this model.
Based upon our model for synaptic firing, we formulate a simplified model that ignores neural processing.For this simplified model and for any choice of parameters ∆t 1 and ∆t 2 , we explain how to obtain transmission-response matrices from actual firing data, by shuffling the firing data appropriately, then applying the algorithm for generating a transmission-response matrix of the previous section.Finally, for each choice of the parameters ∆t 1 and ∆t 2 , we compare the true transmissionresponse matrices for spontaneous activity in the reconstructed microcircuit to those obtained by the simplified generation process.The parameters that we work with in the main body of the paper are the ∆t 1 and ∆t 2 that maximize the difference (measured by the ratio of the numbers of ones in the matrices) between the actual transmission-response matrices and those resulting from the simplified model.SM1.2.Stochastic model with neural processing.Fix time intervals ∆t 1 and ∆t 2 .Let A = (a ij ) denote the structural matrix of a reconstructed microcircuit, and let A(n) = (a n ij ) denote the transmission-response matrix of the n-th time bin, based on firing data from a trial of simulated activity in the microcircuit, for the given intervals ∆t 1 and ∆t 2 .By Condition (1) above, if a n ij = 1 for any n, then a ij = 1.It is reasonable to consider A to be static, at least over the time periods considered here.
We want to compute the probability that a n ij = 1, given that a ij = 1, so we need to determine on which parameters and properties this probability depends.According to the definition of transmission-response matrices, a presynaptic and a postsynaptic spike are required for a n ij to be 1.To simplify the analysis somewhat, we assume that each neuron n i has a time-dependent, instantaneous firing rate F i (t) that determines spiking probability at time t, i.e., spiking can be described as an inhomogeneous Poisson process.Under this assumption, the expected number m i ∆t1 (t 0 ) of spikes of neuron n i in the interval [t 0 , t 0 + ∆t 1 ] can be computed as If K i ∆t1 (t 0 ) denotes the probability that neuron n i spikes at least once in the interval [t 0 , t 0 + ∆t 1 ], then t0) , where P(λ) is the Poisson probability mass function with parameter λ at 0. (Recall that if X is a random variable that counts the number of spikes of neuron n i in the interval [t 0 , t 0 + ∆t 1 ], then P m i ∆t1 (t 0 ) is the probability that X = 0.) If the change in F i (t) is slow compared to ∆t 1 , then m i ∆t1 (t) ≈ F i (t) • ∆t 1 .Moreover, 1 − P(λ) ≈ λ for small values of λ.For small enough ∆t 1 , the expected number m i ∆t1 (t 0 ) of spikes of neuron n i will certainly be small, and change in F i (t) will be slow in compared to ∆t 1 , so that we may assume that For the postsynaptic spike the situation is more complicated.As there is a causal relation between presynaptic and postsynaptic firing, mediated by synaptic transmission, we need to estimate the conditional probability of at least one postsynaptic spike, given that at least one presynaptic spike occured.Let n i and n j denote neurons such that a ij = 1.Let s 0 ∈ [t 0 , t 0 + ∆t 1 ] denote the time of the first presynaptic spike in this interval.Let X j ∆t2 (s 0 ) denote the random variable whose value is the number of times neuron n j spiked in the time window [s 0 , s 0 + ∆t 2 ].Let Y i ∆t1 (t 0 ) denote the random variable whose value is the number of times neuron n i spiked in the time interval [t 0 , t 0 + ∆t 1 ].We need to estimate the conditional probability P X j ∆t2 (s 0 ) > 0 | Y i ∆t1 (t 0 ) > 0 .The nature of this interaction is very intricate and depends on the identities of the presynaptic and postsynaptic neurons, the spiking history of the presynaptic neuron before s 0 , and all other synaptic input the postsynaptic neuron received.It can be described as governed by some function G ij modulating the spiking probability of the postsynaptic neuron n j .This function takes as parameters the expected number of spikes of neuron n j in the interval [s 0 , s 0 + ∆t 2 ], the time t 0 , and the "spiking history" of the presynaptic neuron n i until s 0 , which we write as a function s i * (t) evaluated at s 0 , giving rise to the expression (s0),t0,s i * (s0)) .
Summarizing the analysis above, the following formula provides a good estimate of the probability that a n ij = 1 if a ij = 1, for small enough ∆t 1 and ∆t 2 , where s 0 denotes the time of the first presynaptic spike in the interval [n∆t 1 , (n + 1)∆t 1 ] and t 0 = n∆t 1 . (1) ) .This conditional probability encodes not only the distinctive features of the structural connectivity (via a ij ) but also the potentially stimulus-dependant neuronspecific firing rates (via F i and F j ) and their co-variation.Most crucially, it captures the stimulus-dependent functional modulation of postsynaptic firing by a presynaptic spike as well.We assume that the true transmission-reponse matrices capture the actual transmission of spikes according to the model of synaptic firing described by this formula.SM.1.3.Null hypotheses: no neural processing.We introduce here a simplified model of synaptic spiking that is based upon formula [1] but that ignores pre-post synaptic interaction.We then explain how to obtain transmission-response matrices that correspond to this simplified model from firing data arising from simulated activity.
We begin by setting each G ij to be the projection onto the first component, ignoring the pre-post synaptic interaction.After this simplification, the approximation obtained in the previous section now reads , where s 0 denotes the time of the first presynaptic spike in the interval [n∆t 1 , (n + 1)∆t 1 ] and t 0 = n∆t 1 , as before.Since this drastic simplification neglects the central aspect of neural computation -pre-post synaptic interaction -it gives rise to control cases for each pair of parameters (∆t 1 , ∆t 2 ) and each choice of firing rate functions F i (t).Comparison of the true transmission-response matrices for each pair of parameters to the corresponding control matrices for the same pair and a specific choice of the functions F i (t) will allow us to determine values for ∆t 1 and ∆t 2 for which the true transmission-response matrix optimally reflects neural processing.
We assume moreover that the individual firing rates consist of a neuron-dependent frequency that is up-or down-regulated by a global time series, i.e., that F i (t) = f (i) • F (t), for some function F (t) and some constant f (i) for each neuron n i .Transmission-response matrices corresponding to this simplified model for fixed ∆t 1 and ∆t 2 , which we call simplified transmission-response matrices, can be generated by first shuffling all recorded spikes from simulated activity in the reconstructed microcircuit, while preserving both the number of spikes per neuron and per time bin, then applying the usual transmission-response matrix generation method.SM.1.4.Optimization of parameters.The difference between the true transmissionresponse matrices and the control case described above is a consequence of the prepost synaptic interaction.Comparison with the control case enables us therefore to measure how well that interaction is captured in the true transmission-response matrices.In particular, it is reasonable to optimize the parameters ∆t 1 and ∆t 2 so that the difference between the true transmission response matrices arising from actual simulation data and those arising in the control cases is maximized, as a maximal difference indicates that the effect of the pre-post synaptic interaction is captured optimally by the true transmission-response matrices.
The comparison between the true transmission-response matrices and the control cases was carried out by first producing 20 true transmission-response matrices and 20 simplified transmission-response matrices based on firing data obtained from spontaneous activity in the reconstructed microcircuit for every pair (∆t 1 , ∆t 2 ), where ∆t i ∈ {1, 2, 5, 10, 20, 50, 100} ms for i = 1, 2. The number of ones in each matrix was then computed and the average taken over each set of 20 matrices.Since no stimulus was applied to the microcircuit, the averages computed are meaningful, since the firing data should be fairly homogeneous across the time bins.
The average number of ones in the transmission-response matrix arising from simulated actitivity in the reconstructed microciruit, as a function of ∆t 1 and ∆t 2 , is illustrated in Figure S1. Figure S2 shows the ratio of the average number of ones in the true transmission-response matrices to the average number of ones in the simplified transmission-response matrices, for various values of ∆t 1 and ∆t 2 .
In all cases we find that the maximum lies between ∆t 2 = 5 ms and ∆t 2 = 10 ms, leading us to choose to work with ∆t 2 = 7.5 ms.For ∆t 1 we find a maximum at 50 ms, but we use ∆t 1 = 5 ms (for which the maximum ratio is only slightly lower than for ∆t 1 = 50 ms) instead to avoid more than one spike per neuron per bin.

SM2. Gaussian Bayes classifiers
Suppose there is a distribution ρ over R × {c 1 , c 2 , . . ., c k }, where {c 1 , c 2 , . . ., c k } is a set of class labels.We can project ρ onto each of the coordinates to construct a real-valued random variable X and a class-label-valued random variable Y .We wish to build a classifier C : R → {c 1 , c 2 , . . ., c k } which will, for any real number, choose the most likely class to which it might belong.That is, where P (A|B) is the probability of A conditional on B and argmax a∈A f (a) denotes the element a ∈ A such that f (a) is maximal.This element of A will in practice always be unique.
Bayes' theorem states that A Bayesian classifier picks the class with the highest conditional probability, which using Bayes' theorem is Usually ρ itself is unknown and must be infered from sample data.We then also assume some model distribution to estimate ρ from these samples.The Gaussian Bayes classifier is the Bayes' classifier under the assumption that the distribution of each separate class is Gaussian.
After calculating the means and variances of the sample data within each of the classes separately, we model their respective distributions by the Gaussians N (µ ci , σ 2 ci ).If p(A) denotes the probability density function of A, then A common situation, such as in our analysis, is a uniform prior.A uniform prior over {c 1 , c 2 , . . ., c k } means P (Y = c i ) = 1/k for all i.If we assume a uniform prior, then the factor P (Y =c) p(X=x) is common to all classes and thus does not affect which class achieves the maximum.Thus we get the formula

SM3. Randomization of connection matrices and other control cases
We created four types of random matrices of sizes and connection probabilities similar to the connectivity matrices of the BBP reconstruction.
SM3.1.Generation of Erdős-Rényi random matrices.For this basic control we first computed the overall connection probability in the reconstruction and found it to be 0.8%.We then generated random, binary square matrices of size 3.1 × 10 4 , where 1's were placed at random off-diagonal in the matrix with probability 0.8%.SM3.2.Randomization preserving the distance-dependent connectivity between layers.Input for this randomization method were the structural matrix and the matrix of pairwise soma distances, both generated as part of the standard BBP reconstruction process.The rows and columns of both matrices were first grouped into N = 6 groups according to the layer of the neuron they correspond to.This effectively partitioned both matrices into N * N = 36 submatrices each.For each pair of submatrices, the soma distances were grouped into bins of size 25µm.Next, in the submatrix corresponding to each distance bin, we first replaced all 1's by 0's and then replaced randomly chosen 0's by 1's, so that the total number of 1's was preserved.Creation of autapses, i.e., a connection from a neuron to itself, was avoided by creating a separate bin for distances of 0µm.
The result was a connection matrix with the same number of connections between each pair of layers and the same distance-dependent connection probability between pairs of layers, to within 25µm, as the original matrix.SM3.3.Randomization preserving the distance-dependent connectivity between m-types.This randomization method was identical to the preceding randomization, preserving connectivity between layers, except that the neurons were partitioned initially into N = 55 groups of morphological types instead of only six layers.

SM3.4. Generation of connection matrices according to Peters' Rule.
For this control case, we started with a connection matrix that placed a connection not just where a synaptic connection was found in the reconstructed microcircuit, but between each pair of neurons whose arbors came within close proximity (closer than 3µm).The resulting matrix had approximately 16 times more connections than the structural matrix.These connections were then pruned randomly with a uniform probability until the same number of connections as in the structural matrix was attained.ST1.2.Simplicial complexes.An abstract oriented simplicial complex is a collection S of finite, ordered sets with the property that if σ ∈ S, then every subset τ of σ is also a member of S. A subcomplex of an abstract oriented simplicial complex is a sub-collection S ⊆ S that is itself an abstract oriented simplicial complex.Henceforth, we simplify terminology and usually refer to abstract oriented simplicial complexes merely as simplicial complexes.
The elements of a simplicial complex S are called its simplices.A simplicial complex is said to be finite if it has only finitely many simplices.If σ ∈ S, we define the dimension of σ, denoted dim(σ), to be |σ| − 1, the cardinality of the set σ minus one.If σ is a simplex of dimension n, then we refer to σ as an n-simplex of S. The set of all n-simplices of S is denoted S n .A simplex τ is said to be a face of σ if τ is a subset of σ of a strictly smaller cardinality.A front face of an n-simplex σ = (v 0 , ..., v n ) is a face τ = (v 0 , ..., v m ) for some m < n.Similarly, a back face of σ is a face τ = (v i , . . ., v n ) for some 0 < i < n.If σ = (v 0 , . . ., v n ) ∈ S n , then the i th face of σ is the (n − 1)-simplex σ i obtained from σ by removing the vertex v i .
A simplicial complex gives rise to a topological space by means of the construction known as geometric realization.In brief, one associates a point (a standard geometric 0-simplex) with each 0-simplex, a line segment (a standard geometric 1simplex) with each 1-simplex, a filled-in triangle (a standard geometric 2-simplex) with each 2-simplex, etc., glued together along common faces.The intersection of two simplices in S, neither of which is a face of the other, is a proper subset, and hence a face, of both of them.In the geometric realization this means that the geometric simplices that realize the abstract simplices intersect on common faces, and hence give rise to a well-defined geometric object.A geometric n-simplex is nothing but a (n + 1)-clique, canonically realized as a geometric object.An n-simplex is said to be oriented if there is a linear ordering on its vertices.In this case the corresponding (n + 1)-clique is said to be a directed (n + 1)-clique.
If S is a simplicial complex, then the union S (n) = S n ∪ • • • ∪ S 0 , which is called the n-skeleton of S, is a subcomplex of S. We say that S is n-dimensional if S = S (n) , and n is minimal with this property.If S is n-dimensional, and k ≤ n, then the collection S k ∪ . . .∪ S n is not a subcomplex of S because it is not closed under taking subsets.However if one adds to that collection all the faces of all simplices in S k ∪ . . .∪ S n , one obtains a subcomplex of S called the k-coskeleton of S, which we will denote by S (k) .The computational usefulness of coskeleta will become clear when we discuss homology computation (ST1.3).
Directed graphs give rise to abstract oriented simplicial complexes in a natural way.Let G = (V, E, τ ) be a directed graph.The directed flag complex associated to G is the abstract simplicial complex S = S(G), with S 0 = V and whose nsimplices S n for n ≥ 1 are (n + 1)-tuples (v 0 , . . ., v n ), of vertices such that for each 0 ≤ i < j ≤ n, there is an edge in G from v i to v j .Notice that because of the assumptions on τ , an n-simplex in S is characterised by the (ordered) sequence (v 0 , . . ., v n ), but not by the underlying set of vertices.For instance (v 1 , v 2 , v 3 ) and (v 2 , v 1 , v 3 ) are distinct 2-simplices with the same set of vertices.ST1.3.Homology, Betti numbers, and Euler characteristic.We now recall certain well known invariants of simplicial complexes arising in algebraic topology, which are preserved under a class of morphisms that is relevant in algebraic topology and that includes isomorphisms.These invariants serve to measure the "complexity" of simplicial complexes, from various topological perspectives, leading us to refer to them as metrics.
Homology is an important algebraic invariant of topological spaces.In this paper we use only mod-2 simplicial homology, computationally the simplest variant of homology, which is why we choose to work with it in applications, though other types of simplicial homology may provide deeper information.We do not give a complete account of homology here, but rather an elementary description of what it is and its basic properties.
Let F 2 denote the field of two elements, which we denote by 0 and 1.Let S be a finite simplicial complex.Define the chain complex C * (S, F 2 ) to be the sequence {C n = C n (S, F 2 )} n≥0 , such that C n is the F 2 -vector space whose basis elements are the n-simplices σ ∈ S n , for each n ≥ 0. In other words, the elements of C n are formal linear combinations of n-simplices in S with coefficients in F 2 .For each n ≥ 0, there is a linear transformation called a differential for every n-simplex σ, where σ i is the i-th face of σ, as defined above.Having defined ∂ n on the basis, one extends the definition linearly to the entire vector space C n .
The n-th Betti number β n (S) of a simplicial complex S is the F 2 -vector space dimension of its n-th mod 2 homology group, which is defined by Computing the Betti numbers is conceptually very easy.Let |S n | denote the number of n-simplices in the simplicial complex S. If one encodes the differential ∂ n as a |S n | × |S n+1 | -matrix D n with coefficients in F 2 , then one can easily compute its nullity, null(∂ n ), and its rank , rk(∂ n ), which are the F 2 -dimensions of the null-space and the column space of D n , respectively.The Betti numbers of S are then a sequence of natural numbers defined by The n-the Betti number β n counts the number of "n-dimensional holes" in the geometric realization of S. When S = S(G) is the directed flag complex of a directed graph G, both the simplices of S and these "n-dimensional holes" can be regarded as particularly important "metamotifs" [13] in the graph G.
It is easy to show that the n-th Betti number of a simplicial complex S is equal to that of its (n − 1)-st coskeleton S (n−1) , i.e., β n (S) = β n (S (n−1) ), for all n.This observation turns out to be computationally very useful, since there is no need to store the simplices of dimension less than n − 1 that are not faces of higher dimensional simplices in order to compute β n (S).In this paper it was exactly this trick that allowed us to compute the top dimensional homology of the 42 Ncomplexes we worked with.
Homology actually encodes far more information than what is intimated here, which can potentially be used for analyzing networks, but for the purposes of this article the description above will suffice.
If S is a simplicial complex and |S n | denotes the cardinality of the set of nsimplices in S, then the Euler characteristic of S is defined to be There is a well known, close relationship between Euler characterstic and Betti numbers [5], which is expressed as follows.If {β n } n≥0 is the sequence of Betti numbers for S, then χ(S) = n≥0 (−1) n β n (S).
See Figure 2A for a specific example.
ST1.4.Hasse Diagrams.A Hasse diagram, otherwise known as a directed acyclic graph, is a directed graph H = (V, E, τ ) with no oriented cycles.Hasse diagrams can be used to encode various combinatorial, geometric, and topological structures, such as posets and cubical complexes.Below we explain in detail how Hasse diagrams encode simplicial complexes.We include this discussion here because our computational algorithm (Algorithm 1) is based on this idea.
A Hasse diagram H is said to be stratified if for each v ∈ V , every path from v to any sink has the same length.Thus in a stratified Hasse diagram the vertices are naturally partitioned into disjoint strata, where every directed path from a vertex in the k-th stratum V k to any sink is of length k.In particular, the 0-th stratrum V 0 is the set of sinks of H.Moreover, for all e ∈ E, there exists k > 0 such that τ 1 (e) ∈ V k and τ 2 (e) ∈ V k−1 .Note that if H and H are isomorphic Hasse diagrams, and H is stratified, then so is H .
An orientation ς on a Hasse diagram H consists of a linear ordering < ς,v of the set E v of edges with source v, for every vertex v of H.If H = (V, E, τ ) and H = (V , E , τ ) are Hasse digrams equipped with orientations ς and ς , respectively, then a morphism of oriented Hasse diagrams from (H, ς) to (H , ς ) is a morphism of directed graphs (α, β) : H → H such that for every v ∈ V , the restriction of β to a set map E v → E α(v) preserves the orientation, i.e, if e < ς,v e for some e, e ∈ E v , then β(e) < ς ,α(v) β(e ).A morphism (α, β) of oriented Hasse diagrams is an isomorphism if α and β are bijections.A stratified Hasse diagram equipped with an orientation is called admissible.
Vertices in the k-th stratum of a stratified Hasse diagram H are said to be of level k.If k < n, and v, u are vertices of levels k and n respectively, then we say that v is a face of u if there is a path in H from u to v. If H is also oriented and therefore admissible, and there is a path (e 1 , ..., e n−k ) from u to v such that e i = min E τ1(ei) for all 1 ≤ i ≤ n − k, we say that v is a front face of u.Similarly, v is a back face of u if there is a path (e 1 , ..., e n−k ) from u to v such that e i = max E τ1(ei) for all 1 ≤ i ≤ n − k.We let Face(u) denote the set of all faces of u and Face(v) k the set of those that are of level k, while Front(u) and Back(u) denote its sets of front and back faces, respectively.See Figure S4 for an illustration of the concepts introduced above.
Example 1.If G = (V, E, τ ) is a directed graph, then G can be equivalently represented by an admissible Hasse diagram with level 0 vertices V , level 1 vertices E, and directed edges from each e ∈ E to its source and target.The ordering on the edges in the Hasse diagram is determined by the orientation of each edge e in G.
Every simplicial complex S gives rise to an admissible Hasse diagram H S as follows.The level d vertices of H S are the d-simplices of S.There is a directed edge from each d-simplex to each of its (d − 1)-faces.The stratification on H S is thus given by dimension, and the orientation is given by the natural ordering of the faces of a simplex from front to back.See Figure S5.
The Euler characteristic of a stratified Hasse diagram H = (V, E, τ ) is defined to be the integer It is easy to see that isomorphic stratified Hasse diagrams have the same Euler characteristic.It is also straight forward to show that if H is a stratified Hasse diagram associated to a simplicial complex S, then the Euler characteristic of H coincides with that of S.

ST2. Data structures and algorithms
In this section we describe our basic data structures and provide a detailed overview of the algorithm that constructs the directed flag complex associated to a directed graph.We also indicate briefly how our homology computations were performed.A publicly available C++ implementation of the code will be available on http://neurotop.gforge.inria.fr/.ST2.1.Data structures.We represent an admissible Hasse diagram H corresponding to the directed flag complex of a directed graph G = (V, E, τ ) by a reference-based data structure, using vectors to store the references to the vertices of the diagram.Each vertex v ∈ H stores the following information.
(1) Ver(v): A vector of the vertices of G determining the simplex of the flag complex to which v corresponds.( 2) Tar(v): A vector of references to the vertices that are targets of edges with source v.In particular, the required storage space grows linearly with the number of vertices and with the number of edges.For our complexity analysis below we assume that accessing any vertex, using Tar or Src, takes O(1) time.ST2.2.Creation of the directed flag complex associated to a directed graph.We describe our algorithm that creates a directed simplicial complex given a directed graph G.The output is a Hasse diagram H, stored as the data structure described above.The identifier Ver(v) of a vertex v in H, corresponding to a simplex σ in the directed flag complex, is the vector of vertices in G that represents σ.
For every level n ≥ 1 vertex v in H such that Ver(v) = [v 0 , . . ., v n ], the algorithm additionally records a vector U v of references to level 0 vertices u satisfying the following properties: 2-simplex associated to t.Therefore, the ordering of Ver(e) can be extended to ordering of Ver(t), as in Line 11.Thus all level 2 vertices corresponding to 2-simplices in the directed flag complex of G will be created by the algorithm.Also, since every simplex has a unique 1-dimensional front face, every 2-simplex will be created only once by this process.
Notice also that the if condition in Line 3 ensures that only triangles in G consisting of three edges oriented as (v 1 , v 2 ), (v 2 , v 3 ), and (v 1 , v 3 ) will give rise to level 2 vertices in H.It follows by induction that the analogous condition on orientations is then automatically satisfied for simplices of dimension greater than 2. To see this, fix n ≥ 2, and suppose that all simplices of dimension less than or equal to n have the desired property.Fix an n-simplex S = [v 0 , . . ., v n ] and u ∈ U S .By definition of the set U S , there is an edge from v i to u for every i ∈ {0, . . ., n}.Note that u ∈ U S for any S ∈ Tar(S).The previous iteration of the repeat-until loop (Line 6) created an oriented simplex from S together with u, of which u is the last vertex.Since the ordering of elements in S is a restriction of the ordering of elements in S, the ordering of a n + 1 dimensional simplex [v 0 , . . ., v n , u] restricted to any face yields the orientation of that face.It follows that Algorithm 1 does indeed construct a directed flag complex.
We now discuss the termination of Algorithm 1.If a level n vertex v is a face of a level (n + 1) vertex w, then the last vertex u in Ver(w) is not present in Ver(v), but is listed in U v .From Lines 12 and 21 of the algorithm it is clear that U w ⊂ U v and moreover that u ∈ U v .The cardinalities of the vectors U (−) are therefore decreasing for the newly created vertices.More precisely, for a vertex t and its faces s i , there exist i such that |U t | |U si |.Level n + 1 vertices are created only if there exist a level n vertex t such that U t = ∅.Since the cardinality of the U (−) decreases with each iteration of the repeat loop, the algorithm will terminate.
We remark finally that the size of the directed flag complex corresponding to a given directed graph G may be exponential in the size of G.In that case, the process of creation of a complex is usually stopped at some fixed dimension n.The time complexity of Algorithm 1 is proportional to the size of the output complex H, multiplied by maximal level of a vertex in H (due to the target-source search performed in Line 15) of the algorithm.ST2.4.Homology and Betti numbers.All homology computations carried out for this paper were made with F 2 coefficients, using the boundary matrix reduced by an algorithm from the PHAT [2] library.Comparing randomized and non-randomized transmission-response matrices: average number of ones in a true transmission-response (t-r) matrix divided by the average number of ones obtained when the recorded spikes were randomized before calculating the t-r matrix.Matrices were calculated from simulated spontaneous, ongoing activity with different values for ∆t 1 (in different colors) and ∆t 2 (along the x-axis).For each pair (∆t 1 , ∆t 2 ), matrices for 20 time steps were calculated, and the mean ratio is shown.Spikes were randomized by shuffling the identities of the firing neurons, thus conserving the number of spikes in any given time step and the total number of spikes fired by each neuron.

Figure 1 .
Figure 1.(A) A sparse visualization of the microcircuit (soma and dendrites only).Morphological types are color-coded, with m-types in the same layer having similar colors.(B) Examples of simplices in dimensions 0 through 3. (C) An example of a directed graph and its associated flag complex, in which there is one nsimplex for every directed (n + 1)-clique in the graph.(D) A graph depicting the average number of simplices in each dimension for the flag complexes associated to the reconstructed microcircuit (Ncomplexes) and the four types of random graphs considered, each with the same number of vertices as the reconstructed microcircuit, where shading indicates standard deviation, which was very small for all except the N-complexes.

Figure 2 .
Figure 2. (A) An oriented simplicial complex consisting of eight 2-simplices glued together along their 1-dimensional faces, together with a table of its Betti numbers and numbers of simplices in dimensions 0,1, and 2 and a computation illustrating that the Euler characteristic can be computed as the alternating sum of the Betti numbers or the simplex counts.(B) Graph depicting the average Euler characteristic of the reconstructed microcircuit (N-complexes) and of each of the types of random graph considered, where the whisker indicates standard deviation, which was very small, except for N-complexes and P-complexes.(C) Boxand-whisker plots depicting the Euler characteristics of 35 reconstructed microcircuits, seven for each individual rat.(D) Box-andwhisker plots depicting the 5th Betti number of 35 reconstructed microcircuits, seven for each individual rat.

Figure 3 .
Figure 3. (A) Average firing rate (top-down projection) in the stimulated microcircuit, plotted during the first 35 ms after the first stimulation at t=0 ms in the Point vs. Circle experiment.(B) Raster plots of the same 500 neurons randomly picked from layer 4, for two trials of the circle stimulus.(C) Population PSTH of all neurons in the microcircuit for three trials of the Circle stimulus.(D) Mean firing rate of the Circle and Point stimuli, between t and t + 5 ms, where light shading indicates the standard deviation and dark shading the error of the mean.

2 Figure 4 .
Figure 4. (A) Schematic representation of the transmissionresponse paradigm: there will be an edge from j to k in the graph associated to particular time bin if and only if there is a physical connection from neuron j to neuron k, neuron j fires in the time bin, and neuron k fires at most 7.5 ms after the firing of neuron j.Here, shading indicates the error of the mean.(B)Schematic representation of those firing patterns involving a presynaptic and a postsynaptic neuron that lead to an edge in the transmissionresponse graph, with a red block indicating successful transmission and a white block indicating lack of transmission.(C) Time series plots of the average value of the metrics 1D (number of 1-simplices), 2D (number of 2-simplices), β 0 (the zeroth Betti number, i.e., the number of connected components), β 1 (the first Betti number), β 2 (the second Betti number), and EC (the Euler characteristic) for the Circle and Point stimuli.Here, shading indicates the error of the mean.

Figure 5 .
Figure 5. (A) Times series plot for the first 80 ms of the 40 trials of the percentage of correct classifications performed by a Gaussian Bayes classifier based on each of the metrics FR (sequences of mean firing rates over three consecutive time bins), 1D (number of 1-simplices), 2D (number of 2-simplices), β 0 (the zeroth Betti number, i.e., the number of connected components), β 1 (the first Betti number), β 2 (the second Betti number), and EC (the Euler characteristic).(B) Graphs depicting the percentage of correct classifications performed by a Gaussian Bayes classifier based on each of the metrics in four particularly important time bins: from 0 to 5 ms (immediately after the initial stimulation), from 5 to 10 ms, from 50 to 55 ms (the time bin immediately after the second stimulation), and from 55 to 60 ms.

4. 4 .
Gaussian Bayes classifiers.The Gaussian Bayes classifier minimises the probability of misclassification under the assumption that the distributions are Gaussian.We randomly split the data into training and testing sets.Using the training set we model the distributions of the dot and circle classes by Gaussians N (μ dot , σ2 dot ) and N (μ circle , σ2 circle ) respectively.Assuming a uniform prior and Gaussian distributions, Bayes' theorem provides a classifier Class(x) = argmax

( 3 )
Src(v): A vector of references to the vertices that are sources of edges with target v.The admissible Hasse diagram H is thus represented by an ordered set of d vectors, where d is the maximal level in H, and where the i-th vector contains the references to all level i vertices.Let S int denote the size of integer data types, and for a given graph G = (V, E, τ ), let |V | and |E| denote the cardinalities of the corresponding sets.Each edge of the Hasse diagram is stored in two vertices of the diagram.If each reference requires S int storage, then we require O(|E| • S int ) space to store all references.In addition, each vertex stores the vector of vertices in V of the simplex in the flag complex of G to which it corresponds, which requires an additional O(S int •d) of space per vertex.The total size of a Hasse diagram is thus bounded by O((S int • d) • |V | + |E| • S int ).

Figure S1 .
Figure S1.Average number of ones in the true transmissionresponse matrices for different pairs of parameters (∆t 1 , ∆t 2 ) in a simulation of spontaneous, in-vivo-like activity (Ca 1.2)

Figure S3 .
Figure S3.(A-C) Examples of directed graphs.Graphs (A) and (B) are isomorphic, where the isomorphism is given by the map sending vertex a to 1, b to 2, c to 3, and d to 4. Graphs (A) and (B) are not isomorphic to graph (C).Vertex b in graph (A) is a sink, vertex a in the same graph is a source.Graph (C) has no sources or sinks, which explains the lack of isomorphism to graphs (A) and (B).

Figure
Figure S4.(A) A Hasse diagram that is not stratified, due to the edge from the vertex 1 to 5. (B) A stratified Hasse diagram, where vertices 5, 6, and 7 are the vertices of level 0, vertices 2, 3, and 4 are of level 1, and vertex 1 is of level 2. This is also an admissible Hasse diagram, where the outgoing edges are ordered from left to right.Vertex 2 is a front face of vertex 1, while vertex 3 is neither a front nor a back face of a vertex 1, and vertex 4 is back face of a vertex 1.

Figure S5 .
FigureS5.Top: The geometric realization of a simplicial complex consisting of seven 0-simplices (labeled 1,...,7), ten 1-simplices, and four 2-simplices.The orientation on the edges is denoted by arrows, i.e., the tail of an arrow is its source vertex, while the head of an arrow is its target.Bottom: The Hasse diagram corresponding to the simplicial complex above.Level k vertices correspond to k-simplices of the complex and are labeled by the ordered sets of vertices that constitute the corresponding simplex.Note that, e.g., vertex 23 is a back face of a vertex 123 and a front face of a vertex 234.