## Abstract

The wiring of the brain is organized around a putative unimodal-transmodal hierarchy. Here we investigate how this intrinsic hierarchical organization of the brain shapes the transmission of information among regions. The hierarchical positioning of individual regions was quantified by applying diffusion map embedding to resting-state functional MRI networks. Structural networks were reconstructed from diffusion spectrum imaging and topological shortest paths among all brain regions were computed. Sequences of nodes encountered along a path were then labeled by their hierarchical position, tracing out path motifs. We find that the cortical hierarchy guides communication in the network. Specifically, nodes are more likely to forward signals to nodes closer in the hierarchy and cover a range of unimodal and transmodal regions, potentially enriching or diversifying signals en route. We also find evidence of systematic detours, particularly in attention networks, where communication is rerouted. Altogether, the present work highlights how the cortical hierarchy shapes signal exchange and imparts behaviorally relevant communication patterns in brain networks.

## Author Summary

In the present report we asked how signals travel on brain networks and what types of nodes they potentially visit en route. We traced individual path motifs to investigate the propensity of communication paths to explore the putative unimodal-transmodal cortical hierarchy. We find that the architecture of the network promotes signaling via the hierarchy, suggesting a link between the structure and function of the network. Importantly, we also find instances where detours are promoted, particularly as paths traverse attention-related networks. Finally, information about hierarchical position aids navigation in some parts of the network, over and above spatial location. Altogether, the present results touch on several emerging themes in network neuroscience, including the nature of structure-function relationships, network communication and the role of cortical hierarchies.

## INTRODUCTION

Adaptive behavior requires transmission of information between neuronal populations. The architecture of white matter networks supports an array of signal propagation patterns, linking sensation, cognition, and action (Avena-Koenigsberger, Misic, & Sporns, 2018). Brain networks, reconstructed from multiple species and at multiple spatial scales, possess multiple nonrandom attributes that make such flexible communication possible, including near-minimal path length and high clustering (Hilgetag & Kaiser, 2004; Kaiser & Hilgetag, 2006), as well as assortative community structure (Betzel, Medaglia, & Bassett, 2018) and a densely interconnected core (van den Heuvel, Kahn, Goñi, & Sporns, 2012).

How do signals traverse neural circuits and what types of neuronal populations do they encounter along the way? The sequence of neurons and neuronal populations that a signal passes through presumably transforms the nature of the signal itself and its downstream effect (Amico et al., 2019; Avena-Koenigsberger et al., 2019; Graham & Rockmore, 2011; Mišić, Sporns, & McIntosh, 2014; Mišić, Sporns, & McIntosh, 2014; Seguin, Razi, & Zalesky, 2019; Worrell, Rumschlag, Betzel, Sporns, & Mišić, 2017). For example, signals exchanged between closely clustered and functionally aligned populations may be relatively unchanged, whereas signals exchanged between anatomically and functionally distant populations may be enriched or diversified (Bertolero, Yeo, & D’Esposito, 2017). A simple way to infer potential signal trajectories in a network is the topological shortest path (hereafter simply referred to as a “path”; Avena-Koenigsberger et al., 2017; van den Heuvel et al., 2012). For many classes of networks, including brain networks, decentralized communication mechanisms may also take advantage of shortest paths without any knowledge of the global topology, including diffusion (Goñi et al., 2014) and navigation (Seguin, van den Heuvel, & Zalesky, 2018). Thus, paths connecting pairs of nodes trace out unique motifs along their trajectory, meaning that the nature of communication between any two regions is subject to the underlying structure (Avena-Koenigsberger et al., 2019; Graham & Rockmore, 2011; Mišić et al., 2015).

Signal propagation is likely to be constrained by the hierarchical organization of cortical circuits. Evidence from classical anatomy and modern neuroimaging points to a continuous sensory-transmodal hierarchy, spanning unimodal to transmodal cortex (Goulas, Zilles, & Hilgetag, 2018; Hilgetag & Goulas, 2020; E. Jones & Powell, 1970; Mesulam, 1998). This continuous axis or gradient can be observed in the functional architecture of the cortex (Margulies et al., 2016), running parallel to gradients in intracortical myelin (Huntenburg et al., 2017; Paquola et al., 2019), cortical thickness (Wagstyl, Ronan, Goodyer, & Fletcher, 2015), gene transcription (Burt et al., 2018; Fulcher, Murray, Zerbi, & Wang, 2019), excitation-inhibition ratios (Wang, 2020), and intrinsic temporal timescales (Kiebel, Daunizeau, & Friston, 2008; Murray et al., 2014). The influence of these multimodal gradients on signaling and communication in structural networks is a key question in systems neuroscience (Seguin et al., 2019; Vázquez-Rodríguez et al., 2019).

