A hierarchy index for networks in the brain reveals a complex entangled organizational structure

Authors: Anand Pathak, Shakti N. Menon, Sitabhra Sinha

Abstract

Nervous systems are often schematically represented in terms of hierarchically arranged layers with stimuli in the “input” layer sequentially transformed through successive layers, eventually giving rise to response in the “output” layer. Empirical investigations of hierarchy in specific brain regions, e.g., the visual cortex, typically employ detailed anatomical information. However, a general method for identifying the underlying hierarchy from the connectome alone has so far been elusive. By proposing an optimized index that quantifies the hierarchy extant in a network, we reveal an architectural motif underlying the mesoscopic organization of nervous systems across different species. It involves both modular partitioning and hierarchical layered arrangement, suggesting that brains employ an optimal mix of parallel (modular) and sequential (hierarchic) information processing.

The brain exhibits complexity at multiple levels, from individual neurons interacting with neighboring cells, to the emergence of large-scale spatiotemporal patterns of activity in entire brain areas, ultimately giving rise to behavior and cognition (1). While the size and scale of nervous systems vary widely across species, it is striking that they nevertheless perform the key function of enabling the organism to respond appropriately to an ever-changing environment (2). The nervous system of the nematode Caenorhabditis elegans comprising ∼300 neurons lies at one end of this spectrum, while mammalian brains with tens of billions of neurons straddle the other extreme. Thus, the complexity of the brain does not simply arise from its size alone, e.g., the number of constituent neurons, but is also associated with the connection topology of the wiring between its constituent units (3). Understanding this structural organization, which underpins brain function, requires identifying general design principles that can provide a conceptual scaffolding for describing the connectome. Here, we focus on modularity and hierarchy, attributes that have often been associated with structural and/or functional features of brain organization. However, as these terms have been used with very different connotations depending on the context (e.g., see refs. 4, 5, 6, 7, 8 for instances of distinct ways in which hierarchy has been interpreted), their explanatory power has been limited. While the use of graph theoretic concepts has contributed to a rigorous and widely used framework for understanding modularity (9), even within the specific arena of network neuroscience there has been a multiplicity of approaches that seek to quantitatively characterize hierarchy.

Hierarchical organization has often been inferred, e.g., in the Macaque visual cortex (10–12), by observing how information flows across a sequential arrangement of layers, such that each successive layer integrates the signals obtained from the preceding layer and performs more complex information processing, thereby defining a bottom–up flow. In parallel, top–down feedback connections from higher processing levels to those at lower levels implement control mechanisms that allow adaptation and fine-tuning of responses (13, 14). Similar hierarchical organization has been reported in different species (15, 16) as well as other sensory modalities (17). Indeed, networks in general that are involved in complex information processing appear to be characterized by such an arrangement of reciprocal connections between nodes belonging to successive levels occurring in a sequence (18–20). This can, in principle, be related to the optimal use of computational resources by having successive layers receive information appropriately processed so as not to overwhelm the handling capacity of the constituent nodes. Distributing the processing task across the system, so that various steps are performed sequentially at successive layers, increases the robustness of the network against congestion-driven failure arising from bottlenecks that could result from increased computational load at a key node. In the brain, hierarchical organization is also hypothesized to be crucial for coordinating complex sequential behavior (21, 22), e.g., HVC neurons firing in a precise temporal order during the “song” of the zebra finch (23).

The lack of a universally accepted quantitative measure for hierarchy has meant that almost all earlier reports of hierarchical architecture in brain networks have relied on identifying the distinct layers through additional information about the attributes of the constituent nodes, such as their function and/or anatomical characteristics. However, in instances where such auxiliary knowledge is unavailable or incomplete, we require a procedure for unambiguously reconstructing the hierarchical sequence of layers from the connection topology alone. With rapid progress in methods to determine structural and functional brain networks in recent times, there is rising interest in devising novel approaches for analyzing the resulting abundance of connectome data to identify the inherent organizational features of such networks, in particular, hierarchy (24). As already implied above, hierarchical networks share the structural characteristic that their nodes are sequentially arranged into several layers with the densest connections occurring between successive layers. The presence of such an organization which intuitively has functional implications in terms of directing the flow along the network, especially for systems involved in processing information, can be used to fashion a quantitative metric of hierarchy. Thus, the hierarchical architecture of a network can in principle be disentangled by identifying the arrangement of all N nodes (say) of a network into L levels li (i=1,…,L) such that the number of connections between successive levels is maximized across all possible choices of a) the number of levels L, b) the sequential arrangement of li, and c) the level membership of each individual node.

In this paper, we propose a general framework based on this insight to identify the mesoscopic structure (in particular, the presence of a hierarchical organization) in a network exclusively from information about its connections and apply it to multiple connectomes. We define a metric, the hierarchy index H, such that maximizing it yields the optimal hierarchical decomposition of a network in the sense as described above. This is a nontrivial search problem, as even for a fixed value of L, the number of possible ways in which N nodes can be partitioned among the levels is given by the Stirling number of the second kind, yielding an astronomically high number of possibilities. For instance, a connectome comprising 200 brain regions can be arranged among 10 levels in more than 10193 different ways. We solve this combinatorial optimization problem by introducing a heuristic simulated annealing routine developed specifically for identifying the optimal partitioning of the network components into levels and their corresponding sequential arrangement. Benchmarking was carried out on synthetic networks with embedded hierarchical organization to establish that the hierarchy is correctly identified by the algorithm consistently across realizations.

