Networks of neurons connected by electrical synapses (gap junctions) can be represented as undirected graphs, where information can flow in both directions. Networks with chemical synapses, which impose a single direction of synaptic communication from the pre- to the postsynaptic neuron (Figures 1B2,B3), are more accurately represented as directed graphs. Sub-sampling networks of neurons experimentally has revealed small motifs of synaptic connectivity, but not large cliques of neurons (Song et al., 2005; Perin et al., 2011). Knowing the complete directed network of neurons, as we do in the case of the reconstructed microcircuit, enables us to detect all cliques, directed, and otherwise (Figure 1).

When the direction of connections is not taken into account, a great deal of information is lost. For example, in the undirected case, there is only one possible configuration for a clique of four fully connected neurons (Figure 2A1, left). However, in the directed case, there are 3⁶ = 729 possible configurations, as each of the six connections can be in one of three states (i → j, j ← i, or i ↔ j connection types; Figure 2A1 right).

A clique with reciprocal connections contains two or more cliques consisting only of uni-directional connections (Figure 2A2). When only uni-directional connections are considered, there are 2⁶ possible configurations of four fully connected neurons, which are of two types: those that contain cycles (40 configurations; Figure 2A3 left; Section 4.1.3) and those that do not (24 configurations; Figure 2A3 right). Directed cliques are exactly the acyclic cliques. The net directionality of information flow through any motif can be defined as the sum over all neurons of the squares of the differences between their in-degree and their out-degree (see Equation 2, Figure S1). Directed cliques have the highest net directionality among all cliques (Figure S1; Section 4.1.4). A clique that contains cycles always decomposes into directed cliques with the same number of neurons or fewer, at the very least any single connection between two neurons forms a 2-clique. A cyclical clique of three neurons therefore decomposes into three 2-cliques. Following the conventions in algebraic topology, we refer to directed cliques of n neurons as directed simplices of dimension n-1 or directed (n-1)-simplices (which reflects their natural geometric representation as (n-1)-dimensional polyhedra) (see Figure S2; Section 4.1.3). Correspondingly, their sub-cliques are called sub-simplices.

The presence of vast numbers of directed cliques across a range of dimensions in the neocortex, far more than in null models, demonstrates that connectivity between these neurons is highly organized into fundamental building blocks of increasing complexity. Since the structural topology of the neural network takes into account the direction of information flow, we hypothesized that emergent electrical activity of the microcircuitry mirrors its hierarchical structural organization. To test this hypothesis, we simulated the electrical activity of the microcircuit under in vivo-like conditions (Markram et al., 2015).

Stimuli, configured as nine different spatio-temporal input patterns (Figure 4A), were injected into the reconstructed microcircuit through virtual thalamo-cortical fibers in which spike trains were induced using patterns recorded in vivo (Bale et al., 2015; Figure S4; Section 4.7). These stimuli differed primarily in the degree of synchronous input received by the neurons. As expected, the neurons in the microcircuit responded to the inputs with various spiking patterns (Figures 4B1,B2,B4). We then calculated for each connected pair of neurons the correlation of their spiking activity (Figure 4C) and found a broad distribution of correlation coefficients, with only ~12% of connections where either the pre- or postsynaptic neuron failed to respond during all stimuli.

To avoid redundant sampling when testing the relationship between simplex dimension and activity, we restricted our analysis to maximal simplices, i.e., directed simplices that are not part of any higher-dimensional simplex (Section 4.1.2). A connection can be part of many higher-dimensional maximal simplices, unless it is itself a maximal 1-simplex. Despite the restriction to maximal simplices, we retained all information about the structure of the microcircuit because the complete structure is fully determined by its list of maximal simplices (Section 4.1.2). Correlations were calculated from histograms of the average spiking response (peri-stimulus time histogram, PSTH; bin size, 25 ms) to five seconds of thalamo-cortical input over 30 repetitions of a given input pattern (Figure 4B3). We then calculated the normalized cross-covariance of the histograms for all connections (Figure 4C; Section 4.8) and compared it to the number of maximal simplices associated with each connection in each dimension (see Figure 4D).