Here we investigate how the functional hierarchy shapes the propagation of signals. We reconstruct paths on structural networks and trace their trajectories through the unimodal-transmodal gradient. We find that the hierarchical organization of the cerebral cortex constrains path trajectories, such that many paths follow a canonical bottom-up (ascending the hierarchy) or top-down (descending the hierarchy) trajectory. Importantly, we find that paths may potentially reverse direction in attention networks. Altogether, we find that the hierarchical organization of cortical circuits imposes a communication space on the structural network, potentiating some types of signal propagation patterns while attenuating others.

## RESULTS

The results are organized as follows. We first develop a methodology to trace signal trajectories through the putative unimodal-transmodal hierarchy. We then investigate the extent to which signal flows conform to the hierarchical organization of the cortex, and instances where they diverge. Finally, we consider whether information about hierarchical position is sufficient to sustain a decentralized navigation-like communication process. Data sources include the following (see the Methods section for detailed procedures):

- •
*Structural connectivity*. Structural and functional connectivity were derived from*N*= 66 healthy participants (source: Lausanne University Hospital). Structural connectivity was reconstructed using diffusion spectrum imaging and deterministic streamline tractography. A consistency- and length-based procedure was then used to assemble a group-representative weighted structural connectivity matrix (Betzel, Griffa, Hagmann, & Mišić, 2018; Mišić et al., 2018; Mišić et al., 2015). - •
*Functional connectivity*. Functional connectivity was estimated in the same individuals using resting-state functional MRI (rs-fMRI). A functional connectivity matrix was constructed using pairwise Pearson correlations among regional time courses. A group-average functional connectivity matrix was then estimated as the mean connectivity of pairwise connections across individuals.

We trace path motifsbetween all possible source-target node pairs on the weighted structural network (Figure 1; for a conceptually similar approach, see van den Heuvel et al., 2012). We label nodes according to two different nomenclatures: hierarchical position and intrinsic network affiliation (Yeo et al., 2011). Hierarchical position is defined as the first principal connectivity gradient of the diffusion map embedding over the FC matrix (Margulies et al., 2016; see Methods). The continuous embedding vector spans a putative hierarchy, where lower values correspond to unimodal regions and greater values correspond to transmodal regions. We use the empirical cumulative distribution function of the first gradient to bin nodes into 10 classes of equal size. We enumerate the classes from 1 to 10, where 1 corresponds to unimodal cortex and 10 to transmodal cortex.

### Path Motifs Follow Hierarchies

We first investigate how path motifs map onto the putative unimodal-transmodal hierarchy. For a given source and target class, we consider all possible paths between the constituent source nodes and target nodes. We then compute the mean hierarchical position of every step encountered along the paths. Figure 2 shows the mean path motifs originating from a low-(class 2), intermediate-(class 6), and high-level (class 9) source region, to the same three regions as targets. Colors distinguish paths of different lengths.

In general, path motif shape depends on the relative hierarchical position of the source and target nodes (Figure 2; rows and columns, respectively). For most paths, motifs sequentially transition through the hierarchy, but there exist systematic differences in the nature of the transitions. Paths traversing a larger difference in the hierarchy (e.g., from class 2 to class 9) tend to follow a more monotonic trajectory. Conversely, when the source and the target nodes occupy the same or neighboring positions in the hierarchy, paths are more likely to follow a U-shape, effectively taking detours to intermediate parts of the hierarchy (e.g., when source and target is class 2). These trajectories are in contrast to trajectories recovered from null networks with randomized topology (Supporting Information Figure S1) and randomized hierarchy labels (Supporting Information Figure S2). It is noteworthy that the effect sizes for the latter are considerably smaller, suggesting that path motifs are shaped both by the topology of the network, and by the spatial structure of the unimodal-transmodal hierarchy.

### Inflection Points in Communication Flow

In the previous section we observed that anatomical connectivity fundamentally shapes how the unimodal-transmodal hierarchy is traversed, promoting some transitions while attenuating others. To investigate how these “pushing” and “pulling” forces shape the communication landscape, we next consider path trajectories at the nodal level. Specifically, we study how the orientation of the flows changes along the course of the journey towards the target, which we quantify by the slope of paths through the hierarchy.

For a given node, we calculate the mean slope of all paths as they pass through that node (Figure 3A). If, on average, the slope of the paths when passing through a node is positive, this suggests that the transmission is ascending through the hierarchy, from unimodal towards transmodal cortex. Conversely, nodes with negative slopes direct information flow towards areas lower in the cortical hierarchy. Note that, in general, a node could participate in both ascending and descending paths, and this dependent measure reflects the mean flow of information through that node.