Applying this method to various connectomes, ranging from the macroscale, consisting of tracts linking brain areas, to the microscale, comprising synapses and gap-junctions between neurons, we identify a robust mesoscale feature, viz., modular hierarchy. Indeed, our results suggest that the organization of brain networks is characterized by an interplay between the two prominent mesoscopic structural features, viz., modularity and hierarchy, such that neither can independently explain the trajectory of signals flowing through the nervous system, relayed from layer to layer and module to module. Note that, this concept is distinct from that of hierarchical modularity (25–29), which has been used in the literature to refer to a nested arrangement of modules. We show that the layered structure characterizing each module is not completely independent of that in other modules, suggesting a weak sequential order among the modules themselves rather than a dominant global hierarchy. Such an organization is consistent with the functional requirements of the nervous system which processes information in a segregated manner along specialized streams but eventually requires an overall integration. Taken in conjunction with recent experimental observations indicating that modules that are established initially subsequently get concatenated (30), our results suggest that this structural feature may well be developmentally programmed, pointing again to its potential functional relevance.

Results

The hierarchy index H that we define here (Materials and Methods) attains its highest possible value for a partitioning of the nodes into a sequence of levels that maximizes the density of links between adjacent levels. Thus, uncovering the hierarchical structure of a network is framed as a combinatorial optimization problem that we solve using simulated annealing (Fig. 1; see Materials and Methods for details). To establish the effectiveness of the proposed method in identifying the underlying hierarchical structure of a network, we first apply it to generated ensembles of benchmark random networks where such an organization is present by design and compare the inferred sequence and composition of levels with that known a priori. The process by which links are assigned between different nodes allows us to specify the extent of hierarchical organization, parameterized by the ratio h of the densities of connections between consecutive levels to that between all other levels (as well as, within each level). This allows us to smoothly vary the nature of the constructed networks from ones where hierarchy is completely absent (h=1), with the connections being uniformly distributed throughout the network, to those that are rigidly hierarchical (h=0), with nodes at any level allowed to connect only with those in levels immediately above or below them (see Materials and Methods for details). The benchmark networks can be decomposed into an optimal set of partitions by maximizing H which allows us to recover the mesoscopic topological organization embedded in the network to a remarkable accuracy, as measured by the normalized mutual information between the original and reconstructed hierarchical configurations. As expected, the performance of the algorithm declines as the hierarchical character of the network becomes less pronounced (for h>0.1), with a reduction in the similarity between the partitions identified by the algorithm and those inserted by construction (SI Appendix, Fig. S1).

Having validated the accuracy of the proposed hierarchical decomposition on benchmark networks, we apply this method to uncover any underlying hierarchical organization that may exist in several connectomes that vary in size and complexity as well as the scale of resolution of the network. These include the neuronal network corresponding to the somatic nervous system of the nematode C. elegans (31), an aggregated network of brain regions compiled from numerous tractographic studies of the Macaque brain (32–35), and structural brain networks of multiple human subject obtained through diffusion tensor imaging (36, 37) (see Materials and Methods for details). We describe below the results of the hierarchical decomposition carried out on these networks. We begin with the Macaque macroconnectome, whose modular organization and resulting functional consequences we have recently investigated in detail (38).

A typical hierarchical decomposition of the Macaque connectome obtained by applying our algorithm is shown in Fig. 2 A and B. The specific sequence displayed here comprises 16 levels into which the brain regions are arranged and is taken as the reference sequence (Dataset S1) against which we compare other decompositions generated by multiple realizations of the algorithm (Materials and Methods). As the diameter of the network (=8) is much smaller than the number of layers into which the connectome is seen to be partitioned in these decompositions, it may initially appear counterintuitive that the shortest path length connecting the two regions furthest in terms of the distance measured along the network is not comparable to the separation between the terminal layers in the hierarchical chain. Such a feature suggests a marked deviation from a strict hierarchy (characterized by links existing exclusively between neighboring layers) and has always proved challenging to any effort at describing the brain in terms of a serial arrangement of layers that successively process information (8). This apparent incongruity arises because of a profusion of “short-cuts” linking regions that lie in layers that are far apart along the sequential arrangement. This can be established by observing how the diameter of synthetic networks having a comparable number of nodes and layers as the connectome, decreases as the density of links connecting nonconsecutive layers increases (SI Appendix, Fig. S2). Indeed, it is the presence of these connections which obscures the underlying hierarchical arrangement of brain networks, a problem that has been overcome by the hierarchical decomposition method introduced here.

Almost all the decompositions exhibit a spatially contiguous arrangement in that the sequentially adjacent levels also appear to be spatially adjacent. We observe that the levels exhibit a cyclic progression from the anterior to posterior before eventually turning back. As seen from the sagittal section (Top Right panel of Fig. 2A), the sequence begins at the prefrontal cortex (nodes in layers 1 to 3) and then moves across the parietal lobe (layers 4 to 7) down to the subcortical regions (layers 8 to 10) before proceeding up again to the occipital lobe (layers 11 to 13). The subsequent levels then progress in the reverse direction (see the horizontal section in the Left panel of Fig. 2A) across the temporal lobe (layers 14 and 15) to finally terminate in the prefrontal cortex (layer 16). Thus, the terminal levels of the hierarchy are both located in the frontal lobe. Fig. 2B shows that most of the connections between brain regions tend to be concentrated between consecutive layers (whose nodes occur within the partitions indicated by the red bounding lines) in the hierarchical sequence, consistent with the intuitive notion of hierarchy that we outline earlier. As can be seen, the sizes of the layers, measured by the number of regions that belong to each of them, are highly variable, ranging from 3 (Layer 10) to 36 (Layer 5). However, as the brain regions themselves occupy very different spatial volumes, spanning several orders of magnitude, the size differences between the layers in terms of the number of regions may not easily translate to variation in their spatial scale.