The neurons forming maximal 1-simplices displayed a significantly lower spiking correlation than the mean (Figure 4D), an indication of the fragility and lack of integration of the connection into the network. The mean correlation initially decreased with the number of maximal 2-simplices a connection belongs to, and then increased slightly. We observed that the greater the number of maximal 2-simplices a connection belongs to, the less likely it is to belong to higher-dimensional maximal simplices, with the minimum correlation occurring when the connection belongs to no simplices of dimension higher than 3. In higher dimensions, the correlation increased with the number of maximal simplices to which a connection belongs. While very high mean correlation can be attained for connections belonging to many maximal 3- or 4-simplices, the mean correlation of connections belonging to just one maximal 5- or 6-simplex was already considerably greater than the mean. These findings reveal a strong relationship between the structure of the network and its emergent activity and specifically that spike correlations depend on the level of participation of connections in high-dimensional simplices.

To determine the full extent to which the topological structure could organize activity of neurons, we examined spike correlations between pairs of neurons within individual simplices. These correlations increased with simplex dimension (Figure 4E, blue), again demonstrating that the degree of organization in the activity increases with structural organization. Spike correlation between pairs of neurons is normally an ambiguous measurement of connection strength because it is influenced by the local structure, specifically by indirect connections and/or shared inputs (Palm et al., 1988; Brody, 1999). However, since in our case the local structure is known and described in terms of directed simplices, we could infer how the local structural organization influences spike correlations. We compared the impact of indirect connections and of shared inputs on correlated activity by calculating the average correlation of pairs of neurons at different positions in a simplex when ordered from source to sink (Figure 4E, right panel). The number of indirect connections is highest for the pair consisting of the first (source) and last (sink) neurons (Figure 4E, purple), while the number of shared inputs is highest for the last and second-to-last neurons (Figure 4E, red). The first (source) and second neurons (Figure 4E, green) serve as a control because they have the smallest numbers of both indirect connections and shared inputs in the simplex.

We found that correlations were significantly higher for the last two neurons in the simplex, suggesting that shared input generates more of the pairwise correlation in spiking than indirect connections in directed simplices (p < 8 · 10⁻⁶, all dimensions except 1D). Moreover, the spiking correlation of the source and sink neurons was similar to the correlation of the first and second neurons (Figure 4E, purple and green), further suggesting that spike correlations tend to increase as shared input increases. These results hold for a range of histogram time bin sizes (Figure S5). The specific positions of neurons in local structures such as directed simplices therefore shape the emergence of correlated activity in response to stimuli.

Simplices are the mathematical building blocks of the microcircuitry. To gain insight into how its global structure shapes activity, it is necessary to consider how simplices are bound together. This can be achieved by analyzing the directed flag complex, which is the set of all directed simplices together with the set of all sub-simplices for each simplex (Figure S6, Section 4.1.2). The directed flag complex is a complete representation of the graph, including in particular the cycles neglected when examining directed simplices in isolation. The relationship between any two directed simplices depends on how they share sub-simplices. Just as any simplex can be realized as a polyhedron, a directed flag complex can be realized as a geometric object, built out of these polyhedra. If two simplices share a sub-simplex, the corresponding polyhedra are glued together along a common face (Figure 5A). The “shape” (or, more precisely, the topology) of this geometric object fully describes the global structure of the network.

To analyze directed flag complexes we computed two descriptors, the Euler characteristic and Betti numbers (Section 4.1.5). The Euler characteristic of a flag complex is given by the alternating sum of the number of simplices in each dimension, from zero through the highest dimension (Figure 5A). The Betti numbers together provide an indication of the number of cavities (or more precisely, homology classes) fully enclosed by directed simplices in the geometric object realizing the directed flag complex, where the dimension of a cavity is determined by the dimension of the enclosing simplices. The n-th Betti number, denoted β_n, indicates the number of n-dimensional cavities. For example, in Figure 5A, there is one 2-dimensional cavity (and therefore β₂ = 1) enclosed by the eight triangles; if an edge were added between any two non-connected nodes, then the geometric object realizing the corresponding flag complex would be filled in with solid tetrahedra, and the cavity would disappear. In the flag complexes of the reconstructions, it was not possible to compute more than the zeroth and top nonzero Betti numbers, as lower dimensions were computationally too expensive (Section 4.2.2). We could easily compute all Betti numbers for the C. elegans connectome, however, as it has many fewer nodes and edges (Figure S3).