We first note that the mean slope for a given node is negatively correlated with
the node’s position in the unimodal-transmodal hierarchy (Figure 3B; *r* = −0.72, *P*_{spin} = 0.05). In other words,
nodes that occupy higher positions in the hierarchy tend to direct signal
traffic towards nodes lower in the hierarchy, and vice versa. This is consistent
with the intuition developed in the previous subsection. Figure 3B shows that areas that exhibit mainly positive mean
slope (i.e., direct information to ascend the hierarchy) are the supplementary
motor area, somatomotor cortex, and visual cortex. Areas with negative slope
(i.e., directing information to descend the hierarchy) are prefrontal cortex,
posterior parietal cortex, auditory cortex, and inferotemporal cortex.
Stratifying these nodes by their membership in intrinsic networks (Yeo et al., 2011), we find mean positive slopes
for the visual, somatomotor, dorsal attention, ventral attention and
frontoparietal networks, and mean negative slopes for the limbic and default
mode networks (Figure 3C).

The slope of a path traversing a node also allows us to identify areas that redirect flow direction and promote detours. As we follow a path trajectory, we look for local extremum nodes that reside between slopes with different signs, and tag them as turning points. Depending on the type of extremum, we name the turning points as “turning up” points (local minima) or “turning down” points (local maxima; Figure 3A). For example, a turning down point is a node that occupies a relatively higher position in the hierarchy and connects two lower level brain areas (see Figure 3A). We first stratify nodes by their membership in intrinsic networks and compute the mean turning point probability for each network. Figure 3D shows that networks with the greatest probability of turning up paths (i.e., redirecting them to ascend the hierarchy) are the dorsal attention and ventral attention, whereas the default mode network has the greatest turning down probability. At the regional level, regions with the greatest probability of turning up the paths tend to be in attention-related networks, including the supplementary motor area and posterior parietal cortex (Figure 3E). Conversely, superior and dorsolateral prefrontal, inferotemporal, and lateral temporal cortex are the most probable turning down points (Figure 3F).

### Temporal Evolution of Communication Flow

How does the hierarchical organization of the brain shape communication across time? Given that most communication paths conform to the hierarchical organization of the network, we next ask whether the hierarchy imparts memory on communication processes by exerting influence on the path trajectories. To address this question, we consider the temporal evolution of communication patterns, envisioning the sequence of nodes traversed along a path as a time series.

We explore how the position of a walker traversing the path depends on the
positions it occupied previously in its trajectory. Specifically, we compute the
probability of going from a node that belongs to the hierarchy level *i* to a node that belongs to hierarchy level *j* in one step as a function of the position over the path
(1-hop transitions; Figure 4A, left). To
quantify whether the hierarchical position of a walker depends on its previous
hierarchical positions (multi-hop transitions), we measure the transition
probability of occupying hierarchy level *i* at step *t*, given that the walker occupied hierarchy level *j* at step *t* − *k*,
changing *k* from 1 to the length of the path (multi-hop
transitions; Figure 4A, right).

Figure 4B (top) shows 1-hop transition
probabilities, as in a first-order Markov process. Nodes are stratified
according to their position in the hierarchy, with source nodes in the rows and
target nodes in the columns. Mean transition probabilities are greatest over the
diagonal, favoring the transmission of a walker to neighboring positions in the
hierarchy, repeating the theme from the previous two subsections. Figure 4 (bottom) shows average probabilities
for transitions (memories) over 2 steps or more. For memories 2, 3, and 4 steps
away, transitions become more uniform, meaning that the probability of occupying
current position *j* does not depend greatly on the position it
occupied 2, 3, or 4 steps before. For greater memory values of 5, 6, and 7,
there is an emergence of transitions between lower and intermediate, and higher
and intermediate hierarchy levels, with almost zero probability of a transition
between lower and higher levels. Altogether, we find that most 1-hop transitions
follow the hierarchy, but that communication over longer trajectories is biased
towards some levels of the hierarchy and away from others, particularly if the
starting point is at intermediate levels. In other words, the nodes visited by a
walker earlier in the trajectory may exert influence on transition probabilities
later in the trajectory.

### Navigation via Hierarchies

Given that communication paths closely align with the hierarchy of the network, we finally ask whether it is possible to recapitulate the path architecture of the network by following the hierarchy. We focus on navigation, a decentralized communication mechanism in which a signal is forwarded to the connected neighbor that is closest in some distance to the target. This distance is defined with respect to an underlying metric space, with the simplest such space being the three-dimensional space that nodes are physically embedded in. For example, previous work has demonstrated that it is possible to recapitulate the shortest path architecture by forwarding signals to nodes that are physically closest to the target node (Seguin et al., 2018). Decentralized mechanisms such as navigation are intuitively appealing as they do not assume that signals or nodes possess knowledge of the global topology (Avena-Koenigsberger et al., 2018; Seguin et al., 2018).

*S*

_{R}). Given the importance of spatial embedding, we identify regions for which navigation success improves when hierarchy information is taken into account rather than only spatial factors. To operationalize a node’s proximity to the target we use a linear combination of Euclidean distance in three-dimensional physical space and distance in “hierarchy space” weighted by the parameter

*β*as the following:

*d*(

*i*,

*j*) is the combined distance between nodes

*i*and

*j*, (

*x*

_{i},

*y*

_{i},

*z*

_{i}) gives the position of node

*i*in three-dimensional Euclidean space, and