As noted earlier, it has already been shown that the Macaque connectome has a prominent modular organization of the brain regions, defined by communities characterized by dense intraconnectivity that are spatially localized to a large extent (38). Fig. 2C shows how the two mesoscopic organizational features of the brain, viz., modularity and hierarchy, relate to each other. As can be seen, each module comprises brain regions that largely belong to sequentially adjacent hierarchical layers, such that we can categorize the network as one that is composed of modular hierarchies. In other words, the connections can be partitioned into several modules, each of which can be further decomposed into a series of hierarchical layers. We note that the hierarchy is defined not only in terms of the sequence of layers within each module, but the different modules themselves occur in the decomposition in a specific order.

As the modularity of the network can potentially interfere with the determination of the hierarchical sequence, given that both types of mesoscopic organization are based on the differential connection densities within and across partitions, we have carried out another series of benchmark tests of the decomposition algorithm aimed specifically at such networks whose nodes are arranged into modules, as well as, hierarchical layers. For this purpose, we have constructed an ensemble of synthetic random networks with a hierarchical structure as determined by the parameter h (defined above), and whose modular character is parameterized by the ratio r of the densities of connections between nodes belonging to the same community (ρo) and those belonging to different communities (ρi) (39). Benchmark networks are obtained for given pairs of values of h and r, comprising several modules whose nodes are in turn arranged into multiple sequentially arranged layers, i.e., the networks embody a modular hierarchical architecture (SI Appendix, Fig. S3). The networks are then decomposed by maximizing H and the partitioning thus obtained can be compared with the embedded structure by computing the mutual information between them. We observe that the algorithm uncovers the underlying hierarchical organization in the presence of modules, when both the hierarchical and modular characters of the network are prominent (i.e., for low h and r). We note that the embedded organization can be detected with an accuracy that is comparable to that obtained for the exclusively hierarchical synthetic network ensembles (described above).

As the process for partitioning of the network into hierarchical layers (and modules) is stochastic in nature, different realizations of the decomposition can result in distinct sets of network partitions, and the sequence in which they are arranged may also vary. A network with an inherent hierarchical organization should display broad consistency across the various decompositions, both in terms of the membership of the different layers as well as their sequential order (Materials and Methods). Fig. 2D shows that there is strong agreement between the different sequences of layers in the Macaque connectome obtained from multiple realizations, as indicated by the prominent diagonal (see also SI Appendix, Fig. S4 A–C). This implies that specific levels in the reference sequence (ordinate) occur with high relative frequency at corresponding positions in the sequences obtained from the other decompositions (abscissa). This is complemented by Fig. 2E, which shows that the relative frequency with which a particular brain region occurs in a specific layer across different realizations is strongly localized to a single partition. In other words, the layer memberships of the brain regions are highly correlated, i.e., if a pair of brain regions belong to a particular layer in a sequence (considered as the reference), then they will co-occur in a layer in other hierarchical decompositions of the network with very high probability. Thus, the Macaque connectome shows a very robust hierarchical organization embedded within the inherent modular structure of the network, with the composition of individual layers and the order in which they occur sequentially being largely invariant across realizations of the decomposition algorithm.

By identifying a network architecture that can reproduce observed features arising from the relation between modular and hierarchical characteristics of the connectome, we can get vital clues about its organizational principle. In particular, it will have to demonstrate how a robust, globally sequential ordering of the identified layers can be consistent with the embedding of these layers into prominent modules, which are by definition relatively independent of each other. With this aim in mind, we have considered two classes of synthetic networks possessing both modularity and hierarchy. One of these classes represents networks that are characterized by independent modular hierarchies (Fig. 3A). In these networks, there is a definite sequential ordering of the layers within each module, but not of the modules themselves. As the relatively sparse number of connections between the different modules may connect nodes at any level in a given module to those occurring at any level in another, in principle the modules can be placed in an arbitrary order without disrupting the hierarchy of the network as a whole. To contrast with this, we consider another ensemble of networks in which there is not only a strict ordering of layers within each module, but also across the modules (Fig. 3B). In other words, each module is tethered to a specific position in the sequence relative to the other modules, which arises because most of the connections between consecutive modules in the sequence occur between their respective terminal layers. We note that this implies a rigid sequential arrangement across the modules which would appear to partially contradict the fundamental attribute of relative independence that characterizes modular structure. By decomposing these two classes of networks using the algorithm presented here, we find that while both ensembles are characterized by layer memberships that are consistent across realizations, only networks with independent modular hierarchies exhibit a robust sequential order of the identified hierarchical layers, as is observed in the case of the empirical network. It suggests that the relation between modularity and hierarchy in the connectome is closer to that represented by networks where the hierarchical arrangement in each of the modules is relatively independent.