The Betti number computations showed that there are cavities of dimension 5 (cavities completely enclosed by 5-simplices/6-neuron directed cliques) in all seven instances of each of the reconstructions (Bio1-Bio5, Figure 5B; Bio-M not shown). In contrast, the ER- and PR-control models have no cavities of dimension higher than 3, and the GB-model has no cavities of dimension higher than 4, demonstrating that there are not only non-random building blocks in the reconstruction, but also non-random relationships among them. We found as well that the information encoded in β₅ and the Euler characteristic together captures enough of the structure of the flag complex of a reconstruction to reveal subtle differences in their connectivity arising from the underlying biological data (Figure 5B, different colors).

Thus far we have shown that the structural network guides the emergence of correlated activity. To determine whether this correlated activity is sufficiently organized to bind neurons together to form active cliques and to bind cliques together to form active cavities out of the structural graph, we represented the spiking activity during a simulation as a time series of sub-graphs for which we computed the corresponding directed flag complexes. Each sub-graph in this series comprises the same nodes (neurons) as the reconstruction, but only a subset of the edges (synaptic connections), which are considered active, i.e., the presynaptic neuron spikes in a time bin of size Δt₁ and the postsynaptic neuron spikes within a time Δt₂ after the presynaptic spike (Figure 5C and Figure S7, Section 4.9). By considering subsequent, non-overlapping time bins of constant size Δt₁, we obtain a time series of transmission-response (TR) graphs reflecting correlated activity in the microcircuitry. We converted the time series of TR graphs in response to the different patterns of thalamo-cortical inputs (see Figure 4A) into time series of directed flag complexes. We found significantly more simplices in the TR graphs (Δt₁ = 5 ms, Δt₂ = 10 ms) than would be expected based on the number of edges alone (Figure 5D), indicating that correlated activity becomes preferentially concentrated in directed simplices.

The nine stimuli generated different spatio-temporal responses and different numbers of active edges (Figure 6A). The variation in Betti numbers and Euler characteristic over time indicates that neurons become bound into cliques and cavities by correlated activity (Figure 6A and Figure S8). When we plotted the number of cavities of dimension 1 (β₁) against the number of those of dimension 3 (β₃) (the highest dimension in which cavities consistently occur), the trajectory over the course of ~100 ms (Figure 6B) began ~50 ms after stimulus onset with the formation of a large number of 1-dimensional cavities, followed by the emergence of 2-dimensional (not shown) and 3-dimensional cavities. The decrease in β₁ began while β₃ was still increasing and continued until β₃ reached its peak, indicating that higher-dimensional relationships between directed simplices continued to be formed by correlated activity as the lower dimensional relationships subside.

Different stimuli led to Betti number trajectories of different amplitudes, where higher degrees of synchrony in the thalamic input produced higher amplitudes. The trajectories all followed a similar progression of cavity formation toward a peak level of functional organization followed by relatively rapid disintegration. The center of the projection of each trajectory onto the β₁-axis (its β₁-center) was approximately the same. Together, these characteristics of the trajectories reveal a stereotypical evolution of cliques and cavities in response to stimuli. These observations are consistent with experimentally recorded in vivo responses to sensory stimuli in terms of onset delay, response duration, and the presence of distinct phases of the response (Luczak et al., 2015).

To determine the neurons involved in this robust evolution of functional organization, we recorded the mean levels of spiking activity at different spatial locations within the microcircuit for one exemplary stimulus (Figure 6C). The activity started at depths that correspond to the locations of the thalamo-cortical input (Meyer et al., 2010; Markram et al., 2015), increasing in layer 4 and at the top of layer 6, before propagating downwards, reaching the top of layer 5 and the center of layer 6 as β₁ peaks, consistent with the finding that most directed simplices are in these layers. The transition from increasing β₁ to increasing β₃ coincided with the spread of the upper activity zone deeper into layer 5 and the top of layer 6, consistent with the presence of the highest dimensional directed simplices in these layers. The bottom activity zone also continued moving deeper, until it eventually subsided. As the top activity zone reached the bottom of layer 5, β₃ attained its peak. The zones of activity at the peaks of β₁ and β₃ are highly complementary: zones active at the peak of β₁ were generally inactive at the peak of β₃ and vice versa. The activity zone then remained in layer 5 until the cavities collapsed.