*h*

_{i}gives the position of node

*i*in one-dimensional hierarchical space. For each pair of nodes we measure the success ratio as a function of

*β*, tuning

*β*from 0 to 1, and find the

*β*that maximizes the navigation success. When

*β*is valued close to 1, paths originating from the node are better recovered using spatial proximity compared with hierarchical proximity; the opposite is true when

*β*is valued close to 0 (Figure 5A; see the Methods section for more detail).

Figure 5B shows the distribution of the *β* parameters that maximize navigation success for
each source-target pair. The distribution is bimodal, with prominent peaks at
the extremes, suggesting that most nodes have a strong preference for either
hierarchical or spatial navigation. Stratifying nodes by membership in intrinsic
networks, we find that each network exhibits a unique fingerprint, with some
showing a preference for spatial navigation (frontoparietal), others for
hierarchical navigation (default), and others a mix between the two (dorsal
attention). Figure 5C shows that parts of
the visual system, lateral temporal cortex, and dorsolateral prefrontal cortex
exhibit a strong preference for spatial navigation (red; *β* > 0.8), while medial parietal cortex,
medial prefrontal cortex, and left temporo-parietal cortex exhibit a strong
preference for hierarchical navigation (blue; *β* <
0.2).

### Relation With Simple Measures

In this work we derived four dependent measures based on the concept of path
motifs (slope, tuning point up/down probability, and navigation preference). For
completeness, we assess the extent to which these node-level variables can be
related to simpler measures. Supporting
Information Figure S4 shows linear regressions comparing each of the
four path motif measures (rows) with simpler network measures computed from
structural and functional connectivity matrices (columns). From the structural
network we compute betweenness, closeness, clustering, degree, and mean edge length. All measures except
degree are computed on the weighted network. From the functional network we
compute strength and participation
coefficient (relative to the intrinsic network partition reported by
Yeo et al., 2011). As expected, we find
weak to moderate correlations with path-based measures (betweenness, closeness)
and with degree, consistent with the notion that most centrality measures are
correlated with each other (Oldham et al., 2019). In addition, we find a positive correlation between
participation and mean slope (*r* = 0.52), suggesting that nodes
with more diverse connection profiles are more likely to direct communication
towards the apex of the unimodal-transmodal hierarchy. In sum, we find that the
four path motif-derived measures are correlated with some simpler measures, but
cannot be wholly predicted from any one such measure.

## DISCUSSION

In the present report we asked how signals travel on brain networks and what types of nodes they potentially visit en route. We traced individual path motifs to investigate the propensity of communication paths to explore the putative unimodal-transmodal cortical hierarchy. We find that the architecture of the network promotes signaling via the hierarchy, suggesting a link between the structure and function of the network. Importantly, we also find instances where detours are promoted, particularly as paths traverse attention-related networks. Finally, information about hierarchical position aids navigation in some parts of the network, over and above spatial location. Altogether, the present results touch on a number of emerging themes in network neuroscience, including the nature of structure-function relationships, network communication, and the role of cortical hierarchies.

The most prominent observation is that most paths closely follow the cortical hierarchy, traveling from unimodal to transmodal cortex and vice versa (Figure 2). This is reminiscent of the notion that much of signal traffic follows a sequential bottom-up or top-down trajectory, potentiating direct stimulus-response patterns (Worrell et al., 2017). At the same time, the architecture of the network occasionally serves to redirect signal traffic and promote detours. In particular, the default network appears to be the most likely mediator between areas lower in the cortical hierarchy, while the dorsal attention and ventral attention networks act as mediators between areas higher in the hierarchy (Figure 3). The finding that attention networks are anatomically positioned to redirect signal traffic resonates with modern theories of how attention and control networks shape fluid transitions between segregated and integrated states, guiding adaptive reconfiguration during rest and task (de Pasquale et al., 2012; Fair et al., 2007; Mišić, Fatima et al., 2014; Mohr et al., 2016; Shine et al., 2016).

These results suggest that the architecture of the network promotes behaviorally relevant communication patterns, and that the functional properties of individual areas are fundamentally related to their anatomical network embedding. Indeed, multiple studies point to the idea that the topology of brain networks endows individual regions with specific functional attributes. For example, functional properties depend on local connectional profiles, including motif composition (Gollo, Mirasso, Sporns, & Breakspear, 2014), asymmetry (Knock et al., 2009), length (Vázquez-Rodríguez et al., 2019), and weight distribution (Melozzi et al., 2019). At the global level, anatomical segregation, most prominently observed in the unimodal visual and somatosensory cortices, promotes specialized processing. Conversely, polysensory association cortex is anatomically better integrated in the connectome, potentially allowing information to be sampled from multiple parts of the network (Mišić et al., 2015; Mišić, Goñi, Betzel, Sporns, & McIntosh, 2014; van den Heuvel et al., 2012). A recent report demonstrated that regions at the top of the hierarchy are better positioned to act as “receivers” of information, while regions at the bottom are better positioned to act as “senders” (Seguin et al., 2019). The present results build on this literature, showing that the default and attentional networks present links between parts of the cortical hierarchy. One prediction, to be tested in future studies, is that individual differences in network architecture may potentiate some signaling patterns while attenuating others. Wiring patterns in which the entire sensory-fugal hierarchy is more easily traversed may promote integration between perception, cognition, and action. Thus, individual differences in anatomical connectivity may allow signals to “ride” the hierarchy more easily, resulting in better cognitive performance and more adaptive behavior.