We next investigate the hierarchical organization of a human connectome, obtained from a representative individual subject (Materials and Methods). Fig. 4A shows a specific decomposition that is chosen to be the reference sequence (Dataset S2). The 188 brain regions across the two hemispheres that comprise the connectome are seen to be partitioned into 12 layers. The regions belonging to consecutive layers are also physically adjacent, as can be observed from their spatial locations in the horizontal section of the brain shown in Fig. 4A. As in the case of the Macaque, the number of hierarchical layers is larger than the diameter of the network (=4), which can be attributed to the many connections across nonconsecutive layers functioning as “short-cuts” (see SI Appendix, Fig. S5 which shows the dependence of the diameter of equivalent synthetic networks on the ratio h). Fig. 4B shows that these short-cuts are quite substantial in number. Indeed, they are more numerous in the human connectome compared to that of the Macaque, which can possibly be associated with the much higher overall connection density in the former (40). The network also exhibits modular organization characterized by the existence of four modules, two of which mostly comprise regions from the left hemisphere while the other two have a majority of their members in the right hemisphere. The relation between the compositions of the modules and the hierarchical layers is indicated by the alluvial diagram in Fig. 4C. It suggests that similar to the Macaque connectome, the network can be viewed as possessing modular hierarchies. Again as in the Macaque, we see that the network displays a robust sequential ordering of the layers (Fig. 4D; see also SI Appendix, Fig. S4 D–F), with the identities of the members of each layer being broadly consistent across different realizations of the hierarchical decomposition (Fig. 4E). We have carried out similar hierarchical decompositions of other human connectomes obtained from subjects of different ages, which are partitioned by our algorithm into a similar number of hierarchical layers ranging from 11 to 14 having robust sequential arrangement, as well as, layer membership (SI Appendix, Fig. S6).

As a final demonstration of the proposed hierarchical decomposition method, we consider the network of chemical synapses connecting 279 neurons which belong to the somatic nervous system of C. elegans (Materials and Methods). A typical partitioning of the network, chosen to be the reference sequence (Dataset S3), is shown in Fig. 5A, where the membership of the neurons among the 12 layers that are obtained for this realization are indicated using different colors. As in the case of the networks of brain regions analyzed above, the diameter of this neuronal network (=7) is seen to be lower than the number of layers identified. This can be imputed to short-cut connections spanning nonadjacent layers (SI Appendix Fig. S7). Considering the spatial positions of the neuronal cell bodies, we observe that the initial layers are concentrated around the nerve ring located at the head of the organism. Subsequent layers have neurons that are located in the tail, while the neurons of the ventral cord (laid out along the anterior–posterior axis of the worm body and consisting almost exclusively of motor neurons that coordinate locomotion) occupy the final layers in the sequence. Thus, we see a deviation from the spatial contiguity of consecutive hierarchical layers that marked the connectomes of the macaque and human. This could possibly be a consequence of many of the neurons having long processes that span almost the entire length of the organism, such that their connections are not just confined to the vicinity of the cell body (41, 42). Fig. 5B shows the neurons clustered in the various ganglia that are located in the head of the organism. In this magnified view, we observe that the initial layers broadly appear to be spatially ordered, with their composition being dominated by sensory and interneurons (indicated by the shape of the symbol representing each neuron, see key in Fig. 5A). This is substantiated in Fig. 5C, which indicates the number of neurons of different types that occur in each layer (Top panel), as well as their relative fraction in these layers (Bottom panel), averaged over multiple realizations of the hierarchical decomposition algorithm. The variation in the size of each layer (both in terms of the total number of neurons, as well as, of each functional type) appears to be relatively low across the different partitionings. The initial layers appear to have a larger fraction of sensory neurons, while the interneurons predominate the composition of layers that occur in the middle of the sequence. Motor neurons, on the other hand, constitute the bulk of the last few layers.

The adjacency matrix for C. elegans neuronal connectivity shown in Fig. 5D illustrates the dense interlevel connectivity between consecutive levels, which is characteristic of a strongly hierarchical organization. As in the cases of the networks of brain regions described above, this network also has a modular organization comprising three modules (the association between their neuronal composition and that of the layers is shown in SI Appendix, Fig. S8, Left panel). One of the modules comprises neurons that are mostly located in the ventral cord, while the neurons in the various ganglia are divided among the other two modules. Each module exhibits a distinct hierarchical arrangement of layers within them, suggested by the relatively little overlap in the modular memberships of neurons in each layer. Thus, it appears that the “modular hierarchy” organization principle that we observe at the scale of connections between brain regions could also be operating at the neuronal scale. Also, as in the other networks, the identities of nodes belonging to the different layers are consistent across realizations, indicating the robustness of the decomposition in terms of layer membership of the neurons (Fig. 5E). While the sequential ordering of the layers shows more variability (SI Appendix, Fig. S4 G–I) than that seen in the case of other networks, the reference sequence is largely conserved across the plurality of the obtained decompositions. The robustness of the sequence of hierarchical layers is even more pronounced when we augment the network with additional links (approximately one-third of the number of synapses) corresponding to electrical gap junctions between neurons (SI Appendix, Figs. S9 and S10). Indeed, the network obtained by incorporating these links also exhibits hierarchical organization having a similar number of layers, with the identities of the neurons belonging to each layer remaining consistent across decompositions. Moreover, the relation between the modules of this network and its hierarchical layers also supports the hypothesis that the networks comprise relatively independent modules, each with an embedded set of hierarchically ordered layers (SI Appendix, Fig. S8, Right panel).

Discussion

The mesoscopic organization of a network is expected to reflect its function (43). For instance, the necessity of performing multiple independent tasks in parallel, with relatively low requirement for coordination between them, may favor a modular architecture. A network with such a functional requirement can be partitioned into a number of subnetworks, each characterized by high intraconnection density facilitating recurrent communication between their constituent nodes, while having correspondingly fewer connections between nodes belonging to different subnetworks. On the other hand, a hierarchical network may be preferred if the function typically requires performing several steps in sequence (such that each step needs to be finished before initiating the next), possibly coordinating across many input streams. Such a connection topology would promote efficient serial processing, often in conjunction with feed-back and feed-forward connection across the levels. As we show here, the connection architecture of the brain manifests both of these fundamental organizing principles.