Finally, we applied the same stimulus to the reconstructions based on variations in the underlying biological data (see Figure 5B, Bio-1 to 5) and found similar Betti number trajectories, indicating that the general sequence of cavity formation toward peak functional organization followed by disintegration is preserved across individuals. On the other hand, we observed markedly different amplitudes, indicating that biological variability leads to variation in the number of high-dimensional cavities formed by correlated activity (Figure 6D). We also found that, unlike the case of different stimuli applied to the same microcircuit (Figure 6B), trajectories arising from different biological variations have different β₁-centers. In some cases, we observed reverberant trajectories that also followed a similar sequence of cavity formation, though smaller in amplitude. The general sequence of cavity formation and disintegration, however, appears to be stereotypic across stimuli and individuals.

Specializing basic concepts of algebraic topology, we have formulated precise definitions of cliques (simplices) and cavities (as counted by Betti numbers) associated to directed networks. What follows is a short introduction to directed graphs, simplicial complexes associated to directed graphs, and homology, as well as to the notion of directionality in directed graphs used in this study. We define, among others, the following terms and concepts.

Analyses of connectivity and simulations of electrical activity are based on a previously published model of neocortical microcircuitry and related methods (Markram et al., 2015). We analyzed microcircuits that were reconstructed with layer height and cell density data from five different animals (Bio-1-5), with seven microcircuits per animal forming a mesocircuit (35 microcircuits in total). In addition, we analyzed microcircuits that were reconstructed using average data (Bio-M, seven microcircuits). Simulations were run on one microcircuit each of Bio-1-5 and Bio-M. Each microcircuit contains ~31,000 neurons and ~8 million connections. Data about the microcircuit and the neuron models used in the simulations, as well as the connection matrices, are available on https://bbp.epfl.ch/nmc-portal/ (Ramaswamy et al., 2015).

Additional control models of connectivity were constructed by removing different biological constraints on connectivity. We created three types of random matrices of sizes and connection probabilities identical to the connectivity matrices of the reconstructed microcircuits.

We analyzed part of the C. elegans connectome (Varshney et al., 2011), consisting of 6,393 directed chemical synapses, obtained from www.wormatlas.org/neuronalwiring.html.

We performed simulations of neuronal electrical activity during stimulation with spatio-temporal patterns of thalamic input at the in vivo-like state (as in Markram et al., 2015), in the central microcircuit of Bio-M. Additionally, we repeated the same simulations in the central microcircuits of the Bio-1-5 reconstructions. We ran simulations using nine different organizations of thalamic input spike trains (see below).

We constructed postsynaptic time histograms (PSTHs) for each neuron for each stimulus, using the mean response to 30 trials of 5 s of thalamic stimulation (with bin size of 25 ms; for additional control, bin sizes of 10, 50, 100, 250, and 500 ms were also used). We then computed the normalized covariance matrix of the PSTHs of all neurons

where C_ij is the covariance of the PSTHs of neurons i and j. PSTHs of simulations with different thalamic stimuli were concatenated for each neuron to yield an average correlation coefficient for all stimuli. In total, correlations are based on the response of all neurons during 30 trials of nine stimuli for 5 s of activity (22.5 min).

The temporal sequence of transmission-response matrices associated to a simulation of neuronal activity of duration T is defined as

where the n-th matrix, A(n), is a binary matrix describing spiking activity in the time interval [n·Δt₁, (n + 1)·Δt₁ + Δt₂], and where N = T/Δt₁. The (j, k)-coefficient of A(n) corresponding to the n-th time bin is 1 if and only if the following three conditions are satisfied, where sij denotes the time of the i-th spike of neuron j.

In other words, a non-zero entry in a transmission-response matrix denotes a presynaptic spike, closely followed by a postsynaptic spike, maximizing the possibility of a causal relationship between the spikes. Based on firing data from spontaneous activity in the reconstructed microcircuit, we optimized Δt_i such that the resulting transmission-response matrices best reflect the actual sucessful transmission of signals between the neurons in the microcircuit (see Supplementary Methods). Unless noted otherwise, Δt₁ = 5 and Δt₂ = 10 ms were used throughout the study.

Analysis of the model and simulations was performed on a Linux computing-cluster using Python 2.7, including the numpy and scipy libraries (Jones et al., 2001), and custom Python scripts. We calculated p-values using Welch's t-test (scipy.stats), unless noted otherwise.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.