In addition to network topology, we also highlight the contribution of spatial structure in shaping path motifs. Comparison with rewired networks with geometric constraints yielded large deviations from empirically derived path motifs (Supporting Information Figure S1), but comparison with label-rotated null models yielded more modest effect sizes (Supporting Information Figure S2). This suggests that path trajectories through the unimodal-transmodal hierarchy are driven not only by network topology, but also by the strong spatial structure of the hierarchy. Indeed, multiple studies have demonstrated that the probability and strength of connectivity between areas is anticorrelated with their spatial proximity (Betzel & Bassett, 2018; Ercsey-Ravasz et al., 2013; Horvát et al., 2016; Mišić, Fatima et al., 2014; Roberts et al., 2016; Stiso & Bassett, 2018). As a result, the connectional fingerprint of a given area tends to follow the distribution of spatial distances with its neighbors, such that hierarchical position and spatial separation are closely intertwined (Oligschläger et al., 2017; Oligschläger et al., 2019).

In shaping communication patterns, network architecture may also impart memory on signal traffic, such that transitions depend not only on the current position of the signal, but also on positions they previously occupied in the hierarchy in their journey (Figure 4). We find that most transitions between nodes that occupy neighboring positions in the hierarchy are memoryless, but that transitions across disparate levels are not, with areas intermediate in the hierarchy helping to mediate longer memory communication. The phenomenon of network structure imposing non-Markovian network flows is also observed in other complex systems, such as air passenger flows and journal citation flows (Rosvall, Esquivel, Lancichinetti, West, & Lambiotte, 2014). In the brain, the functional consequences of this phenomenon may be the well-studied hierarchical organization of timescales and temporal receptive windows. Numerous studies point to the idea that information accumulates over time across the cortical hierarchy, such that the temporal dynamics of higher order regions unfold over slower timescales (Kiebel et al., 2008; Murray et al., 2014), manifesting as a preference for long-range contextual information (Honey et al., 2012). Understanding the precise link between network structure and the hierarchy of intrinsic timescales remains a major question for future research (Demirtaş et al., 2019; Gollo, Roberts, & Cocchi, 2017).

More generally, these results open new questions about how broad spatial gradients in synaptic connectivity, cytoarchitecture, and molecular composition interact with macroscale network topology (Margulies et al., 2016; Fulcher et al., 2019; Vázquez-Rodríguez et al., 2019; Wang, 2020). Current graph models of brain networks assume that all nodes are the same, but as signals propagate through the network, they pass through a series of heterogeneous neural circuits and populations (Amico et al., 2019; Suarez, Markello, Betzel, & Misic, 2020). Each stage in the trajectory may entail transformations that modern methods in network neuroscience do not consider. For example, the majority of path-based metrics consider the total lengths of paths between areas, but not the identity of nodes traversed during the path. By drawing path motifs through maps annotated by molecular and cellular data, the present methodology permits closer investigation into how local attributes of nodes may influence communication in the network.

Note that the present method only traces out the shortest paths in the network, but this does not preclude the possibility that communication takes place via mechanisms that are unaware of the shortest paths in the network (Avena-Koenigsberger et al., 2018; Avena-Koenigsberger et al., 2017; De Domenico, 2017; Suarez et al., 2020). Several recent reports point to diffusion-like and navigation-like processes as potentially more biologically realistic alternatives (Goñi et al., 2014; Gulyás, Bíró, Kőrösi, Rétvári, & Krioukov, 2015; Mišić, Sporns, & McIntosh, 2014; Seguin et al., 2018), as they do not assume that signals possess knowledge of the global topology. At the same time, multiple studies suggest that shortest paths in brain networks are readily accessible by both diffusion (Goñi et al., 2014) and spatial navigation (greedy routing; Seguin et al., 2018), without knowledge of the global topology. In other words, there are likely to be some similarities between random walk motifs and path motifs, but this needs to be investigated in more detail. We envision that future studies will consider diffusion and navigation trajectories, analogous to the approach we took with shortest paths.