Indeed, our analysis of the connectomes suggests a structural feature at the mesoscopic level in these networks that we term modular hierarchies. These are characterized by the brain regions being segregated into distinct communities, while at the same time being arranged in a specific sequence of levels within their own community. The robustness of the modular partitioning, as well as, the hierarchical sequence, suggests that both of these features are fundamental attributes of the network organization. In fact, while there have been previous attempts to identify signatures of hierarchy in the brain, we venture that it is the simultaneous presence of a strong modular arrangement that has made such an attempt particularly challenging. The method proposed here is particularly suited for identifying the interplay of these two kinds of mesoscopic organization. We note that a similar architecture is known to be extant in the visual cortex, where information coming from different parts of the visual field are processed by different microcolumns (hence, representing a modular partitioning) with each microcolumn being composed of a sequential arrangement of neurons that process the information from the specific part of the visual field received upstream (44). Thus, it appears that a common organizing principle may be operating at both the micro- and macroscale of the connectome.

It is of interest that the various hierarchical modules themselves appear to be sequentially arranged but is consistent with a model postulating that the hierarchical organization in each of the modules is independent of the others (Fig. 3). We note that this parallels recent evidence in mice about the manner in which the visual network develops, with independent modules being initially established and subsequently concatenated (30). It suggests a process by which the modular hierarchies may be developmentally programmed, lending further support to the organizational principle suggested by the network structural analysis presented here. In addition, a plausible functional relevance of this architecture is hinted at by the robust sequential relation between the various modules (depicted in SI Appendix, Fig. S4). The observations are consistent with the possibility that intermodular connections may have a preference for connecting a node in a given layer a of a particular module with that in layer b of a different module, where these two layers belong to an overall hierarchical sequence. We note that if such a global arrangement was strictly enforced, nodes would have been observed to be organized into layers following a global hierarchical arrangement, which would have been independent of the specific modules to which the individual nodes belong. On the other hand, in the absence of any such preference, the hierarchical level to which a node belongs will have no effect outside of its module and intermodular links should be equally likely between any pair of levels. Our results, while not supporting a global hierarchy independent of the modules, do appear to suggest a certain preference for sequential arrangement across modules. Functionally, this may provide a basis for systematic integration of information and hence allow for distributive processing in a network that otherwise has a markedly modular organization and hence would have appeared to support a segregated (or specialized) mode of processing (45). Indeed, such a relation between modularity and hierarchy has been hypothesized in the specific case of the visual system, as a possible solution to the ill-posed problem of feature binding (46).

Our analysis of the modular hierarchy present in three connectomes, viz., the network of white matter tracts linking brain areas in the macaque and the human, and that of neurons in the nematode C. elegans, using the hierarchy detection algorithm we introduce here, reveals various facets of the intriguing interplay between hierarchy and modularity. In the first brain network we investigate, viz., that of a single hemisphere of the macaque, our results suggest the existence of two distinct streams along which signals could propagate parallel to the anteroposterior axis. One of these extends from the frontal lobe to the occipital via the parietal lobe, whereas the other extends from the occipital lobe to the frontal via the temporal lobe. As the occipital lobe comprises the primary visual cortex, belonging to one of the first layers in which sensory information is processed, we can distinguish the two pathways identified above tentatively as “downstream” and “upstream,” respectively. While the latter may correspond to sensory stimuli being successively analyzed in brain areas that perform higher-level processing, the former can plausibly be involved in sending feedback signals back to the initial layers. We note that an analogous model of bottom–up and top–down processing working in conjunction has been proposed in the context of vision (13, 14, 47–49). Indeed, the macaque visual system has been reported to comprise several stages of hierarchical processing, the number of layers estimated being comparable to that obtained using our algorithm (which is distinct from the procedures used by these earlier studies) (8, 50).

For the individual human connectomes, as in the case of the Macaque, we observe a clearly identifiable arrangement of brain regions into modular hierarchies, that is consistent across subjects in terms of the number of layers and their sequential ordering. Thus, across different individuals we observe that regions occurring in the terminal layers are located either in the frontal or in parietal lobes. However, the progression from the initial layers can either be from the left to the right hemispheres or the reverse, depending on the subject whose connectome is being analyzed, suggesting that the connectivity from the fronto-parietal regions to the left and right hemispheres is asymmetric. More generally, the structural organization of the human brain uncovered using the techniques employed here closely resembles the hierarchically layered architectures employed by artificial neural networks implementing deep learning (51). These in turn have been inspired by ideas that, in the brain, information is processed through several sequentially arranged stages of transformation and representation. Indeed, recent studies using MEG and fMRI suggest a strong correspondence between object representation in the various layers of designed artificial deep neural networks and the hierarchical topography of visual representations in the human brain (52). It has long been assumed that the hierarchical organization of the brain is reflected functionally in a sensory stimulus being represented at different levels of abstraction in the successive layers, e.g., in the context of the visual system, the various stages respond successively to edges, primitive shapes, etc., and eventually to complex forms (53). The mammalian brain is, from the context of such layered neural network architectures, believed to implement a deep architecture with many layers, where depth corresponds to the number of sequentially arranged stages of nonlinear operations, applied on the output of the immediately preceding layers (54). Intriguingly, our results show that connectomes exhibit aspects of “deep,” as well as, “shallow” architecture. Specifically, the various modules that we identify typically comprise 3 to 4 layers, while the network viewed in totality can be seen to possess a much higher number of layers, e.g., 15 to 16 in the macaque, as the modules are themselves arranged in a sequential order. Thus, it appears that there may not be a strict dichotomy between shallow and deep architectures, but rather they are associated with the scale at which one analyzes the network organization of the brain.

