Structural Features of the Human Connectome That Facilitate the Switching of Brain Dynamics via Noradrenergic Neuromodulation

The structural connectivity of human brain allows the coexistence of segregated and integrated states of activity. Neuromodulatory systems facilitate the transition between these functional states and recent computational studies have shown how an interplay between the noradrenergic and cholinergic systems define these transitions. However, there is still much to be known about the interaction between the structural connectivity and the effect of neuromodulation, and to what extent the connectome facilitates dynamic transitions. In this work, we use a whole brain model, based on the Jasen and Rit equations plus a human structural connectivity matrix, to find out which structural features of the human connectome network define the optimal neuromodulatory effects. We simulated the effect of the noradrenergic system as changes in filter gain, and studied its effects related to the global-, local-, and meso-scale features of the connectome. At the global-scale, we found that the ability of the network of transiting through a variety of dynamical states is disrupted by randomization of the connection weights. By simulating neuromodulation of partial subsets of nodes, we found that transitions between integrated and segregated states are more easily achieved when targeting nodes with greater connection strengths—local feature—or belonging to the rich club—meso-scale feature. Overall, our findings clarify how the network spatial features, at different levels, interact with neuromodulation to facilitate the switching between segregated and integrated brain states and to sustain a richer brain dynamics.


INTRODUCTION
The human brain generates a rich repertoire of spatiotemporal dynamics characterized by the integrated and segregated functional states (Tononi, 2004). Information processed in parallel by domain-specific systems (segregated) is brought together (integrated) to guide adaptive behavior (Dehaene and Changeux, 2011). The balance between segregation and integration is essential to coordinate the local and global communication of neural information, it is needed to support a wide variety of cognitive functions, and has been proposed as a prominent organizational principle of the brain (Sporns, 2013;Cohen and D'Esposito, 2016;Shine, 2019;Wang et al., 2021). The dynamics and flexibility of brain activity, necessary for the coherent global functioning of the brain, enables the coexistence of segregated and integrated brain states (Kelso, 2012;Tognoli and Kelso, 2014;Wang et al., 2021).
Neuroimaging recording techniques such as electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) allow the characterization of functional connectivity (FC) of the brain, from which the functional integration and segregation can be quantified using network theory tools (Bullmore and Sporns, 2009;González et al., 2016). The observed patterns of FC reflect the diversity of neuronal dynamics that emerge, among others, from the nonlinear dynamics of brain regions interconnected through structural connectivity (SC) (Deco and Jirsa, 2012;Lord et al., 2017;Guan et al., 2020). FC continuously evolves even in resting conditions (Allen et al., 2014;Hansen et al., 2015;Cabral et al., 2017), moreover, it changes across several tasks, highlighting the flexible nature of brain dynamics (Cohen and D'Esposito, 2016;Shine et al., 2016Shine et al., , 2019Wang et al., 2021).
A plausible mechanism to facilitate-and regulatetransitions between different FC patterns are neuromodulatory systems. Neuromodulators do not directly excite neurons. Instead, they change their excitability and response to neurotransmitters, increasing or decreasing the probability of firing action potentials (Thiele and Bellgrove, 2018). The role of the cholinergic and noradrenergic systems in managing the segregation/integration balance has been evidenced in experimental (Shine et al., 2016(Shine et al., , 2018bPfeffer et al., 2020), and theoretical frameworks (Shine et al., 2018a;Pfeffer et al., 2020;Coronel-Oliveros et al., 2021).
The noradrenergic system is involved in arousal when subjects engage in high-load cognitive tasks (Aston-Jones and Cohen, 2005;Shine et al., 2016Shine et al., , 2018b. For example, in fMRI recordings during an N-back task (for assessing working memory), the pupil diameter-a marker of noradrenergic tone (Reimer et al., 2016)-increases (Shine et al., 2016(Shine et al., , 2018b. The principal source of noradrenaline in the cerebral cortex comes from the locus coeruleus (LC) (Fuxe et al., 2010). The GANE model of gain modulation (Mather et al., 2016;Lee et al., 2018), proposes that the noradrenergic system modulates neural response through an excitatory feedback loop between glutamate receptors on varicosities of LC projections and adrenergic β receptors on presynaptic glutamatergic neurons. At the same time, less activated neurons are suppressed through the action of adrenergic α 2 autoreceptors expressed on the varicosities. The overall result comprises an increase of the neuron responsivity above a threshold, and a decrease of the responsivity below this threshold. This is equivalent to increasing the slope of the input-output sigmoid function, also named filter gain, as proposed in Servan-Schreiber et al. (1990) and Aston-Jones and Cohen (2005).
In a recent article (Shine, 2019), the noradrenergic system was considered to promote an integrated functional network configuration increasing the filter gain (Servan-Schreiber et al., 1990;Aston-Jones and Cohen, 2005;Mather et al., 2016;Thiele and Bellgrove, 2018). Computational studies (Shine et al., 2018a;Coronel-Oliveros et al., 2021) have also shown how the interplay between cholinergic and noradrenergic systems can regulate the segregation/integration balance. While recent theoretical articles point out that a non-uniform neuromodulation can explain better the effects of neuromodulatory systems in brain dynamics (Deco et al., 2018;Kringelbach et al., 2020), most studies so far have considered homogeneous neuromodulation, i.e., acting in all nodes in the same way.
There is evidence about the importance of network properties of the human connectome (Cabral et al., 2014;Zamora-López et al., 2016;Wang et al., 2019;Castro et al., 2020). For example, its hierarchical modular organization is needed to sustain a richer brain dynamics (Zamora-López et al., 2016;Wang et al., 2019). Then, the repertoire of network configurations, as a way to conceptualize the dynamical richness, can be affected by neuromodulation. Using a neural mass model to simulate neural activity, Shine et al. (2018a) showed that rich club regions were strongly neuromodulated compared with non-rich club members, especially between the transition from functional segregation to integration. This work notably suggests that some particular brain regions play a key role in the switching between different functional states via neuromodulation. Here, instead of quantifying what regions would be strongly neuromodulated, we studied how much the impact would be on integration and segregation when neuromodulating specific subsets of nodes, and analyzed the structural features that define the nodes that, upon modulation, have the largest effect on the network dynamics as a whole.
To investigate this issue, we built a whole-brain model based on the Jansen and Rit equations (Jansen et al., 1993;Jansen and Rit, 1995) coupled to a human SC matrix, that allows us to simulate the effect of the noradrenergic system on the functional integration and segregation features of the network (Coronel-Oliveros et al., 2021). The interaction between neuromodulation and structural connectivity was studied at three levels: at the global-scale, we used random surrogate connectomes that preserve the number and strength of connections but disrupt the global patterns. At the meso-scale, we determined whether the modulation of a node subset containing the anatomical rich club (Opsahl et al., 2008;Van Den Heuvel and Sporns, 2011) or the critical s-core (Garas et al., 2012;Eidsaa and Almaas, 2013), is optimal to produce a change in network dynamics, compared to randomly chosen subsets. At the localscale, we explored which local properties define the set of nodes that, when being neuromodulated, maximize the effect on network dynamics.
We found that when we selectively neuromodulated the brain regions by the rich club (meso-scale property) or the high strength criteria (local-scale) the whole-brain network dynamics is most effectively modified. Additionally, we observed that surrogate connectomes reduced FC richness, compared with human SC, when neuromodulated. Overall, our findings clarify how the neuromodulation interacts with the anatomical network features at local-, meso-, and macro-scale levels in a whole-brain FIGURE 1 | Whole-brain neural mass model with neuromodulation. (A) The Jansen and Rit model is composed of a population of pyramidal neurons with excitatory and inhibitory feedback mediated by interneurons. A series of constants C i connect each population. The outputs are transformed from average pulse density to average postsynaptic membrane potential by an excitatory (inhibitory) impulse response function h E (t) [h I (t)]. Then, a sigmoid function S performs the inverse operation. Pyramidal neurons project to distant brain areas and receive both uncorrelated Gaussian-distributed inputs p(t) and inputs from other regions z(t), scaled by a global coupling parameter α. (B) Each node represents a cortical region, whose dynamics are ruled by the Jansen and Rit equations. The structural connectivity matrix, M, is the map of the connections (and their weights in the color bar) between cortical regions (row and columns of the matrix). The noradrenergic system increments pyramidal neuron responsivity to relevant stimuli with respect to noise, as a filter, by increasing the slope r 0 of their sigmoid function. (C) The whole-brain model comprises 90 cortical and subcortical regions linked by a human connectome. For each region, the model produces both EEG-like and BOLD-like signals. The brain figure was obtained using the BrainNet Viewer (Xia et al., 2013). model to facilitate switching between segregated and integrated brain states.

RESULTS
To study the effect of neuromodulatory systems on the integrative/segregative capacities of the human connectome, we used a whole-brain model of brain activity (Coronel-Oliveros et al., 2021). In this model, each node corresponds to a brain area represented by a neural mass, which consists of three populations (Jansen et al., 1993;Jansen and Rit, 1995): pyramidal neurons, excitatory, and inhibitory interneurons ( Figure 1A). We used the same parameters as in Jansen et al. (1993) and Jansen and Rit (1995), except the connectivity constant from inhibitory interneurons to pyramidal neurons C 4 , which we modified to C 4 = 0.5C, being C the original intra-area connectivity constant of the model. The nodes are connected through a weighted undirected structural connectivity matrix derived from human MRI data (Deco et al., 2018), parcelated in 90 cortical and sub-cortical regions with the automated anatomical labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002; Figure 1B). Pyramidal neurons connect regions (or nodes) because it is considered that long-range projections are mainly excitatory (Gilbert et al., 1990;McGuire et al., 1991). The simulations generate firing rates in each node of the network, which was used as an input to a generalized hemodynamic model (Stephan et al., 2007). In this way, we obtained fMRI BOLD-like signals ( Figure 1C) from which we built the FC matrices.
We modeled the influence of the noradrenergic system through the manipulation of the filter gain (Aston-Jones and Cohen, 2005;Shine, 2019; Figure 1B). The filter gain r 0 modifies the sigmoid function slope of pyramidal neurons, increasing their responsivity to relevant stimuli, decreasing the response to low amplitude stimuli, and boosting the signal-to-noise ratio.

Human Structural Connectivity Enhances Dynamical Richness
First, we analyzed how neuromodulation depends on the connectivity pattern of the human connectome by using different randomized surrogate connectomes. We employed a degree-and strength-preserving randomization (DSPR), which randomizes the structural connectivity while preserving original degree and strength distributions ( Figure 2B); in this way we can study the effect of disrupting the global connectivity without altering the local nodal properties. In addition, we employed a complete randomization of the structural connectivity ( Figure 2C), which does not preserve the degree and strength distributions. Finally, a homogenization (binarization) of the connectome was considered ( Figure 2D); this surrogate preserves the topology, disrupting the non-uniform weight distribution. We simulated EEG-like and fMRI BOLD-like signals from the Jansen and Rit model at different combinations of α ∈ [0, 1] and r 0 ∈ [0, 1] FIGURE 2 | Effects of network structure in neural synchronization and integration. From top to bottom: structural connectivity matrix, phase synchronyR, global efficiency E w (a measure of integration), and modularity Q w (measure of segregation), obtained in the model with different structural connectivities.R was obtained from EEG-like simulated activity, while E w and Q w were calculated using the FC obtained from the corresponding fMRI BOLD-like traces. parameters. Here, the value of the parameters is equal for all the nodes, and we refer to this case as uniform neuromodulation. We computed the mean of the Kuramoto order parameter, also known as phase synchronyR, the global efficiency, E w , and the modularity, Q w , as measures of global phase synchronization, integration, and segregation, respectively. Global efficiency is a measure of integration defined as the inverse of the characteristic path length (Rubinov and Sporns, 2010). High values of E w represent an efficient coordination between all pairs of nodes in the network, a signature of integration. Modularity is a measure of segregation based on the detection of network communities, or modules (Rubinov and Sporns, 2010); high modularity Q w is associated with segregation and vice-versa. Figure 2A shows how neuromodulation of the human connectome causes a shift of the model toward a synchronized and integrated state, with maximum integration observed in an intermediate region of the parameter space, as previously reported in Shine et al. (2018a) and Coronel-Oliveros et al. (2021). The synchronyR has an upper bound of 0.76, that is, the network never fully synchronizes. The transition is gradual, with many regions showing an intermediate behavior characterized by higher metastability and richer dynamics (Zamora-López et al., 2016;Shine et al., 2018a;Coronel-Oliveros et al., 2021). Moreover, the region of the parameter space whereR increases matches the increment in global efficiency, E w , verifying a link between the fast dynamics of EEG and the slower one of fMRI-BOLD.
We repeated the same exploration using the DSPR surrogate connectome ( Figure 2B; Rubinov and Sporns, 2011). The area of synchronized activity in the parameter space (r 0 , α) is reduced, and a spot of over-synchronized activity can be appreciated. Most importantly, the area of intermediate values of synchrony and integration is largely reduced, suggesting a reduction of dynamical richness. When the connectivity matrix is completely randomized (Figure 2C), or made homogeneous by assigning equal weights to all connections ( Figure 2D), neuromodulation produces a large area of over-synchronized activity in the parameter space and fewer regions with intermediate behavior.
The dramatic decrease in E w in Figures 2C,D is a consequence of the over-synchronization (R ≈ 1) triggered by randomization. When signals are highly synchronized in our model, the envelope in the alpha band of the EEG (between 8 and 13 Hz) becomes flat, and so does the BOLD-like signal calculated with the hemodynamic model (Foster et al., 2016). For this reason, this drop in E w should not be interpreted as a reduction of integration but a limitation of the hemodynamic model we employed in input simulations. Nevertheless, an over-synchronized regime of activity is a feature never found in the healthy brain (Miron-Shahar et al., 2019).
Thus, in line with several previous reports (Cabral et al., 2014;Zamora-López et al., 2016;Wang et al., 2019;Fukushima and Sporns, 2020), disrupting the organization of the human connectome (or the weight relationships between nodes) causes over-synchronization, and highly metastable regimes can not be easily reached employing neuromodulation.
In the following, we will study which local-or meso-scale organization features are determinant in the effect of neuromodulation of human connectome by evaluating the network behavior when changing the r 0 parameter in subsets of network nodes.

Neuromodulation of High-Strength Nodes Promotes Better Functional Integration
In this section, we investigate the impact on functional integration when an increasing number of nodes are neuromodulated. The order in which nodes are modulated is defined considering nodal measures obtained from the structural matrix M. We calculated, for each node i ∈ [1 . . . n]: node strength, K w i , nodal efficiency, E w i , and clustering coefficient, C w i (Rubinov and Sporns, 2010). The superscript w indicates the use of the weighted versions of the measures. Then, for each metric, nodes were ordered either from high to low or from low to high. We fixed the global coupling α = 0.65, and swept r 0 ∈ [0.33, 1] and the number of nodes being neuromodulated in [0, 90] in steps of three. As before, we used the EEG-like and BOLD-like signals to extract synchrony, integration, and segregation. A particular example of partial neuromodulation is shown with some detail in Figure 3. The (α, r 0 ) parameter space is shown in Figure 3A depicting global phase synchronyR, global efficiency E w and modularity Q w in a uniform neuromodulation scenario (all nodes identical). Figure 3B shows sample BOLD traces, the functional connectivity (FC) and the functional connectivity dynamics (FCD) matrices for α = 0.65, r 0 = 0.33 (red dot in Figure 3A). The FCD matrix visually represents the dynamical richness of the network, by computing timedependent FCs using sliding windows (Cabral et al., 2017;Orio et al., 2018). Then, FCs are vectorized and compared to each other using the Clarkson distance (Clarkson, 1936), resulting in a matrix of time vs. time. At the bottom, Figure 3D shows the same analysis for α = 0.65, r 0 = 0.67 (blue dot in Figure 3A). In the middle, Figure 3C shows the results when only half of the nodes have been neuromodulated to r 0 = 0.67 while the rest remain with r 0 = 0.33. As the number of nodes with r 0 = 0.67 increases, the FC matrices become more integrated (high E w Blue curves for neuromodulation of nodes with low strength, orange the opposite, and green for a random ordering of the nodes. Shaded areas correspond to 95% confidence intervals, for 10 realizations. and low Q w values). Similarly, the FCD matrices change from incoherence (red FCD), to exhibit multi-stable behavior (FCD with yellow-green patches), and finally to show correlated FC patterns (blue FCD). In summary, the increment of the number of neuromodulated nodes increases phase synchrony, functional integration, and the time correlation of FCs captured by the FCD. Figure 4 shows the result of neuromodulating r 0 with a target value in the [0.33, 1] interval and with the number of neuromodulated nodes ranging from 0 to 90. The order in which nodes are neuromodulated is either from low to high K w i ( Figure 4A) or viceversa ( Figure 4B). When the number of neuromodulated nodes is large,R and E w raise markedly in both cases; the opposite can be observed for Q w . However, picking the nodes of high strength first ( Figure 4B) has greater impact in the change of those metrics. The difference is best appreciated in Figure 4C, where we selected a target r 0 value of 0.67. There, the curves for the high to low K w i sorting (in orange) present a larger effect at the beginning, compared with the low to high K w i sorting (in blue). We can conclude that nodes with higher strength have a greater impact on functional integration, and inspection of the colormaps of Figures 4A,B reveals that this is true for almost all values of target r 0 . The results were also compared with a random selection of nodes for neuromodulation (green curves). As the blue curve is mainly below the green curve, neuromodulation of nodes with low K w i produces less synchronized and integrated dynamics than expected by a random neuromodulation. In contrast, choosing high K w i nodes is not different from random selection, when looking at the measures of integration and segregation. In consequence, there is a range (or possibly a set) of nodes that produce a robust integration when neuromodulated, compared to a random choice of nodes.
We compared the results of sorting the nodes based on strength K w i , with ranking the nodes based on nodal efficiency E w i or clustering coefficient C w i (Figure 5). A node with a high E w i has many short paths to the rest of the nodes of the network, while a high C w i is expected for nodes whose neighbors are FIGURE 5 | Incremental neuromodulation based on nodal efficiency and clustering coefficient. (A,B) Phase synchronyR, global efficiency E w , and modularity Q w as a function of the number of neuromodulated nodes, for r 0 = 0.67 as target value. In (A), nodes were sorted according to nodal efficiency, E w i , and in (B) according to the clustering coefficient, C w i . Blue curves for sorting of nodes from low to high (of the particular metric), orange is from high to low, and green for a random sorting of the nodes. (C) Difference between the area under the curve (AUC) of high vs. low sorting, averaged over the 10 realizations. Shaded areas in (A,B) correspond to 95% confidence intervals, and barplots were built using the mean ± standard deviation. ***p < 0.001. also connected between them. Figure 5A shows the result of modulating an increasing number of nodes from a basal r 0 = 0.33 to a target r 0 = 0.67, when the nodes are ordered from low to high or high to low E w i . The results are similar to the ones obtained using the strength K w i : neuromodulation of nodes of high E w i has a greater impact in synchronization and integration, compared with the nodes of low E w i . Here, the random sorting of nodes is similar to the high to low E w i . However, when the nodes are sorted according to their clustering coefficient C w i (Figure 5B), there is little difference compared to random sorting.
When comparing the results in Figures 5A,B with the neuromodulation of a random subset of nodes (green curves), there is no clear advantage of selecting the nodes with high E w i or C w i . Despite the increase inR being slightly higher for the orange curves, compared with the green curves, the difference in E w is unnoticeable, except in Q w when ordering the nodes from high to low E w i . These results contrast with the ones in Figure 4C, where the neuromodulation of nodes with high K w i produced an increase in synchronization and integration higher than random neuromodulation.
To summarize these results, we computed the difference between the area under the curve (AUC) for the high-to-low minus low-to-high (orange minus blue AUCs; Figure 5C). A larger difference implies a higher impact of neuromodulating first the nodes with a higher value of the chosen metric in synchronization, integration, and segregation. The mean difference in the AUCs for the measures R , E w , and − Q w (note that the sign is inverted for visualization purposes), is lower for C w i than for K w i and E w i (p < 0.001 for all comparisons using Student's t-test).
FIGURE 6 | Neuromodulation based on the rich club organization. (A) Normalized and weighted rich club coefficient φ w norm (K) (blue curve), as a function of the degree-based threshold K. This coefficient was calculated as the ratio between the human rich club coefficient (purple solid curve) and the mean coefficient for random surrogates (DSPR, purple dashed curve). The red arrow marks the point in which φ w norm (K) is maximal; rich club nodes were found at that point. (B) Schematic depiction of the rich club, feeders (not belonging to the rich club, but connected with it) and local nodes, (connected only to feeders). At the right, we show a glass brain with the nodes identified as rich club (n = 17 nodes), feeders (n = 60), and local (n = 13). (C) Changes in synchronyR, global efficiency E w and modularity Q w when neuromodulating 24-node sets containing the rich club (blue), local nodes (green), or only feeders (orange). The results are shown as the difference with respect to a random subset of nodes of equal size (null case). The bottom row summarizes the area under the curve (AUC) for each metric and nodal category, averaged over the 10 realizations. Shaded areas correspond to 95% confidence intervals, and bar plots were built using the mean ± standard deviation. **p < 0.01, ***p < 0.001.

Neuromodulation of Rich Club Nodes Strongly Impacts Functional Integration
Node strength, nodal efficiency and clustering coefficient are considered local-scale properties, i.e., they belong to each node. Several meso-scale network properties have been described as being determinant for network dynamic too, such as the rich club organization (Van Den Heuvel and Sporns, 2011) and the s-core (Hagmann et al., 2008;Garas et al., 2012;Eidsaa and Almaas, 2013;Castro et al., 2020). We identified the nodes belonging to the "rich club, " using the weighted rich club coefficient φ w (K), where K is a threshold based on degree (Opsahl et al., 2008). The rich club comprises a subset of the graph, thresholded at K, in which nodes are more strongly interconnected than expected by chance (Van Den Heuvel and Sporns, 2011). The coefficient is normalized using random surrogates φ w rand (K), in our case DSPR surrogates (Rubinov and Sporns, 2011). If the normalized coefficient φ w norm (K) is greater than 1, the network has a rich club organization at threshold K. Figure 6A shows a plot of φ w norm (K) (blue) as a function of K. The red arrow marks the point in which the normalized coefficient is maximal [φ w norm (K) = 1.367, p < 0.002]. Then, we identified feeder nodes-nodes that do not belong to the rich club but are connected to at least one of its members-and local nodesconnected to feeders but not to the rich club ( Figure 6B). We found 17 nodes belonging to the rich club, 60 feeders and 13 local nodes ( Figure 6B). The rich club members are the brain regions in Table 1.
As the analysis of the rich-club properties of the human SC defines sub-networks, instead of sorting the nodes, we chose a different approach than the neuromodulation of increasing subsets of nodes. Here, we simulated neuromodulation of a fixed-size subset of nodes, that included all nodes belonging to a certain category (rich club, feeders, or local). Because the categories differ in size, we complemented the rich club and local nodes with 7 and 11 nodes, respectively, selected randomly from the feeders. For the last one, we randomly chose 24 feeder nodes. Also, we had a null case, composed of 24 nodes randomly selected from the complete set of nodes. We repeated the random selection of nodes with 10 realizations, always using subsets of 24 nodes. The nodes started with a basal r 0 value of 0.33, and r 0 was swept up to 1 but only in the designated subset of nodes. For each r 0 increment, we measured R, E w , and Q w . Then, we subtracted to each measurement the result of the null case. The results are shown in Figure 6C. Neuromodulation of the rich club nodes produces an increase in synchronization and integration, and a decrease in modularity, above chance. The difference becomes more pronounced with further increments of r 0 . Opposite results were observed for the subsets containing local nodes. Finally, neuromodulation of subsets containing only feeder nodes produce no difference compared to random selection of nodes. As a summary index, we calculated the AUC for each nodal category ( Figure 6C). Considering the three measurements, the AUC is higher (lower in the case of modularity Q w ) for the rich club respect to feeders and local, and higher (lower in the case of Q w ) between feeders and local (p < 0.001 for all comparisons using Student's ttest). Our results show that noradrenergic neuromodulation of a subset including the rich club nodes has a greater impact on integration compared to the feeders, locals, and a random selection of nodes.
As previously shown, functional integration is also achieved by neuromodulation of highest strength nodes. To highlight the difference between the local and meso-scale approaches, we quantified the overlap between the rich club nodes and the 17 nodes with higher strength. We found that only 8 members of the rich club belong to the subset of 17 nodes with higher strength (Table 1). Thus, there are some high-strength key nodes that do not belong to the rich club, that promote functional integration via neuromodulation.
To explore a second meso-scale network organization, we performed a s-core decomposition (Garas et al., 2012;Eidsaa and Almaas, 2013) that classifies nodes according to their coreperiphery organization (Hagmann et al., 2008; Figure 7A). We defined three categories considering a range of critical s-core values: S 3 with 10 nodes (1.54 < s < 1.78), S 2 with 56 nodes (1.48 < s < 1.54), and S 1 with 24 nodes (s < 1.48; Figure 7B). The critical s-core is defined as the maximal value of s at which nodes are still connected to the network. Thus, S 3 are nodes connected within them with highest strength, S 2 middle-strength nodes, and S 1 the nodes with the lowest strength. The S 3 subset comprises the brain regions shown in Table 1.  We simulated the neuromodulation in subsets of 24 nodes, containing either the S 3 or the S 1 category, and complementing S 3 with 14 random nodes from S 2 as done with the rich club. A third group was built with 24 nodes randomly selected from S 2 , and all groups were compared to a random selection of 24 nodes from the whole set. As shown in Figure 7C, the selection of S 2 nodes for neuromodulation shows the largest effect in R , E w , and Q w , compared with S 1 nodes (p < 0.001 using Student's ttest) and compared to the selection of S 3 nodes (p < 0.001, except for E w with p = 0.106). Thus, nodes belonging to the highest s-core (nodes of the highest within-strength sub-network) do not behave like the rich club, as their neuromodulation does not have the highest impact on network synchronization and integration/segregation properties.

DISCUSSION
In this work, we sought to identify the relationship between structural features of the human connectome and the specific set of regions that, when neuromodulated in a biologically realistic whole-brain model, produce a significant increase in functional integration. We found that the global organization of the connectome sustains rich metastable and partially synchronized states, essential to the effects related to neuromodulation. At the meso-and local-scales, nodes belonging to the anatomical rich club, and those having high nodal strength, produce a marked increase in functional integration (and a decrease in segregation) when neuromodulated.
Our results show that the whole-brain model exhibits oversynchronized behavior when using surrogate connectomes, restricting the dynamic features of the model. This result is in the FIGURE 7 | Neuromodulation based on the s-core decomposition. (A) s-core decomposition. At the left, example based on degree (k-core). Nodes are recursively removed based on a degree threshold. The remaining nodes form a subgraph or core where all nodes have a within-degree above the threshold. For example, the three-core of the figure corresponds to a subgraph where all nodes have a degree of three or more. The numbers on the circles correspond to nodes' degree. The right plot shows the number of remaining nodes in s-core after applying the strength-based threshold s. (B) Brain regions identified using the s-core decomposition. We defined three categories, considering a range of s values: S 3 (n = 10 nodes), S 2 (n = 56), and S 1 (n = 24). (C) Changes in synchronyR, global efficiency E w and modularity Q w when neuromodulating 24-node sets containing the S 3 nodes (blue), S 1 nodes (olive green), or only S 2 nodes (pink). S 1 and S 3 sets were complemented with random nodes from S 2 to obtain sets of 24. Results are shown as the difference with respect to a random subset of nodes of equal size (null case). The bottom row summarizes the area under the curve (AUC) for each metric and nodal category, averaged over the 10 random seeds. Shaded areas correspond to 95% confidence intervals, and bar plots were built using the mean ± standard deviation. ***p < 0.001. same line as other previous findings (Cabral et al., 2014;Zamora-López et al., 2016;Wang et al., 2019;Fukushima and Sporns, 2020). Here, we show this behavior in the (α, r 0 ) parameter space, where simulations with randomized connectomes show either incoherent or over-synchronized activity. Using a wholebrain model to simulate and fit magnetoencephalography (MEG) resting-state recordings, Cabral et al. (2014) found not only that randomized and homogenized versions of the human structural connectivity did not fit empirical data; moreover, they found that the fit was maximal in the metastable region of the parameter space, when unsynchronized (segregated) and synchronized (integrated) regimes of activity coexist. In the same way, Fukushima and Sporns (2020) using more complex surrogate data in the context of whole-brain models, found features of the human connectome that better capture the dynamic fluctuations in fMRI resting-state recordings. Additionally, computational studies conducted by Zamora-López et al. (2016) showed that the human connectome better maximizes functional complexity in fMRI recordings, compared with different surrogate connectomes. Finally, Wang et al. (2019) analyzed how the hierarchical modular structure of the human connectome enables the coexistence of segregated and integrated functional states, also with the use of network surrogates in which hierarchical levels were controlled. Our study interpret the explorations of the parameter space as levels of neuromodulation, that allow the brain to tune its integration or segregation levels to environmental demands. However, neuromodulation cannot bring back a dynamically rich regime to a network without a structural connectome that sustains it.
At the local level, the effects of neuromodulation strongly depend on the characteristics of the nodes in the human connectome. In our model, the nodes with high strength are the ones that better facilitate functional integration when neuromodulated. This result resonates with a recent work by Herzog et al. (2020), who studied a whole-brain model fitted to reproduce the effects of lysergic acid diethylamide (LSD) in resting-state brain dynamics. In their model, the serotonergicinduced changes in nodal entropy correlated positively with node strength. Notably, the correlation disappears when the human connectome was randomized without preserving the strength distribution, emphasizing the importance of the specific organization of the human connectome in shaping brain dynamics. Interestingly, the entropy changes described by Herzog et al. (2020) are poorly explained by the 5HT 2A receptor density map, obtained by PET (Beliveau et al., 2017), and depends on both node strength and receptor density. Thus, the interaction between the structural connectivity, receptor density, and neuromodulation is not straightforward. A similar complex picture arises when our results are contrasted with receptor maps of noradrenergic receptors (see below).
Network hubs, or nodes belonging to the rich club or network's ignition core, can be critical elements for binding information of segregated brain regions, that is, to integrate information across brain areas (Griffa and Van den Heuvel, 2018;Castro et al., 2020). Considering the relevance of integration for the brain function (Tononi, 2004), and the noradrenergic influence on integration (Shine, 2019), we hypothesized that anatomical network hubs are pivotal elements for promoting functional network integration. Our results confirmed this hypothesis, being the neuromodulation of rich club nodes the one that most effectively facilitates functional integration and synchronization of brain activity. This result agrees with findings reported in a fMRI resting-state model of the brain by Deco et al. (2017), where removing the rich club nodes causes a larger decrease in integration compared to the removal of the non-rich club members. Similar results have been found in computational models of noradrenergic neuromodulation where rich club nodes are strongly neuromodulated causing functional networks to switch from segregation to integration (Shine et al., 2018a).
Notably, neuromodulation of nodes belonging to the critical s-core (the maximally inter-connected core) does not promote integration as the rich club nodes do. Both meso-scale analyses rely on sets of nodes organized with strong connection weights. However, they do it differently. The rich club coefficient threshold is based on degree, and rich club members are highly connected between them as well as with the non-rich club members. In contrast, the s-core decomposition find subsets of nodes highly interconnected at strength s, but not necessarily well connected to the rest of the network. Thus, the whole-network changes are more easily achieved if the set of nodes to be neuromodulated is highly connected both between them and with the rest of the network. The rich-club organization captures additional information that is missing in the local (weight) analysis. For example, the 17 rich club nodes have an overlap of ≈50% with the 17 highest strength nodes. In contrast, nodes belonging to the S 3 category are the nodes of the highest strength in the network; however, they cannot boost functional integration to the same extent as the rich club nodes.
Part of the brain regions we found in the rich club support high order brain functions. For example, frontoparietal regions play an important role in cognition, and are markedly activated when subjects engage in cognitive tasks (Cavanna, 2007). Precuneus has been associated with consciousness, and a decrease in its activity was reported in sleep, anesthesia, and vegetative states (Lückmann et al., 2014). Thalamus, the brain "relay station, " strongly connects several networks that comprise multiple cortical regions (Hwang et al., 2017). Multi-task fMRI recordings in humans suggest a robust role of the anatomical rich club as facilitating elements of functional integration in overall tasks . An extended analysis and discussion about the role of the rich club, in both health and disease, can be found in Griffa and Van den Heuvel (2018).
The non-uniform expression of receptors across several brain areas suggests that the brain uses selective or partial neuromodulation. In this way, the effect of the noradrenergic system on filter gain may be modeled as proportional to adrenergic receptor expression. Experimentally, the optogenetic activation of the LC in mice increased average functional connectivity, which correlates with the expression of α 2 , α 1 , and β 1 adrenergic receptors (Zerbi et al., 2019). Thus, a future research avenue in computational models may include a densitydependent noradrenergic neuromodulation with the addition of some receptors maps, obtained by positron emission tomography (PET), or even gene expression maps (Shen et al., 2012) that correlate with receptor density maps (Komorowski et al., 2017). Surprisingly, using the Allen Human Brain Atlas database (Shen et al., 2012) we found that some adrenergic receptor genes, i.e., the ADRA2A and ADRB1, are less expressed in the rich club nodes than in feeders and local nodes (Figure 8). As a consequence, the noradrenergic-mediated increase in filter gain could have a lower impact on rich club nodes. It is possible that this reduced expression constitutes a compensation for the high connectivity of rich club nodes, specially considering the higher metabolism of rich regions that exposed them to oxidative stress and neuroinflammation (Griffa and Van den Heuvel, 2018). On the other hand, receptor expression can itself be compensated by a specific sub-cellular localization or other excitability factors that may enhance the effect of noradrenaline.
It has been suggested that the effect of noradrenaline in functional connectivity is context-dependent (Shine et al., 2018b;Pfeffer et al., 2020). In that line, modeling the effect of noradrenaline in resting-state and task conditions could untangle the mechanisms behind this context-dependent effect of noradrenaline. The anatomical backbone and other dynamical parameters of this model can be substituted to study the mouse or monkey brain and to any other species for which the whole-brain white-matter connectivity is available.
Our work considers an arbitrary basal value of r 0 . Despite this, we reported a clear effect of the selective noradrenergic neuromodulation on functional integration, that is, some brain FIGURE 8 | Expression of some noradrenergic receptors genes in brain regions. Genes ADRA2A, ADRA2C, and ADRAB1 are related to noradrenergic receptors α 2A , α 2C , and β 1 , respectively. The normalized expression was obtained from the Allen Human Brain Atlas using the AAL parcelation. Bar plots were built using the mean ± standard deviation. ***p < 0.001, **p < 0.01, *p < 0.05, ∼p < 0.1.
regions have a greater impact in the noradrenaline-mediated effect on brain function. A further improvement to our approach constitutes the use of a different benchmark, e.g., fitting the model to reproduce the empirical FC in resting-state, and then apply a homogenoeus or selective neuromodulation. Furthermore, the addition of receptors maps may be considered, as commented above.
Overall, our results offer new insights into the key regions of the human brain that, when neuromodulated via the noradrenergic system, promote transitions to integrated functional states. Our results highlight the importance of the rich club and high-strength connections in producing changes related to neuromodulation. We hope that our theoretical framework inspires new research toward clinical applications or treatments of human brain disorders caused by or associated with changes in functional and structural brain connectivity.

Whole-Brain Neural Mass Model
We simulated neuronal activity using the Jansen and Rit neural mass model (Jansen et al., 1993;Jansen and Rit, 1995). In this model, a brain area consists of three populations of neurons: pyramidal neurons, excitatory and inhibitory interneurons. The dynamical evolution of the three populations within the brain area is modeled by two blocks each. The first block transforms the average pulse density in average postsynaptic membrane potential (which can be either excitatory or inhibitory; Figure 1A). This block, denominated post synaptic potential (PSP), is represented by an impulse response function. For the excitatory outputs: and for the inhibitory ones The constants A and B define the maximum amplitude of the PSPs for the excitatory (EPSPs) and inhibitory (IPSPs) cases respectively, while a and b represent the inverse time constants for the excitatory and inhibitory postsynaptic action potentials, respectively. The second block transforms the postsynaptic membrane potential in average pulse density, and is given by a sigmoid function of the form S(ν, r) = ζ max 1 + e r(ν th −ν) , with ζ max as the maximum firing rate of the neuronal population, r the slope of the sigmoid function, ν th the half maximal response of the population, and ν their average PSP. Additionally, pyramidal neurons receive an external stimulus p(t), whose values are taken from a Gaussian distribution with mean µ = 2 impulses/s and standard deviation σ = 2. In this model ( Figure 1A), each node i ∈ [1 . . . n] represents a single brain area. The global coupling is scaled by a parameter α, and nodes are connected by a normalized structural connectivity matrix M ( Figure 1B). This matrix is derived from a human connectome (Deco et al., 2018) parcelated in n = 90 cortical and subcortical regions with the automated anatomical labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002); the matrix is undirected and takes values between 0 and 1. Because long-range connections are mainly excitatory (Gilbert et al., 1990;McGuire et al., 1991), only links between the pyramidal neurons of a node i with pyramidal neurons of a node j are considered. We applied a global normalization procedure to the structural connectivity matrix M. The normalization consisted of dividing all the values of the matrix by the mean strength of the nodes. The resulting normalized matrix M is defined as The set of equations, for a node i, includes the within and between nodes activitẏ where x 0 , x 1 , x 2 correspond to the outputs of the PSP blocks of the pyramidal neurons, and excitatory and inhibitory interneurons, respectively, and x 3 the long-range outputs of pyramidal neurons. The constants C 1 , C 2 , C 3 , and C 4 scale the connectivity between the neural populations (see Figure 1A). The first pair of equations, x 0 and y 0 , are related to the outputs of pyramidal cells to both interneurons; the second pair, x 1 and y 1 , represent all the local excitatory inputs that the pyramidal neurons receive; x 2 and y 2 constitute the inhibitory contribution to pyramidal cells. An additional pair of equations (x 3 and y 3 ) are introduced to represent long-range (inter-area) connections, as they target the apical dendrites of pyramidal neurons and thus their EPSP have a larger characteristic time constant. We used the original parameter values of Jansen and Rit (Jansen et al., 1993;Jansen and Rit, 1995), except for C 4 : ζ max = 5 s −1 , ν th = 6 mV, r 0 = r 1 = r 2 = 0.56 mV −1 , a = 100 s −1 , b = 50 s −1 , A = 3.25 mV, B = 22 mV, C 1 = C, C 2 = 0.8C, C 3 = 0.25C, C 4 = 0.5C, and C = 135. Changing C 4 from 0.25 C to 0.5 C allowed the model to sustain oscillations in a wider range of α values. The parameters A, B, a, and b were selected to produce IPSPs longer in amplitude and latency in comparison with the EPSPs. The inverse of the characteristic time constant for the long-range EPSPs was defined asā = 0.5a. This choice was based on the fact that longrange excitatory inputs of pyramidal neurons target their apical dendrites, and consequently this slows down the time course of the EPSPs at the soma (Branco and Häusser, 2011).
The input from brain areas j = i to the region i is given by The average PSP of pyramidal neurons in region i characterizes the EEG-like signal in the source space; it is computed as (Jansen et al., 1993;Jansen and Rit, 1995) The firing rates of pyramidal neurons ζ i (t) = S[ν(t) i , r 0 ] were used to simulate the fMRI-BOLD signals.

Neuromodulation
The effect of the noradrenergic system was simulated controlling the parameter r 0 (filter gain; Figure 1B), which represents the sigmoid function slope of the pyramidal population, and increases the signal-to-noise ratio of pyramidal cells (Servan-Schreiber et al., 1990;Thiele and Bellgrove, 2018). Details about the relationship between the noradrenergic system and filter gain can be found in the Introduction section. Further analysis about this relationship has been presented previously in Mather et al. (2016) and Shine (2019). We analyzed the effect of the noradrenergic neuromodulation in three scenarios: Macro-scale: Noradrenergic neuromodulation was studied in interaction with the cholinergic system, represented by the parameter α. The parameters were the same for all nodes. We changed the features of the connectivity matrix M (see Figure 2) to study the combined effect in neural coordination.
Meso-scale: nodes were classified in different categories, either according to the rich club organization (Van Den Heuvel and Sporns, 2011) or s-core decomposition of the network (see section 4.5; Garas et al., 2012;Eidsaa and Almaas, 2013). Global coupling was fixed in α = 0.65, and the basal value of r 0 was 0.33 for all nodes (Figure 3). We incremented r 0 in a subset of 24 nodes belonging to a particular category, and compared the results with the neuromodulation of a equal-length random subset of nodes.
Because the categories differ in the number of nodes, a fair comparison must considered subsets of equal size. To achieve that, we complemented the rich club with seven randomly selected feeder nodes, while the local nodes were complemented with 11 randomly selected feeders. Likewise, we complemented the S 3 category with 14 randomly selected S 2 nodes. From both the feeders and S 2 nodes we selected 24 nodes randomly. All subsets consisted on 24 nodes, were generated 10 times with different random seeds and the results averaged.
Local-scale: Nodes were sorted using one of three metrics: node strength K w i , nodal efficiency E w i , or clustering coefficient C w i (Rubinov and Sporns, 2010). We neuromodulatedincreasing r 0 -node by node in increments of three, considering the metric from high to low and vice-versa (Figure 3).

Simulations
Following Birn et al. (2013), we ran simulations to generate the equivalent of 11 min real-time recordings, discarding the first 60 s. The system of stochastic differential equations (5) was solved with the Euler-Maruyama method, using an integration step of 1 ms. We used 10 random seeds (realizations) which controlled the initial conditions and the stochasticity of the simulations. We simulated neuronal activity sweeping the parameters α ∈ [0, 1] and r 0 ∈ [0, 1], for the macroscale scenario. In the local-and meso-scale scenarios, we swept r 0 ∈ [0.33, 1] for a susbset of nodes, considering a basal value of r 0 = 0.33 and a fixed α = 0.65. All the simulations were implemented in Python and the codes are freely available at: https://github.com/vandal-uv/Structural_ Neuromod_2021.git. The graph analysis was performed using the Brain Connectivity Toolbox for Python (https://github.com/ fiuneuro/brainconn; Rubinov and Sporns, 2010).

Simulated fMRI-BOLD Signals
We used the firing rates ζ i (t) to simulate BOLD-like signals using a generalized hemodynamic model presented in Stephan et al. (2007). In this model, an increment in the firing rate ζ i (t) triggers a vasodilatory response s i , producing blood inflow f i , changes in the blood volume v i and deoxyhemoglobin content q i . The corresponding system of differential equations iṡ where τ s = 0.65, τ f = 0.41, τ v = 0.98, τ q = 0.98 represent the time constants for the signal decay, blood inflow, blood volume, and deoxyhemoglobin content, respectively. The stiffness constant (resistance of the veins to blood flow) is given by κ, and the resting-state oxygen extraction rate by E 0 . We used κ = 0.32 and E 0 = 0.4. The BOLD-like signal of node i, denoted B i (t), is a non-linear function of q i (t) and v i (t) where V 0 = 0.04 represents the fraction of venous blood (deoxygenated) in resting-state, and k 1 = 2.77, k 2 = 0.2, k 3 = 0.5 are kinetic constants.
The system of differential equations (8) was solved using the Euler method with an integration step of 1 ms. The signals were band-pass filtered between 0.01 and 0.1 Hz with a 3rd order Bessel filter. These BOLD-like signals were used to build the functional connectivity (FC) matrices from which the subsequent analysis of functional network properties was performed.

Macro-Scale
To compare different Macro-scale features of the connectome we used four connectivity matrices (see Figure 2). The first matrix corresponds to the original human connectome matrix (Human, Figure 2A) (Deco et al., 2018). The second to a degree and strength preserving randomization of the matrix (DSPR, Figure 2B; Rubinov and Sporns, 2011). The third to a randomization, which only preserves the weight distribution of the original matrix (Random, Figure 2C). The fourth matrix was built setting to 0 all entries of M ij < 0.05, and 1 otherwise (Homogeneous, Figure 2D).

Meso-Scale
We identified the nodes belonging to the "rich club" sub-network of the graph (Van Den Heuvel and . Nodes were ranked according to degree, and then a subgraph was built using a threshold K, retaining the nodes with a degree greater than K. For each K value the weighted rich-club coefficient was computed as (Opsahl et al., 2008).
where W >K is the sum of the weighted edges of the subgraph of nodes with a degree greater than K, E >K represent the total number of edges of the subgraph, and w rank l a vector that contains all the weighted edges of the entire network sorted from high to low values. If φ w (K) = 1, then the sum of the weights of the "rich nodes" is maximal. Otherwise, φ w (K) < 1 indicates the proportion of the weighted edges of network that are into the sub-network, and then some of the stronger connections were missed when applying the threshold K. The rich club coefficient was normalized in relation to DSPR surrogate graphs.
being φ w norm (K) the normalized rich club coefficient, and φ w rand (K) the mean rich club coefficient for a set of 1,000 random surrogates graphs. Values of φ w norm (K) > 1 indicates a rich-club organization, and nodes retained at K are defined as "rich club" nodes ( Figure 6A). The nodes that do not belong to the rich club, but are connected with these nodes are called "feeders." The remaining nodes correspond to "local" nodes ( Figure 6B). For a maximum φ w norm (K) = 1.367 (p < 0.002), we identified 17 "rich club" nodes, 60 feeder nodes and 13 local nodes ( Figure 6B). Because the high density of the structural matrix M (≈ 40%) hindered the discerning of the local nodes from feeders, we identify these nodes applying an absolute threshold of 0.05 to M. We selected this value as the maximum threshold that, when applied, preserves the fitting of the model to the empirical resting-state FC matrix.
The core-periphery organization (Hagmann et al., 2008) was analyzed performing a s-core decomposition (Garas et al., 2012;Eidsaa and Almaas, 2013), which identifies the cores of densely interconnected nodes in the network. The method consists in removing recursively a shell of nodes with strength less than s to obtain the network core nodes. The nodes were assigned to a category that corresponded to the maximal s value at which they are still connected to the network, defined as the critical s-core ( Figure 7A). We defined three categories for different s values: s 1 with 24 nodes (s < 1.48), s 2 with 56 nodes (1.48 < s < 1.54), and s 3 with 10 nodes (1.54 < s < 1.78; Figure 7B).

Local-Scale
We employed three different metrics to characterize individual nodes. Node strength (weighted degree) was computed as where N is the set of nodes and w ij the weighted edge of the matrix M (Rubinov and Sporns, 2010). We computed the nodal efficiency as where d w ij is the shortest path between the nodes i and j. Shortest paths were calculated from the sum of the inverse of the weights of M; the shortest path between two nodes (i, j) is the path that minimizes this sum (the distance). Using the shortest paths, nodal efficiency E w i was computed. Nodes with high values of E w i are those with high proportion of short paths to the rest of the nodes of the network (Rubinov and Sporns, 2010). Finally, we calculated the clustering coefficient for each node (Rubinov and Sporns, 2010) where t w i is the proportion of triangles around the node i, calculated as A node with a high C w i is highly connected with adjacent (local) nodes.

Phase Synchronization
As a measure of global synchronization, we calculated the Kuramoto order parameter R(t) (Acebrón et al., 2005) of the EEG-like signals ν(t) derived from the Jansen and Rit model. First, the raw signals were filtered with a 3rd order Bessel bandpass filter using their frequency of maximum power (usually between 4 and 10 Hz) ±3 Hz. Then, the instantaneous phase θ (t) was obtained using the Hilbert transform. The global phase synchrony is computed as: where θ i (t) is the phase of the oscillator i over time, j = √ −1 the imaginary unit, |•| denotes the module, N denotes the average over all nodes, and t the average over time.

Functional Integration and Segregation
Functional Connectivity (FC) matrices were built from Pearson correlations of the entire BOLD-like time series. Instead of employing an absolute or proportional thresholding, we thresholded the FC matrices using Fourier transform (FT) surrogate data (Lancaster et al., 2018) to avoid the problem of introducing spurious correlations (Fornito et al., 2013). The FT algorithm uses a phase randomization process to destroy pairwise correlations, preserving the spectral properties of the signals (the surrogates have the same power spectrum as the original data). We generated 500 surrogate time series of the original set of BOLD-like signals, to obtain the surrogate sFCs matrices. For each one of the (n 2 − n)/2 possible connectivity pairs (with n = 90) we fitted a normal distribution of the surrogate values.
Using these distributions, we tested the hypothesis that a pairwise correlation is higher than chance (that is, the value is at the right of the surrogate distribution).
To reject the null hypothesis, we selected a p-value equal to 0.05, and corrected for multiple comparisons with the FDR Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995) to decrease the probability of making type I errors (false positives). The entries of the sFC matrix associated with a p > 0.05 were set to 0. The result is a thresholded, undirected, and weighted (with only positive values) sFC matrix.
Integration was evaluated over the thresholded FC matrices. We employed the weighted version of the global efficiency (Latora and Marchiori, 2001). This measure of integration is based on paths over the graph: it is defined as the inverse of the average shortest path length. This metric is computed as being N the set of all nodes, n number of nodes, and d w ij the shortest path between the nodes i and j.
Segregation was quantified using modularity Q w , a metric for the detection of the network's communities (Rubinov and Sporns, 2010). The detection of so-called communities or network modules in the thresholded FC matrix, was based on the Louvain's algorithm (Newman, 2006;Blondel et al., 2008). We used the weighted version of the modularity (Newman, 2004) defined as where w ij is the weight of the link between the nodes i and j, l w is the total number of weighted links of the network, m i (m j ) the module of the node i (j), and k w i (m j ) the weighted degree (named also strength) of i (j). The algorithm assigns a module to each node in a way that maximizes the modularity (18). The Kronecker delta δ m i ,m j is equal to 1 when m i = m j (that is, when two nodes belongs to the same module), and 0 otherwise.
Because the Louvain's algorithm is stochastic, we employed the consensus-clustering algorithm (Lancichinetti and Fortunato, 2012). We ran the Louvain's algorithm 200 times with the resolution parameter set to 1.0 (this parameter controls the size of the detected modules; larger values of this parameter allows the detection of smaller modules). Then, we built an agreement matrix G, whose entries G ij ∈ [0, 1] indicates the proportion of partitions in which the pairs of nodes (i, j) share the same module. Then, we applied an absolute threshold of 0.5 to the matrix G, and ran the Louvain's algorithm again 200 times using G as input, producing a new consensus matrix G ′ . This last step was repeated until convergence to a unique partition.

Functional Connectivity Dynamics
The FCD matrix captures the evolution of FC patterns and, consequently, the dynamical richness of the network (Hansen et al., 2015;Cabral et al., 2017). We used the sliding window approach (Hansen et al., 2015;Orio et al., 2018), with windows of 100 s length and a displacement of 2 s between consecutive windows. The length was chosen on the basis of the lower limit of the band-pass filter (0.01 Hz), in order to minimize spurious correlations (Leonardi and Van De Ville, 2015). For each window, a FC matrix was calculated from the Pearson correlation of BOLD-like signals. We obtained 251 weighted and undirected FC matrices from the 600 s simulated BOLD-like signals. The upper triangular of each FC matrix is unfolded to make a vector, and the FCD is built by calculating the Clarkson distance λ(x, y) = 1 √ 2 x ||x|| − y ||y|| between each pair of FCs (Clarkson, 1936).

Gene Expression Maps
To quantify the expression of some noradrenergic receptor genes in brain regions, we used the microarray expression data of the Allen Human Brain Atlas (Shen et al., 2012). The dataset was processed and normalized employing the Abagen library for Python (https://github.com/rmarkello/abagen/tree/0. 1; Arnatkevicute et al., 2019), and then parcellated using the AAL atlas (Tzourio-Mazoyer et al., 2002). We compared the expression of the ADRA2A, ADRA2C, and ADRB1 genes in rich club, feeders and local nodes. Statistical comparison was performed with a Student's t-test for independent samples.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://github.com/vandaluv/Structural_Neuromod_2021.

AUTHOR CONTRIBUTIONS
CC-O, SC, RC, and PO contributed to conception and design of the study and wrote the manuscript. CC-O performed the simulations and statistical analysis. All authors contributed to manuscript revision, read, and approved the submitted version.