We close by noting important methodological considerations. Although the present networks are derived from a consensus of 66 participants with high-quality imaging (Betzel, Griffa, Hagmann, & Mišić, 2018), there are several limitations. First, structural networks were reconstructed using diffusion-weighted MRI and computational tractography, a technique that results in systematic false positives and false negatives (de Reus & van den Heuvel, 2013; Maier-Hein et al., 2017; Thomas et al., 2014). Second, the technique cannot be used to resolve the direction of white matter projections, which means that some paths recovered from the network may not exist. Third, the present reconstruction only includes cortical regions, leaving out important topological contributions from the subcortex and cerebellum. Network communication is undoubtedly shaped by both sets of structures, and future studies should consider subcortical-cortical and cerebellar-cortical signal traffic.

In summary, we develop a simple framework to trace communication patterns in brain networks. We show that the putative unimodal-transmodal hierarchy shapes the propagation of signals across the connectome. The present work highlights the importance of considering sequences of nodes encountered during signaling, and the role they might play in network-wide communication.

## METHODS

### Data Acquisition

A total of *N* = 66 healthy young adults (16 females, 25.3
± 4.9 years old) were scanned at the Department of Radiology, University
Hospital Center and University of Lausanne. The scans were performed in 3-Tesla
MRI scanner (Trio, Siemens Medical, Germany) using a 32-channel head coil. The
protocol included (a) a magnetization-prepared rapid acquisition gradient echo
(MPRAGE) sequence sensitive to white/gray matter contrast (1mm in-plane
resolution, 1.2mm slice thickness), (b) a diffusion spectrum imaging (DSI)
sequence (128 diffusion-weighted volumes and a single b0 volume, maximum b-value
8,000/mm^{2}, 2.2 × 2.2 × 3.0 mm voxel size), and (c)
a gradient echo EPI sequence sensitive to BOLD contrast (3.3mm in-plane
resolution and slice thickness with a 0.3mm gap, TR 1,920 ms, resulting in 280
images per participant). Participants were not subject to any overt task demands
during the fMRI scan.

### Structural Network Reconstruction

Gray matter was parcellated into 68 cortical nodes according to the Desikan-Killiany atlas (Desikan et al., 2006). These regions of interest were then further divided into 1,000 approximately equally sized nodes (Cammoun et al., 2012). Structural connectivity was estimated for individual participants using deterministic streamline tractography. The procedure was implemented in the Connectome Mapper Toolkit (Daducci et al., 2012), initiating 32 streamline propagations per diffusion direction for each white matter voxel. Structural connectivity between pairs of regions was defined as the number of streamlines normalized by the mean length of streamlines and mean surface area of the two regions, termed fiber density (Hagmann et al., 2008). This normalization compensates for the bias towards longer fibers during streamline reconstruction, as well as differences in region size.

To mitigate concerns about inconsistencies in reconstruction of individual
participant connectomes (Thomas et al., 2014; D. Jones, Knösche, & Turner, 2013), as well as the sensitive dependence of network
measures on false positives and false negatives (Zalesky et al., 2016), we adopted a group-consensus
approach (Betzel, Griffa, Hagmann, & Mišić, 2018; de Reus & van den Heuvel, 2013; Roberts, Perry, Roberts,
Mitchell, & Breakspear, 2017).
In constructing a consensus adjacency matrix, we sought to preserve (a) the
density and (b) the edge length distribution of the individual participant
matrices (Betzel et al., 2016; Betzel,
Griffa, Hagmann, & Mišić, 2018; Mišić et al., 2015). We first collated the extant edges in the individual
participant matrices and binned them according to length. The number of bins was
determined heuristically, as the square root of the mean binary density across
participants. The most frequently occurring edges were then selected for each
bin. If the mean number of edges across participants in a particular bin is
equal to *k*, we selected the *k* edges of that
length that occur most frequently across participants. To ensure that
interhemispheric edges are not underrepresented, we carried out this procedure
separately for inter- and intrahemispheric edges. The binary density for the
final whole-brain matrix was 2.17%. The weight associated with each edge
was then computed as the mean weight across all participants.

### Functional Network Reconstruction

Functional MRI data were preprocessed using procedures designed to facilitate subsequent network exploration (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012). FMRI volumes were corrected for physiological variables, including regression of white matter, cerebrospinal fluid, as well as motion (three translations and three rotations, estimated by rigid body coregistration). BOLD time series were then subjected to a low-pass filter (temporal Gaussian filter with full width half maximum equal to 1.92 s). The first four time points were excluded from subsequent analysis to allow the time series to stabilize. Motion “scrubbing” was performed as described by Power and colleagues (Power et al., 2012). The data were parcellated according to the same atlas used for structural networks (Cammoun et al., 2012). Individual functional connectivity matrices were defined as zero-lag Pearson correlation among the fMRI BOLD time series. A group-consensus functional connectivity matrix was estimated as the mean connectivity of pairwise connections across individuals.

### Diffusion Map Embedding