The connectomes of the macaque and human discussed above differ in a fundamental manner, in that the former is directed while the latter is undirected. The method proposed here is nevertheless able to detect broadly similar hierarchical organizations in both, which points to the robustness of the technique to the occurrence of directed links. This is important in the context of our analysis of the neuronal network of C. elegans, where we focus on the directed network comprising chemical synapses. However, the neurons are also connected by electrical gap junctions, that in principle allow bidirectional communication, and hence can be viewed as constituents of an undirected network. Considering these networks of synapses and gap junctions together brings up additional challenges as the co-occurrence of directed and undirected links in the same network can obscure the hierarchical nature of the connections. Furthermore, the size of the organism is small relative to the typical length scales of the components of its nervous system, with the longest neuronal processes spanning almost the entire body. This implies that there may not be a strict correspondence between the sequential arrangement of the hierarchical layers and the spatial proximity of neurons occurring in neighboring levels, unlike in the macroscale networks for the macaque and human.

The hierarchical architecture reconstructed from the neuronal network (around three-quarters of whose links comprise synapses) agrees with our intuitive notion of how signals are processed through the nervous system following the stimulation of specific sensory organs, encountering in turn, sensory, inter-, and motor neurons, the latter serving as actuators for possible muscle activity. The fact that our algorithm is able to do this despite the specific challenges of analyzing the nematode nervous system, not only highlights its effectiveness in determining hierarchical organization across scales (from macroconnectome of brain areas to microconnectome of neurons) but also underscores the ubiquity of the hierarchical architectural plan of the nervous system. To conclude, the algorithm that we present here provides a comprehensive method for uncovering the hierarchical organization of networks appearing in very different species, effective across scales, nature of links (viz., directed or undirected), and the existence of other structural features such as modules. Indeed, our results provide a perspective on debates concerning the extent that processing in the brain is sequential (as in a hierarchically layered system) as opposed to being compartmentalized (as in modular systems), by suggesting an architecture, viz., modular hierarchies, combining aspects of both these mesoscopic organization principles.

Materials and Methods

Data.

We have considered empirical connectomes comprising unweighted links corresponding to anatomical tracts between brain regions in the macaque and human, and synapses and gap junctions between neurons in the nematode C. elegans. In cases where the original data have weights associated with each link, we consider only the adjacency matrix of the network.

Macaque.

Connectivity.

We have used a reconstructed macaque structural connectome comprising 266 cortical and subcortical brain regions, with 2,602 directed links between them, as described in ref. 38. It is a revised version of an earlier database (35), compiled from more that 400 separate tract tracing studies cataloged in CoCoMac, a comprehensive neuroinformatics electronic archive (32–34).

Spatial information.

The stereotaxic coordinates and the volume of each brain region in the connectome have been obtained from several sources, including the website https://scalablebrainatlas.incf.org/macaque/PHT00 associated with the Paxinos Rhesus Monkey Atlas (55), as well as, manual curation from the relevant research literature (for details, see ref. 38).

Human.

Human brain structural connectomes were chosen from those of subjects in the Nathan Kline Institute (NKI)/Rockland Sample (36) repository of diffusion tensor imaging (DTI) data, made publicly available by the UCLA multimodal connectivity database at http://umcd.humanconnectomeproject.org/ (37) as undirected connectivity matrices. The three-dimensional coordinates locating each brain region in a standardized space have been obtained from the above-mentioned database.

C. elegans.

Connectivity.

Information about the directed connections (corresponding to synapses), as well as undirected gap-junctions, between the 279 connected neurons of the C. elegans somatic nervous system has been obtained from the dataset published in ref. 31.

Functional type.

Information about the functional type of each neuron, i.e., whether it is a sensory, inter-, or motor neuron, has been obtained from the database provided in ref. 56.

Spatial information.

Coordinates of each neuronal cell body projected on a two-dimensional plane defined by the anterior–posterior axis and the dorsal–ventral axis, were obtained from the database associated with ref. 57, accessible online from https://www.dynamic-connectome.org/.

Hierarchy Index.

In analogy with the intuitive notion of hierarchy as a sequential ordering of items (58), we consider a network to be hierarchically organized if its nodes can be partitioned among multiple levels that are arranged in a specific sequence, with a discernible preference for nodes in neighboring levels to be connected to each other. Thus, specifying the hierarchical organization of a network not only requires the partitioning of the nodes into different levels but the sequential order of the levels that maximizes the connectivity between adjacent levels must also be identified. To quantify the extent to which a given network exhibits hierarchical organization, we introduce a hierarchy index H, which for a directed, unweighted network whose nodes have been partitioned into a number of sequentially arranged levels is defined as[1]H=1L∑i,jAij−kiin·kjoutL·(δli,lj+1+δli+1,lj),

where Aij represents an element of the adjacency matrix (=1 if there is a directed link from node j to node i, and 0 otherwise), L (=ΣijAij) is the total number of connections in the network, kiin(=ΣjAij) is the in-degree of node i, i.e., the total number of connections received by it and kjout(=ΣiAij) is the out-degree of node j, i.e., the total number of its outgoing connections. The largest magnitude of H is obtained for a partitioning that maximizes the number of connections between adjacent layers in the hierarchical arrangement. This is ensured by performing the summation over only those pairs of nodes that occur in immediately neighboring layers in a given partition via the introduction of the Kronecker delta functions δli,lj+1+δli+1,lj (=1 if levels li and lj to which i and j belong, are adjacent to each other in the sequential arrangement, and =0, otherwise). As kiinkjout/L is the probability of a connection from j to i in a homogeneous random network with the same degree sequence as the network under consideration, the difference with the occurrence frequency in the empirical network (given by the adjacency matrix A) provides a measure of the excess number of links between adjacent layers over that expected by chance. Note that it is expected that H∼0 for a homogeneous, unstructured network. We would like to point out that the concept of hierarchy that we quantify by H is distinct from that of hierarchical modularity, i.e., network structure characterized by internested communities (25–28).

Maximization of the Hierarchy Index.

Having defined the hierarchical arrangement of a network to be a partitioning of the nodes into L sequentially arranged levels that maximizes the hierarchy index H, we require a procedure by which to obtain this arrangement given only the adjacency matrix of a network comprising N nodes. Noting that the analogous task of uncovering the community organization of a network by determining its optimal partitioning that yields the maximum possible value of the modularity measure Q is known to be NP-hard (59, 60), we proceed to obtain an approximate solution to the problem of finding the configuration that maximizes H by using a probabilistic algorithm, specifically, simulated annealing (61) (see Fig. 1 which summarizes the procedure described below). As the configuration specifies not only the node membership of each level but also the sequence in which these levels occur, the heuristic routine for searching the configuration space needs to explore different partitionings of nodes, as well as, alternate arrangements of levels. This is achieved by beginning with an initial configuration that comprises an arbitrary number (typically 5) of sequentially arranged levels with the nodes randomly partitioned between them and then iteratively altering the configuration by performing any one of the following operations at each i) randomly select any one of the N nodes and move it from the level it is occupying to any of the other L−1 levels, ii) create a new level, placed at the end of the existing sequence, comprising any one of the N nodes extracted at random from the L existing levels, iii) merge two levels that are chosen at random from the L levels, iv) exchange the positions in the sequence of any two levels chosen at random from the L levels, and v) split any one of the L levels chosen at random into two, placed adjacent to each other in the sequence. The total number of possible ways in which these operations can be carried out is (N(L−1)+N+LC2+LC2+L=)NL+L2 and at each iteration of the algorithm we choose any one of these with equal probability.

If the hierarchy index computed for the new configuration resulting from the operation carried out at a particular step is higher than the H of the existing configuration, the alteration to the hierarchical structure is accepted. On the other hand, if the change in the hierarchy index ΔH<0, the new configuration is accepted with a probability P∼e−|ΔH|/T, where the parameter T is referred to as temperature in analogy with thermal annealing. In the course of annealing, the temperature is gradually reduced according to a cooling schedule, viz., Tn=T0e−λn, where n is the number of iterations and T0 is the initial temperature. For the results shown in our paper, we have chosen T0=10, λ=2×10−6 and have carried out the annealing for nmax=2×107 iterations which was sufficient for convergence of the process for networks having N∼300 nodes (corresponding approximately to the sizes of the empirical networks we have considered). As the temperature decreases, the system tends to spend longer times in a particular configuration until a new configuration is accepted and we terminate the algorithm if the configuration has not altered in the preceding ncutoff iterations [we have set ncutoff equal to five times the total number of possible operations, viz., 5 (NL+L2)].

The solutions resulting from applying the algorithm on a given network need not be unique, as there could be multiple optimal hierarchical configurations characterized by high values of H. Given this degeneracy, multiple realizations of the hierarchy index maximization process have been carried out ( 62) to construct an ensemble of optimal partitionings for each network (e.g., 103 realizations for the Macaque connectome and 200 for the human and C. elegans connectomes). A network is determined to possess a robust hierarchical organization if the solutions comprising the ensemble are mutually consistent in terms of both the node membership of the different levels, as well as, the sequence in which these levels are arranged. The process by which the level of agreement between the different optimal partitionings of a network is quantified is described below (see section on determining the robustness of the partitions and their sequence). The impact of randomized rewiring that, by design, reduces the hierarchical character of a network, on the hierarchy index H is shown in SI Appendix, Fig. S11.

Benchmarking the Performance of the Algorithm.

In order to show that the maximization of hierarchy index using the procedure outlined above does indeed uncover the inherent hierarchical structure of a network (if any), we test the algorithm on ensembles of random networks whose connection topology has a hierarchical organization by design. To construct a network having a desired extent of hierarchy, we use the ratio h=ρnc/ρcon as a tuning parameter, where ρcon and ρnc represent the density of connections between nodes occurring in consecutive levels in the hierarchical sequence and that between nodes occurring in all other levels (including the same level), respectively. The parameter h can vary over the interval [0,1], with h=1 corresponding to a homogeneous network without any hierarchy, while for h=0, nodes at each level connect only to those at the levels immediately above or below, corresponding to a rigidly hierarchical organization. Apart from the ensemble of hierarchical random networks, we have also considered an additional ensemble of modular hierarchical random networks wherein the network comprises multiple modules or communities, each containing an embedded hierarchical structure.

For the ensemble of hierarchical random networks (SI Appendix, Fig. S1), we generate a benchmark network for a given value of h, having N nodes equally distributed among L levels by linking nodes occurring at different levels with the connection probabilities[2]ρcon=ρnch=k·L2(1−h)(1−L)(N/L)+hLN,

where k is the average degree of the network. For the simulation results reported here, we chose L=4, with N=272 and k=10 (similar to the corresponding values for the empirical networks).