Diffusion map embedding is a nonlinear dimensionality reduction algorithm
(Coifman et al., 2005). The algorithm
seeks to project a set of embeddings into a lower dimensional Euclidean space.
Briefly, the similarity matrix among a set of points (in our case, the
correlation matrix representing functional connectivity) is treated as a graph,
and the goal of the procedure is to identify points that are proximal to one
another on the graph. In other words, two points are close together if there are
many relatively short paths connecting them. A diffusion operator, representing
an ergodic Markov chain on the network, is formed by taking the normalized graph
Laplacian of the matrix. The new coordinate space is described by the
eigenvectors of the diffusion operator. We set the diffusion rate *α* = 1 and the variance of the Gaussian used in
affinity computation *σ* = 1. The procedure was
implemented using the Dimensionality Reduction Toolbox (https://lvdmaaten.github.io/drtoolbox/; van der Maaten, Postma,
& van den Herik, 2009).

### Shortest Path Retrieval

Structural connectivity was encoded as an undirected weighted graph *G* ≡{*V*,*W*} composed of nodes *V* =
{*v*_{1},*v*_{2},
… *v*_{n}} and a matrix of
fiber density values *W* =
[*w*_{ij}], valued on the
interval [0,1]. To recover shortest paths, we first define a topological
distance measure. The weighted adjacency matrix was transformed from a
connection weights to connection lengths matrix using the transform *L* = −*log*(*W*), such
that connections with greater weights are mapped to shorter lengths (Goñi
et al., 2014). Note that other
transformations are also possible, including *L* =
1/*W*. The drawback of this inverse transform is that it
generates highly skewed distributions of lengths *L*. As a
result, a small number of connections are valued much more than the rest, and
they are disproportionately represented in shortest paths. The logarithmic
transform controls for this, yielding log-normal distributions of *L* (Avena-Koenigsberger et al., 2017). Weighted shortest paths were recovered using
the Floyd-Warshall algorithm (Floyd, 1962; Roy, 1959; Warshall, 1962; *Brainconn* Python Toolbox). Note that in many types of networks there may exist multiple
shortest paths between two nodes (edge-disjoint or not); in our network this was
not the case as we computed weighted shortest paths, yielding unique paths
between all source-target pairs.

### Slopes and Inflection Points

A shortest path is defined as a sequence of nodes *v*_{1},*v*_{2},…,*v*_{n} where *v*_{i} is the node occupied at
the *i*th step of a path of length *n*. Each node
has an associated hierarchical position value *h*(*v*_{i}), so the
slope at node *v*_{i} over one path is $svi=h(vi+1)\u2212h(vi)$. The average slope will be the average of $svi$ across all paths that node *v*_{i} participates in, except
in cases where the node is the initial source or final target (i.e., the most
extreme point in the path). As an example, let us consider two paths, 1 :
{*A*,*B*,*C*,*D*},
and 2 :
{*D*,*C*,*B*,*A*}.
The slope of the node C for the first path will be $sC1=h(D)\u2212h(C)$, while for the second path it will be $sC2=h(B)\u2212h(C)$. The average slope of the node C will be *s*_{C} =
(*h*(*D*) + *h*(*B*) −
2*h*(*C*))/2.

### Transition Probabilities

**T**as a function of the position

*t*, where

*t*= 1,…,

*Diam*(

*G*

_{SC}) − 1. For each position there was one matrix

**T**(

*t*) defined as

*i*and

*j*represent the hierarchy bins, thus (

*i*,

*j*) ∈{

*h*

_{1},…,

*h*

_{10}}. The expression #{

*sp*: ∘} represents the number of shortest paths, from the set of all shortest paths, that satisfy the condition ∘.

*k*, denoted as

**M**(

*k*), where

*k*= 1,…,

*Diam*(

*G*

_{SC}) − 1. For each hop length

*k*the matrix

**M**(

*k*) was defined as

*i*,

*j*∈{

*h*

_{1},…,

*h*

_{10}},

*θ*=

*length*(

*sp*) −

*k*,

*θ*∈ [1,

*Diam*(

*G*

_{SC}) − 1].

### Network Navigation

*i*to target node

*j*. The agent has no knowledge of the global topology; instead, they hop towards neighbors who are closest to the target node in some underlying metric space. Across all source-target pairs, we measure the proportion of paths that are successfully recovered (success ratio;

*S*

_{R}; Seguin et al., 2018). To operationalize a node’s proximity to the target we use a linear combination of the Euclidean distances in three-dimensional physical space and in hierarchy space weighted by the parameter

*β*as

*d*

_{i,j}is the combined distance between nodes

*i*and

*j*, (

*x*

_{i},

*y*

_{i},

*z*

_{i}) is a node’s position in Euclidean space, and

*h*

_{i}is a node’s position in hierarchy space. Both values were normalized to lie in the interval [0,1]. In the context of neural communication, spatial navigation could simply reflect a cost-minimization strategy, as signals are forwarded to the neighbor proximal to the target. Since hierarchical position is estimated using diffusion map embedding applied to functional connectivity, hierarchy distance can be interpreted as distance in functional similarity. More generally, as the unimodal-transmodal hierarchy is observed for a number of local properties, hierarchy distance could represent a mechanism where signals are forwarded to populations with similar intrinsic electrophysiological rhythms, similar receptor distributions, or similar molecular, or cytoarchitectonic properties.

Given that navigation is a deterministic model, we calculated, for each node as a
source, the success ratio as a function of *β*. The
resulting curves showed a global trend preferring Euclidean distance, with
optimal *β* values close to 0.8, consistent with previous
reports (Seguin et al., 2018). In each
node, however, there exists substantial variance across *β* values, showing a changing preference deviating
from global trend. To better capture this preference for individual nodes, we
detrend the mean success ratio and select for each source node the *β* that maximizes the detrended success ratio (Supporting Information Figure S5).

### Null Models

The critical question underlying all reported analyses is the link between structural connectivity and hierarchical position. To assess this question, we used two null models that systematically disrupt the relationship between network topology and hierarchical labels. The first null model rewires the structural network while preserving the edge lengths and hierarchical labels (Betzel & Bassett, 2018; Gollo et al., 2018; Roberts et al., 2016). Edges were first binned according to Euclidean distance. Within each length bin, pairs of edges were then selected at random and swapped (Betzel & Bassett, 2018). The procedure was repeated 2,000 times, generating a population of rewired structural networks that preserve the degree sequence of the original network and approximately preserve the edge length distribution of the original network.

The second null model permutes hierarchical labels, but preserves their spatial
autocorrelation (Alexander-Bloch et al., 2018). We first created a surface-based representation of the
Cammoun atlas on the FreeSurfer *fsaverage* surface using the
Connectome Mapper Toolkit (https://github.com/LTS5/cmp; Daducci et al., 2012). We used the spherical projection of the *fsaverage* surface to define spatial coordinates for each
parcel by selecting the vertex closest to the center-of-mass of each parcel. The
resulting spatial coordinates were used to generate null models by applying
randomly sampled rotations and reassigning node values based on the closest
resulting parcel (2,000 repetitions). The rotation was applied to one hemisphere
and then mirrored to the other hemisphere.

## AUTHOR CONTRIBUTIONS

Bertha Vázquez-Rodríguez: Methodology; Formal analysis; Writing - Original Draft. Zhen-Qi Liu: Methodology; Formal analysis; Writing - Original Draft. Patric Hagmann Data curation. Bratislav Misic: Conceptualization; Supervision; Writing Review & Editing.

## FUNDING INFORMATION

Bratislav Misic, Natural Sciences and Engineering Research Council of Canada (NSERC Discovery Grant), Award ID: 017-04265. Bratislav Misic, Canada Research Chairs Program. Bratislav Misic, Fonds de recherche du Qubec - Sant (Chercheur Boursier).

## DATA AVAILABILITY

The processed dataset (structural and functional matrices) is available at https://doi.org/10.5281/zenodo.2872624 (Griffa, Alemán-Gómez, & Hagmann, 2019).

## SUPPORTING INFORMATION

Supporting information is available at https://doi.org/10.1162/netn_a_00153.

## ACKNOWLEDGMENTS

The authors thank Dr. Alessandra Griffa for collecting, preprocessing, and sharing the neuroimaging dataset. We acknowledge the Department of Psychiatry of Lausanne University Hospital and particularly Professor Philippe Conus, Professor Kim Do Cuenod, Raoul Jenni, and Martine Cleusix for having helped with the recruitment process of the study volunteers. The authors thank Ross Markello, Golia Shafiei, Vincent Bazinet, Laura Suarez, and Justine Hansen for helpful comments on the manuscript. This research was undertaken thanks in part to funding from the Canada First Research Excellence Fund, awarded to McGill University for the Healthy Brains for Healthy Lives initiative.

## TECHNICAL TERMS

- Shortest path:
The minimum contiguous set of edges between two nodes.

- Diffusion:
A decentralized communication mechanism in which signals are randomly forwarded among nodes in a network.

- Navigation:
A decentralized communication mechanism in which signals are forwarded to nodes that are closest to the target node in some underlying metric space, such as Euclidean space or hierarchical space.

- Unimodal-transmodal hierarchy:
A continuous axis separating unimodal sensory-motor cortex from transmodal association cortex.

- Path motif:
The sequence of node annotations encountered on a shortest path between two nodes.

- Inflection point:
Points along a path where the path motif changes slope, indicating a change of direction, such as first descending and then ascending unimodal-transmodal hierarchy, or vice versa.

- Intrinsic networks:
Subnetworks of brain areas with coherent time courses, identified by clustering, independent component analysis, or community detection.

- Degree:
The number of edges or connections in which a node participates.

- Participation coefficient:
A measure of how evenly distributed a node’s connections are among the network’s communities (1 = maximum, 0 = minimum).

- Transition probability:
The probability that a signal is forwarded from a source node to a target node.

## REFERENCES

## External Supplements

## Author notes

Competing Interests: The authors have declared that no competing interests exist.

These authors contributed equally to this work.

Handling Editor: Alex Fornito