To generate hierarchical modular random networks (SI Appendix, Fig. S3), we assume that the N (=272) nodes of the network are clustered into m (=4) modules. Within each module Lm (=4) hierarchical layers are embedded, such that the network has L=mLm hierarchical layers in total. Apart from the tuning parameter h for the hierarchy, we also use an additional parameter r that specifies the extent of modularity or community organization in the network. It is defined as the ratio of the density of connections between nodes belonging to different modules (ρo) to those occurring in the same module (ρi) (39). The nodes of a network are equally distributed between the modules and are linked according to the connection probabilities[3]ρi=ρor=k(N/m)(1−r)+Nr,

where the average degree k=10. Having determined the mean connection density ρi within a module, we can obtain the connection probability between nodes occurring in the different hierarchical levels that are embedded within each module as[4]ρcon=ρnch=Lm(N/m)ρi2(1−h)(Lm−1)(N/L)+hLm(N/m).

For each class of benchmark networks (hierarchical and modular hierarchical) we generate 20 adjacency matrices for each value of h which logarithmically spans the interval [10−2,1]. The algorithm for maximizing the hierarchy index is applied after randomly permuting the order of the nodes in the adjacency matrix so that its hierarchical structure is no longer apparent. The optimal hierarchical configurations of the network that are obtained from multiple realizations of the H maximization process can then be compared with the original partitioning of the nodes into levels that have been embedded by design. We can quantify how close two different hierarchical decompositions A and B (comprising LA and LB number of levels, respectively) of a network are to each other by computing a similarity score between the sequence of levels {liA}i=1LA and {ljB}j=1LB describing the two decompositions. For this purpose we use the normalized mutual information (62), viz.,[5]Inorm(A,B)=2∑i∑jP(liA,ljB)ln[P(liA,ljB)/P(liA)P(ljB)]−∑iP(liA)lnP(liA)−∑jP(ljB)lnP(ljB),

where P(liX) is the probability that a randomly chosen node lies in level li in partition X∈{A,B}, while P(liA,ljB) is the joint probability that a randomly chosen node belongs to level liA in partition A but occurs in level ljB in partition B (i=1,…,LA, and j=1,…,LB).

Establishing Robust Hierarchical Structure in the Empirical Networks.

As mentioned above, the different hierarchical configurations of a network obtained from multiple realizations of the H maximization algorithm should be similar both in terms of the node composition of their levels, as well as, the sequence in which the levels occur, for any hierarchical organization identified in the network to be robust. The distribution of normalized mutual information Inorm (Eq. 5) between every pair of hierarchical configurations obtained can give us a gross measure for the variability between the solutions obtained from the different realizations. However, in order to quantify the extent to which each node in the empirical network occupies a consistent position in the hierarchical sequence of levels we need to first identify a reference sequence R with which to compare all configurations. To this end, for each hierarchical decomposition A, we compute the mean of the normalized mutual information between it and the decompositions X obtained from all other realizations, viz., I¯norm(A)=⟨Inorm(A,X)⟩X, and choose the configuration that has the maximum value of I¯norm as the reference.

Robustness of nodal composition of the levels.

Upon numbering the levels of the reference decomposition as they occur in sequence from 1,…,LR, a mapping Fl:{liX}i=1LX→{ljR}j=1LR is established between the levels in any hierarchical configuration X obtained from the different realizations and those occurring in R. This is achieved by identifying for each level i∈X the corresponding level in R with which it has maximum overlap. Thus, it allows us to express the identity of the level that a particular node belongs to in a given realization in terms of a standard numbering convention common across all realizations, viz., that of the levels of the reference configuration R. We then compute the fraction of realizations (or trials) ftrls(q,p) in which node p occurs in level q of the reference sequence (see Figs. 2E, 4E, and 5E). For a hierarchical organization identified by the algorithm to be considered robust, the nodes should occur consistently in the same hierarchical level, implying that the distribution of ftrls(q,p) is highly localized. We have ensured that for each of the empirical networks that we have investigated here, the majority of nodes p∈{1,2,…N} satisfy maxqftrls(q,p)>0.5.

Robustness of the sequence of hierarchical levels.

While the composition of each of the network partitions that correspond to the different levels may be consistent across multiple realizations of the hierarchical decompositions, it is possible that the order {1,…,LX} in which the levels occur sequentially in a given realization X may be drastically different from that of the reference sequence R. Therefore, we need to ensure that the sequential arrangement of the partitions is also consistent between the different realizations. For this purpose, we construct another mapping Fs:{liR}i=1LR→{ljX}j=1LX which relates the levels in the reference sequence R to those occurring in the configuration X. In contrast to the mapping Fl, this is achieved by identifying for each level i∈R the corresponding level in X with which it has maximum overlap. Subsequently, we reorder the LR layers of R according to the rank of the layers in X that they map to. This allows us to identify the extent to which the sequential arrangement of levels in R gets rearranged in a realization X. The robustness of the sequential order of the hierarchical levels is quantified by computing the fraction of realizations (or trials) ftrls(i,j) in which layer i of the reference sequence R occurs in position j of the sequence obtained upon reordering according to the mapping with X (see Figs. 2C and 4C, as well as, SI Appendix, Figs. S4 and S10). In the ideal situation, where the hierarchical levels consistently occur in exactly the same sequence across all realizations, ftrls(i,j)=δij, i.e., the Kronecker delta function for i,j=1,…,LR. Note that for the empirical networks investigated here, the matrices representing ftrls have most diagonal entries close to 1 with off-diagonal entries ≪1 indicating that the sequential arrangement of the levels is robust (see SI Appendix, Text and Figs. S12–S14 for alternative measures of the robustness of hierarchical partitioning).

Supplementary Material