The ascending arousal system promotes optimal performance through mesoscale network integration in a visuospatial attentional task

Cognitive processes emerge from the dynamic interplay between diverse mesoscopic brain systems (Shine, 2021; Shine et al., 2016). Thus, the neural activity supporting cognition does not exist in a vacuum, but instead is deeply embedded within the ongoing dynamics of the physiological networks of the body (Varela et al., 2001). In particular, the neural processes underlying cognition are shaped and constrained by the ascending arousal system, whose activity acts to facilitate the integration between internal states and external contingencies (Parvizi & Damasio, 2001). Timely and selective interactions between the ascending arousal system and the network-level configuration of the brain are thus likely to represent crucial constraints on cognitive and attentional processes. Yet, despite these links, we currently have a relatively poor understanding of how the ascending arousal system helps the brain as a whole to functionally reconfigure during cognitive processes, such as attention, in order to facilitate effective cognitive performance.

Recent evidence has linked higher order cognitive functions in the brain to the intersection between whole-brain functional network architecture and the autonomic arousal system (Alnæs et al., 2015; Alnæs et al., 2014; Munn et al., 2021; Shine et al., 2016; Shine, Hearne, et al., 2019). Central to these relationships is the unique neuroanatomy of the ascending noradrenergic system. For instance, the pontine locus coeruleus, which is a major hub of the ascending arousal system, sends widespread projections to the rest of the brain (Samuels & Szabadi, 2008). Upon contact, adrenergic axons release noradrenaline, which acts as a ligand on three types of post- and presynaptic adrenergic receptors (i.e., α1, α2, and β). The functional effects of each of these receptors depend on their differential sensitivities to noradrenaline (affinities for the ligand differ across receptors: α2 > α1 > β) and intracellular cascades, as well as their neuronal and regional distributions (Aston-Jones & Waterhouse, 2016; Bouret & Sara, 2005; Robbins & Arnsten, 2009; Samuels & Szabadi, 2008; Sara, 2009; Shine, 2019). By modulating the excitability of targeted regions, the locus coeruleus can effectively coordinate neural dynamics across large portions of the cerebral cortex (Shine et al., 2021; X. J. Wang, 2020). However, it is challenging to noninvasively track the engagement of the locus coeruleus during whole-brain neuroimaging and cognitive task performance.

Fortunately, it has been widely shown that the pupil diameter directly responds to changes in the activity of the locus coeruleus, and thus serves as an indirect, noninvasive measure of the noradrenergic system (Aston-Jones & Cohen, 2005; S. Joshi et al., 2016). Specifically, pupil diameter has been shown to indirectly monitor the neuromodulatory influences of the ascending arousal system on a variety of different brain regions (Alnæs et al., 2014; Liu et al., 2017; Sara, 2009; van den Brink et al., 2016). Moreover, noradrenergic-mediated dilations in pupil diameter have been shown to effectively track the allocation of attentional resources (Gilzenrat et al., 2010; Kahneman & Beatty, 1966; Wainstein et al., 2017), in addition to both physical and mentally effortful processes (Mulder, 2012; Varazzani et al., 2015). Fast, phasic changes in pupil diameter have also been shown to directly relate to changes in the activity of the locus coeruleus (S. Joshi et al., 2016; Murphy et al., 2016; Reimer et al., 2014). While there is some evidence that pupil diameter covaries with other subcortical systems (S. Joshi & Gold, 2020), such as the cholinergic and serotoninergic system (Cazettes et al., 2020), the physiological mechanism for these effects is more opaque, and there is also clear causal evidence linking stimulation of the locus coeruleus to dilation of the pupil (Liu et al., 2017; Zerbi et al., 2019). Despite these insights, several questions remain unanswered regarding how these processes are related to the complex architecture of the brain (Shenhav et al., 2017). For instance, the processes by which the ascending arousal system modulates the functional dynamics of brain networks to facilitate attention, decision-making, and optimal behavioral performance have only begun to be explored (de Gee et al., 2017; Shine, Breakspear, et al., 2019; Shine et al., 2018; Zerbi et al., 2019).

To examine these relationships in more detail, participants performed a motion-tracking task (top panel of Figure 1A) involving four levels of increasing attentional load, which was modulated by manipulating the number of items required to covertly attend to over an 11-s tracking period. Specifically, subjects were instructed to covertly track the movement of several preidentified targets (two to five) in a field of nontarget stimuli (10 in total, including targets; see Figure 1). To investigate the network topological signatures of performing this task, we collected concurrent BOLD fMRI and pupillometry data. We hypothesized that, if increasing mental effort led to the reconfiguration of large-scale network architecture via the ascending arousal system, then the number of items required to be tracked over time (i.e., the attentional load) should relate to (a) increased pupil diameter; (b) heightened BOLD activity within attentional networks; and (c) augmented topological integration. Also, we predicted that individual differences in pupil diameter should track individual differences in effective attentional performance and decision processes (de Gee et al., 2017; de Gee et al., 2014; Donner et al., 2000). Finally, we tested whether the regional patterns of network configuration were predicted by the distribution of a predefined adrenergic receptor density atlas (Fornito et al., 2019; Richiardi et al., 2015; Shine, Breakspear, et al., 2019; Zerbi et al., 2019). Our results confirm these predictions, and hence provide a mechanistic link between network topology, ascending noradrenergic arousal, and attentional load.

Consistent with previous work (Alnæs et al., 2014), our two-level analysis—linear regression within each subject, and a two-tailed t test between subjects—found that task performance (i.e., correct responses) decreased with attentional load (mean β_Acc = −6.66; t₍₁₇₎ = −5.19, p = 7.2 × 10⁻⁵; Supplementary Figure S1B) while the reaction time (RT) increased with attentional load (mean β_RT = 0.06, t₍₁₇₎ = 5.10, p = 8.8 × 10⁻⁵). We expanded on this result by translating performance into EZ-diffusion model parameters. Roughly, this approach uses the accuracy and reaction time distribution to estimate three latent parameters (de Gee et al., 2014): drift rate, a marker of the accumulation of decision evidence (Equation 1); boundary criteria, the amount of evidence required to make a decision (Equation 2); and non-decision time, the epoch spent processing the tasks perceptually (Equation 3). The advantages of using this model are twofold: first, there are well-known links between the parameters to decision-making processes (Ratcliff et al., 2016; Ratcliff et al., 2015), pupil diameter (Murphy et al., 2016; Murphy et al., 2014) and network reconfiguration (Shine et al., 2016); second, drift rate accounts for the accuracy–reaction time trade-off, as it takes into consideration both accuracy and the variability in reaction time into its calculation. In this way, our approach offers a better approximation of the ongoing computational processing during the task than does accuracy and RT (Ratcliff et al., 2015; Wagenmakers et al., 2007). Using this approach, we observed a decrease in both the boundary criteria (β_Bound = −0.01, t₍₁₇₎ = −2.70, p = 0.015) and the drift rate (mean β_Drift = −0.03, t₍₁₇₎ = −7.43, p = 9.7 × 10⁻⁷; Figure 1B), and an increase in the non-decision time (mean β_nd = 0.07, t₍₁₇₎ = 5.32, p = 5.5 × 10⁻⁵) with increasing attentional load.

By calculating the linear effect of load on pupil size across a moving average window of 160 ms (see Methods), we observed a main effect of increased pupil diameter across both the tracking and the probe epochs (β_pupil > 0, p_FDR < 0.01; light gray area in Figure 1A depicts significant epochs of time during the task; and in Supplementary Figure S1A shows the group average β_pupil time series). We also observed a positive correlation between mean β_pupil during the significant period (for simplicity we will refer to this value as βpupil) to the mean drift rate, mean boundary criteria, and accuracy across all loads (Pearson’s r_drift = 0.8, p = 1.0 × 10⁻⁴; Figure 1C; r_acc = 0.68, p = 1.5 × 10⁻³, Supplementary Figure S1C; r_Bound = 0.71, p = 9 × 10⁻⁴). The same relationships were not observed with non-decision time (Pearson’s r_nd = −0.31, p = 0.19). Additionally, we analyzed whether this effect was present both within and between subjects in a trial-by-trial manner. To this end, we created a logistic linear mixed model (Equation 6) to test whether pupil diameter was a predictor of performance (i.e., correct or incorrect response), as we would expect that incorrect responses should relate to decreased pupil diameter in difficult trials. We used the average pupil diameter within each trial of Load 4 and 5 (to account for the ceiling effect of Load 2 and 3) as regressors and subject as a grouping variable. We found a statistically significant fixed effect of pupil diameter on performance within each trial (β = 0.0127 ± 5 × 10⁻⁴; t₍₂₈₆₎ = 2.48; p = 0.013). Furthermore, we analyzed the random-effect coefficients, which are the dispersion of the regressor across the grouping variable from the fixed regressor (in this case there is one value per subject), to assess the role of average across task performance. We found that the random effect covaried with the average performance and drift rate of each subject (Accuracy: Pearson’s r = 0.73, p = 8 × 10⁻⁵; Drift: Pearson’s r = 0.73, p = 5 × 10⁻⁵), suggesting that trial-by-trial pupil diameter was a better predictor of performance (i.e., correct or incorrect) on subjects with higher average performance in comparison to subjects with lower performance across the task. In conclusion, these results suggest that attentional load manipulation and pupil dilation covaried with performance on this attentionally demanding task both within and between subjects.

Based on previous studies, we hypothesized that an increase in attentional load should recruit a distributed functional network architecture (Alnæs et al., 2014), heightening network integration (Shine, 2019; Shine et al., 2016; Shine, Breakspear, et al., 2019). To test this hypothesis, we implemented a hierarchical topological network analysis (Bassett et al., 2010; Meunier et al., 2010; Meunier et al., 2009) on the average time-resolved functional connectivity matrix calculated across the tracking period of the task. Our analysis identified a subnetwork of tightly interconnected regions that were part of attentional, somatomotor, and cerebellar network (red in Figure 2) that increased its BOLD activity after the tracking onset (Figure 2F). The tightly integrated regions were diversely connected to a separate frontoparietal submodule (blue in Figure 2) that was less active during the trial. Two remaining submodules (yellow and green in Figure 2) showed a negative BOLD response during the tracking period and were part of a diverse set of networks. Interestingly, 81% of the frontoparietal network (FPN) and all the default mode network (DMN) were found to be within this less active group (see Supplementary Table S2 for the complete list of regions and submodule assignments).

Contrary to expectations, we did not observe significant parametric topological change (i.e., modularity, Q) at the macroscopic level as a function of attentional load (p > 0.05 for all TRs, Supplementary Figure S2A). However, when analyzing the correlation between modularity and performance measures (i.e., accuracy, drift rate, and pupil diameter), we observed that an increase in the large-scale modularity load effect (i.e., higher modularity with load, β_QL) positively correlated with higher mean drift rate (Pearson’s r = 0.53; p = 0.022; Figure 2D), mean accuracy (Pearson’s r = 0.61; p = 0.007; Supplementary Figure S3A), but was independent from βpupil (Pearson’s r = 0.43; p = 0.073). These results suggested that the system reconfigured during tracking towards increasing modularity, which in turn affected the efficient encoding of the ongoing task during tracking and hence, the decision-making process during the task probe.

Upon closer inspection of the data (Figure 2C), we observed a substantial number of nodes that were playing an integrative role during task performance, albeit at a finer resolution than the initial analysis suggested. We performed the modularity assignment within each large-scale module. The hierarchical analysis resulted in two pairs of submodules at the mesoscale level with a significant modularity (compared with 100 random graphs with preserved signed degree distribution; Q_M1 = 0.137, p = 0; Q_M2 = 0.137, p = 0; Figure 2E). Specifically, the red submodule was found to selectively increase its participation coefficient (PC) at the mesoscale level (i.e., by increasing the connection weights to the blue submodule in comparison with intramodular connections; Equation 5) as a function of increasing attentional load (β_PC = 2.4 × 10⁻³, t₍₁₇₎ = 3.57; p = 0.002; Figure 3A). Additionally, the extent of integration in the red submodule was positively correlated across subjects with βpupil (Pearson’s r = 0.62, p = 0.006; Figure 3B), drift rate (Pearson’s r = 0.66, p = 0.002; Figure 3C), and accuracy (Pearson’s r = 0.57, p = 0.012, Supplementary Figure S3B). Importantly, these relationships were found to be specific to the red submodule (blue: Pearson’s r = −0.02, p = 0.936; yellow: Pearson’s r = −0.011, p = 0.965; green: Pearson’s r = −0.12, p = 0.617).

Based on these results, we implemented a linear mixed model (Equation 7), using the subjects’ average pupil response within each load as a regressor and the average participation of the red submodule as the dependent variable, with grouping by subject. Using this approach, we observed a significant fixed effect of pupil diameter on PC (β = 7.6 × 10⁻³ ± 3 × 10⁻³, t₍₇₀₎ = 2.60, p = 0.011). Furthermore, the random-effect coefficients (i.e., the between-subject variation of the regressor value) correlated positively with accuracy (Pearson’s r = 0.47, p = 0.048) and drift rate (Pearson’s r = 0.62, p = 0.005), suggesting that subjects with a strong relationship between red module integration and pupil diameter have better behavioral outcomes. We then correlated the red βPC to the load effect on large-scale modularity (βQ_L, Figure 2D) and observed a significant positive correlation (Pearson’s r = 0.59, p = 0.009). Finally, given that both topological parameters were correlated with drift rate and also with each other, we performed a partial correlation between drift rate and βPC controlling by βQ_L (r = 0.51, p = 0.034), and the partial correlation between drift rate and βQ_L controlling by βPC (r = 0.36, p = 0.145). This suggests that drift rate is correlated with the mesoscale integration of the red submodule, but less so with increases in large-scale modularity. Thus, although the macroscale network did not demonstrate increased integration per se, the relative amount of mesoscale integration within the red community was associated with increased performance (i.e., drift rate) and sympathetic arousal (i.e., pupil diameter), both between and within subjects. In this way, these results provide a direct relationship between the effect of attention load on pupillometry, drift rate, and a trade-off between large-scale segregation and mesoscale network integration.

Given the relationship between mental effort, noradrenergic tone, and pupil dilation (Alnæs et al., 2014; S. Joshi et al., 2016; McGinley et al., 2015; Reimer et al., 2014; Varazzani et al., 2015), the results of our analyses strongly suggested that the adrenergic system is involved in the mesoscale network reconfiguration observed during attentional tracking. The locus coeruleus can impact the cortical system in multiple ways, both through direct release of noradrenaline onto cortical neurons, and through the modulation of subcortical regions (such as the thalamic nuclei) with concurrent impact on the cortical dynamic. Importantly, in either case, the modulation is dependent on the noradrenergic receptors subtypes, which have different sensitivities to noradrenaline (Robbins & Arnsten, 2009; M. Wang et al., 2007) and variable expression in the cerebral cortex (Santana & Artigas, 2017; Zilles & Palomero-Gallagher, 2017), and also belong to distinct classes (i.e., α1, α2, and β receptors). In particular, the α2a has been previously associated with working memory, adaptive gain, and effective attention (Arnsten et al., 2012; Robbins & Arnsten, 2009; M. Wang et al., 2007). To gain a deeper insight into the role of α2a receptors in mesoscale integration during attentional tracking, we extracted the regional expression of the ADRA2A gene (which codes for α2a adrenoceptors) from the Allen Human Brain Atlas repository (Gryglewski et al., 2018; Hawrylycz et al., 2012), and compared the cortical regional expression of this gene with the brain activity patterns identified in our network analysis (Figure 2E).

Based on the relationships between pupil diameter (Figure 1), topological signatures (Figure 2), and task performance (Figure 3), and the known link between these variables and engagement of the noradrenergic system, we hypothesized that the different modules and submodules that we observed should have different densities of neuromodulatory receptors to account for the differential patterns across the network. To test this hypothesis, we conducted a two-tailed t test in each hierarchical level comparing the density of the ADRA2A expression between modules. To account for spatial autocorrelation, we generated 5,000 surrogate maps with the same spatial autocorrelation of the ADRA2A map, calculated a t statistic for each surrogate, and evaluated the probability of finding the observed t statistic against the null distribution (Burt et al., 2020; Markello & Misic, 2021). We indeed observed significant differences between modules at the mesoscale level. Specifically, we found significant differences between the blue and yellow submodules (t₍₁₉₄₎ = 3.82, p = 2 × 10⁻⁴, p_SA = 0.02) and the differences between green and yellow submodules (t₍₁₇₇₎ = −4.47, p = 1.3 × 10⁻⁵, p_SA = 0.004), while the other differences did not survive the spatial autocorrelation test (green-red: t₍₁₅₂₎ = 0.47; p = 0.635, p_SA = 0.590; yellow-red: t₍₁₅₆₎ = −3.02, p = 0.003, p_SA = 0.121; green-blue: t₍₁₇₃₎ = −0.68, p = 0.496, p_SA = 0.324; red-blue: t₍₁₃₅₎ = −1.30, p = 0.195, p_SA = 0.237; Supplementary Figure S5A).

The modulatory effects of noradrenaline have been argued to depend directly on ongoing glutamatergic activity in target regions (Mather et al., 2016; Shine, 2021). Moreover, it has been shown that the main source of the BOLD activity is the neurovascular response caused by pyramidal neurons containing cyclo-oxygenase-2 (Lecrux & Hamel, 2016). Importantly, this evoked response following noradrenergic activation is dependent on the ongoing activity of the pyramidal neurons (Bekar et al., 2012). Thus, the role of noradrenaline on brain dynamics and BOLD response depends critically on ongoing glutamatergic activity, which putatively represents pooled neural spiking activity (Logothetis, 2003). Given the differential task-related BOLD activity of the different submodules (i.e., Figure 2F, Supplementary Figure S4, and Figure 4A), and the observed regional variability and specificity of integration across the network, we hypothesized that network-level integration would be explained by the combined effect of ongoing BOLD activity and the distribution of the adrenergic receptor expression. Finally, we predicted that the role of the α2a receptor atlas in shaping brain activity and topology should be dependent of the subjects’ pupil diameter, such that higher βpupil should rely on a stronger relationship between network topology and α2a receptor expression.

To evaluate between these different hypotheses, we created three linear mixed models in order to better disentangle the different plausible interactions between the variables (see Methods), while still controlling for between-subject variability as a grouping variable. Additionally, to control for spatial autocorrelation, we used 5,000 surrogate maps that maintained the spatial autocorrelation of the α2a while permuting the density values. In the first model (Equation 8), we tested the hypothesis that the parametric BOLD effect (i.e., βBOLD, Supplementary Figure S4) is shaped by the distribution of α2a receptors. We found significant evidence for a positive fixed effect of α2a on βBOLD activity; however, this effect did not survive correction for spatial autocorrelation (β_α2a = 0.037 ± 0.016; t₍₅₉₉₂₎ = 2.29; p = 0.022; p_SA = 0.106; Supplementary Table S1). Furthermore, we correlated the random-effect coefficients (from the original and the surrogate maps) to both βPC and βpupil and observed a significant positive correlation between the participation coefficient and both pupils (Pearson’s r = 0.54, p = 0.02, p_SA = 0.036; Figure 4B) and mean drift rate (Pearson’s r = 0.70, p = 0.001, p_SA = 0.001; Figure 4C). This result shows the manner in which pupil diameter linearly shapes βBOLD cortical map through the engagement of the α2a receptor expression map. Importantly, although the fixed effect of α2a on βBOLD didn’t survive the spatial autocorrelation correction, the linear correlation of this effect with both βpupil and drift rate (between subjects) did survive the correction.

To further analyze the between-subject differences in the role of α2a receptor atlas in shaping the βBOLD map, we ran a separate linear model within each subject with α2a as a regressor and βBOLD of each region as the dependent variable (while also correcting for spatial autocorrelation using 5,000 surrogate maps). As can be seen in Figure 4B–C, we observed a dependency between the p_SA value, βpupil, and drift rate, in which the respective within-subject effects that survived the spatial autocorrelation correction are shown (p_SA < 0.05; marked circles in Figure 4B–C). Despite these results, there was no significant effect of α2a on βPC (Equation 9; β_α2a = 0.001 ± 0.003; t₍₅₉₉₂₎ = −0.51; p = 0.6), and no significant Pearson’s correlations were found between the random effects and both βpupil or drift rate (r = −0.24, p = 0.33; and r = −0.23, p = 0.341, respectively). However, we did find a significant effect of βBOLD on βPC (Equation 10; β = 0.0259 ± 0.006; t₍₅₉₉₂₎ = 3.96; p = 7.55 × 10⁻⁵; Supplementary Table S1 and Figure 4D). Together these results propose a closer link between pupil diameter, ascending neuromodulation, and the cortical neuromodulation dependent on α2a receptor density.

Finally, we observed a differential relationship between βPC and βBOLD depending on the large-scale module to which the regions were assigned. We expanded the former result by measuring, within each subject, the Pearson correlation between the βBOLD and βPC separately in each large-scale module (M1 being the modules assigned as red and blue, and M2 assigned as yellow and green; Figure 2). The results demonstrated a significant difference between modules, meaning that M1 has a higher correlation with βPC, in comparison with M2 (t₍₁₇₎ = −12.99, p = 2.93 × 10⁻¹⁰, Supplementary Figure S5C). These results provided evidence that the adrenergic receptor distribution of α2a shapes the βBOLD activation map in proportion to the subject’s pupil diameter. Additionally, βBOLD activation map modulates (i.e., was related to) mesoscale integration, and mesoscale integration is related to pupil diameter. Based on these results, we hypothesize that the adrenergic system shapes the BOLD activity, which in turns shapes the topology of the network towards integration. However, future work is required in order to test this hypothesis more directly, for instance by combining optogenetic approaches with neuronal recordings in awake animals.

Here, we leveraged a unique dataset to simultaneously track pupil diameter and network topology during an attentional demanding task with increasing attentional load. Our results provide integrative evidence that links the ascending arousal system to the mesoscale topological signature of the functional brain network during the processing of an attentionally demanding cognitive task. Pupil diameter was tracked with attentional load (Figure 1A) and was related to the speed of information accumulation as estimated by a drift diffusion model (Figure 1B–C). Additionally, we observed concurrent pupil dilations and adaptive mesoscale parametric topological changes as a function of task demands (Figures 2 and 3). Finally, we found evidence that topological reconfiguration was dependent on the regional activity and the genetic expression of the adrenergic receptors in the brain (Figure 4). Together, these results provide evidence for the manner in which the ascending arousal noradrenergic system reconfigures brain network topology so as to promote attentional performance according to task demands.

The relationship between performance and pupil diameter is consistent with the predictions of adaptive gain theory (Aston-Jones & Cohen, 2005). Within this framework, the locus coeruleus is proposed to adaptively alter its activity according to the demands imposed on the system. More specifically, the theory proposes that performance follows an inverted U-shaped relationship with arousal, such that maximal operational flexibility in the noradrenergic system is associated with optimal task performance (Arnsten et al., 2012; Robbins & Arnsten, 2009). We observed that load-related increases in pupil diameter, presumably due to increased activity in the ascending arousal system (Aston-Jones & Cohen, 2005; S. Joshi et al., 2016; Liu et al., 2017), relates closely with the activity and topology of the broader brain network (Figure 2), in a manner that is reflective of effective task performance (Figure 3). Similar effects have been described in animal models after a chemogenetic activation of the locus coeruleus, which strongly alters the large-scale network structure towards large-scale integration, specifically in regions with heightened adrenergic receptor expression (Zerbi et al., 2019). How these changes, which are likely related to the modulation of the neural gain that mediates effective connections between distributed regions of the brain (Shine et al., 2021; Shine et al., 2018), are traded off against requirements for specificity and flexibility remains an important open question for future research.

The addition of attentional load was found to alter the integration of mesoscale submodules, but not the higher level modular organization. This topological result is somewhat more targeted than those described in previous work (Shine, Breakspear, et al., 2019; Shine et al., 2016). While these differences may be related to disparities in the way that the data were analyzed, the results of our study do demonstrate that alterations in the cerebral network topology at a relatively local (i.e., submodular) level are crucial for effective task performance (Akiki & Abdallah, 2019). Additionally, our results replicate and expand upon a previous study (Mohr et al., 2016), in which the authors found that short-term practice on an attentional task was related to increased coupling between attentional networks and segregation among task-negative (DMN) and frontoparietal network (FPN). Our study replicates the graph theoretical results of that study, while also directly relating the findings to the architecture of the ascending neuromodulatory system. One potential explanation for these results comes from animal studies, in which rapid changes in pupil diameter have been compared with changes in neural population activity at the microscale (S. Joshi et al., 2016; McGinley et al., 2015; Reimer et al., 2014). These studies suggest that the ascending arousal system may be able to alter the topology of the network in a hierarchical manner that is commensurate with the spatiotemporal scale of the arousal systems’ capacity (Shine et al., 2016). Future work that integrates results across spatiotemporal scales is required to appropriately adjudicate the implications of this hypothesis.

Importantly, our approach is not without limitations. For one, the participation measures used in our linear mixed model were estimated at the mesoscale level, and hence derived from different modular partitions. Furthermore, the specificity of the pupillary response as a correlate of locus coeruleus (LC) activity is currently under active debate. For instance, in addition to the strong empirical links between the noradrenergic system and pupil dilation, there is also evidence that the pupil is dilated in concert with activity in the basal forebrain cholinergic system (Reimer et al., 2016), however it bears mention that both peripheral (Kaymak et al., 2018) and central cholinergic tone (Yüzgeç et al., 2018) are associated with pupillary constriction. There are more plausible physiological routes for the serotonergic system to dilate the pupil (via the excitation of the intermediolateral cell column), and in keeping with this, there is evidence that the serotonergic system is linked with pupil dilation (Cazettes et al., 2020). Nevertheless, it is important to take into account that the neuromodulatory arousal system is replete with complex interconnections (Avery & Krichmar, 2017; Briand et al., 2007; A. Joshi et al., 2017; Smiley et al., 1999). In addition, based on the current lack of a specific mechanism involving pupillary changes through the cholinergic system, it is highly probable that those correlations are due to indirect modulation of pupillary responses (e.g., via indirect neuromodulation mediated by the LC system). On the other hand, we acknowledge the limitations of the atlas receptor analysis and the linear model used in our study. More specific neurobiological properties of the receptor distributions are needed to make better inferences, and hence provide more accurate answers of their role in brain dynamics. For instance, it would be ideal to compare receptor distributions that incorporated layer-specific expression, as there are well-known cellular and circuit differences across layers in the cerebral cortex (Douglas & Martin, 2004; Palomero-Gallagher & Zilles, 2019). Importantly, taking into consideration the strong correlation between different genetic expression maps (Fornito et al., 2019), it is possible that the current correlation between ADRA2A expression and brain activity is a false positive caused by another neuroanatomical gradient strongly correlated to the ADRA2A. Therefore, future work studying the interaction between genetic expression of the neuromodulatory receptors, pupil diameter, and brain activity is needed. In spite of this limitation, we believe in the importance of integrating pupil diameter and receptor distribution in the analysis as the relationships between noradrenergic tone, brain activity, and network topology will help us to disentangle the mechanistic steps connecting the locus coeruleus system to both pupil diameter and brain dynamics.

In summary, we provide evidence linking mesoscale topological network integration, hierarchical organization, and BOLD dynamics in the human brain that increases in attentional load, thus providing further mechanistic clarity over the processes that underpin the adaptive gain model of noradrenergic function in the central nervous system.

Eighteen right-handed individuals (age 19–26 years; five male) were included in this study. Exclusion criteria included standard contraindications for MRI; neurological disorders; and mental disorders or drug abuse. All participants gave written informed consent before the experiment.

Each trial of the task involved the same basic pattern (Figure 1A): The task begins with a display presenting the objects (i.e., blue colored disks); after a 2.5-s delay, a subset of the disks turn red for another 2.5 s; all of the disks then return to blue (2.5 s) before they start moving randomly inside the tracking area. The participants’ job is to track the “target” dots on the screen while visually fixating at the cross located at the center of the screen. After a tracking period of ∼11 s, one of the disks is highlighted in green (a “probe”) and the subject is then asked to respond, as quickly as possible, as to whether the green probe object was one of the original target objects. The number of objects that subjects were required to attend to across the tracking period varied across trials. There were five trial types: passive viewing (PV), in which no target is assigned; and four load conditions, in which two to five targets were assigned for tracking. We operationalized attentional load as the linear effect of increasing task difficulty (i.e., the number of targets to be tracked).

The experiment was conducted using a blocked design, in which each block included the following: instruction (1 s); fixation (0.3 s, present throughout the rest of trial); object presentation (all objects were blue; 2.5 s); target assignment (i.e., the targets changed color from blue to red; 2.5 s); object representation (objects back to the original blue color; 2.5 s); object movement/attentional tracking (moving blue dots; 11 s); object movement cessation (0.5 s); and a final probe (color change to green and response; 2.5 s). The total duration of each trial was 22.8 s. Each condition was repeated four times in one fMRI run, which also included four separate fixation periods of 11 s each between five consecutive trials. All participants completed four separate runs of the experiment, each of which comprised 267 volumes. The order of the conditions was pseudorandom, such that the different conditions were grouped in sub-runs of triplets: PV, pseudorandom blocks of Loads 2 through 5, and a fixation trial. All objects were identical during the tracking interval and standard object colors were isoluminant (to minimize incidental pupillary responses during the task).

Fluctuations in pupil diameter of the left eye were collected using an MR-compatible coil-mounted infrared EyeTracking system (NNL EyeTracking camera, NordicNeuroLab, Bergen, Norway), at a sampling rate of 60 Hz and recorded using the iView X Software (SensoMotoric Instruments, SMI GmbH, Germany). Blinks, artifacts, and outliers were removed and linearly interpolated (Wainstein et al., 2017). High-frequency noise was smoothed using a second-order 2.5-Hz low-pass Butterworth filter. To obtain the pupil diameter average profile for each level of attentional load (Figure 1B), data from each participant were normalized across each task block (corresponding to the five consecutive trials between fixations). This allowed us to correct for low-frequency baseline changes without eliminating the load effect and baseline differences due to load manipulations (Campos-Arteaga et al., 2020; Rojas-Líbano et al., 2019). Following this, a linear regression was performed in each time point using the task load as regressor and resulting in a “load effect” time series for each subject.

Imaging data were collected on a Philips Achieva 3 Tesla MR-scanner, equipped with an eight-channel Philips SENSE head coil (Philips Medical Systems, Best, Netherlands) at the Intervention Centre, Oslo University Hospital, Norway. Functional data were collected using a BOLD-sensitive T2*-weighted echo-planar imaging sequence (36 slices, no gap; repetition time (TR), 2,2 s; echo time (TE), 30 ms; flip angle, 80°; voxel size, 3 × 3 × 3; field of view (FOV), 240 × 240 mm; interleaved acquisition). Anatomical T1-weighted images consisting of 180 sagittal-oriented slices were obtained using a turbo field echo pulse sequence (TR, 6.7 ms; TE, 3.1 ms; flip angle 8°; voxel size 1 × 1.2 × 1.2 mm; FOV, 256 × 256 mm).

After realignment (using FSL’s MCFLIRT), we used FEAT to unwarp the EPI images in the y-direction with a 10% signal loss threshold and an effective echo spacing of 0.333. Following noise-cleaning with FIX (custom training set for scanner, threshold 20, included regression of estimated motion parameters), the unwarped EPI images were then smoothed at 6-mm FWHM, and nonlinearly coregistered with the anatomical T1 to 2-mm isotropic MNI space. Temporal artifacts were identified in each dataset by calculating framewise displacement (FD) from the derivatives of the six rigid-body realignment parameters estimated during standard volume realignment (Power et al., 2014), as well as the root mean square change in BOLD signal from volume to volume (DVARS). Frames associated with FD > 0.25 mm or DVARS > 2.5% were identified; however, as no participants were identified with greater than 10% of the resting time points exceeding these values, no trials were excluded from further analysis. There were no differences in head motion parameters between the four sessions (p > 0.500). Following artifact detection, nuisance covariates associated with the six linear head movement parameters (and their temporal derivatives), DVARS, physiological regressors (created using the RETROICOR method), and anatomical masks from the cerebrospinal fluid and deep cerebral white matter were regressed from the data using the CompCor strategy (Behzadi et al., 2007). Finally, in keeping with previous time-resolved connectivity experiments (Gu et al., 2015), a temporal band pass filter (0.0071 < f < 0.125 Hz) was applied to the data.

Following preprocessing, the mean time series was extracted from 375 predefined regions of interest (ROIs). To ensure whole-brain coverage, we extracted the following: (a) 333 cortical parcels (161 and 162 regions from the left and right hemispheres, respectively) using the Gordon atlas (Gordon et al., 2016); (b) 14 subcortical regions from the Harvard-Oxford subcortical atlas (bilateral thalamus, caudate, putamen, ventral striatum, globus pallidus, amygdala, and hippocampus; https://fsl.fmrib.ox.ac.uk/); and (c) 28 cerebellar regions from the SUIT atlas (Diedrichsen et al., 2009) for each participant in the study.

To estimate functional connectivity between the 375 ROIs, we used the jackknife correlation (JC) approach (Thompson et al., 2018). Briefly, this approach estimates the static correlations between each pair of regions, and then recalculates the correlation between each pair after systematically removing each temporal “slice” of data (i.e., each TR). By subtracting the jackknifed correlation matrix from the original “static” matrix, the difference in connectivity at each slice from the static connectivity value can be used as an estimate of time-resolved functional connectivity between each pair of regions at each TR in a way that does not require windowing.

For each subject, we calculated the mean adjacency matrix from 1 TR before tracking until the end of the tracking period. Afterwards, a consensus partition was estimated across subjects. Finally, to identify multilevel structure in our data, we repeated the modularity analysis for each of the modules identified in the first step (Meunier et al., 2010; Meunier et al., 2009). With this final module assignment, we were afforded an estimate of the time-resolved, multilevel modularity (Q_T) within each temporal window for each participant in the study.

To investigate the potential correlates of mesoscale integration, we interrogated the neurotransmitter receptor signature of each region of the brain. We used the Allen Brain Atlas microarray atlas dataset (https://human.brain-map.org/; Hawrylycz et al., 2012) to identify the regional signature of genetic expression of the α2a subtype of the adrenergic receptor (ADRA2A). This receptor has been a priori related to cognitive function and attention (Arnsten & Haven, 2013), and is one of the most abundant adrenergic subtypes expressed in the cerebral cortex (Perez, 2020). This atlas contains postmortem samples of six donors that underwent microarray transcriptional characterization. The spatial map of α2a mRNA expression was obtained in volumetric 2-mm isotropic MNI space, following improved nonlinear registration and whole-brain prediction using variogram modeling (Gryglewski et al., 2018). We used these data instead of the native sample-wise values in the AHBA database to prevent bias that could occur because of spatial inhomogeneity of the sampled locations. We projected the volumetric α2a expression data onto the Gordon atlas with linear interpolation and calculated the mean value within each parcel using custom MATLAB codes.

We analyzed the between-subject effect of load on the behavioral, pupillometric, and fMRI-related variables by performing a two-level linear model analysis. In the first level, we used attentional load as a regressor (2 to 5) and—in independent models—the mean accuracy, reaction time, standard deviation of reaction time, drift rate, boundary criteria, and non-decision time as dependent variables (i.e., four values per subject). From this, we ran a two-tailed t test on the statistical effects (i.e., the β value from the regression, one for each subject; N = 18). Similarly, to calculate the load effect on pupil diameter, we calculated the average pupil diameter on each load condition within each subject. Then, we performed a first-level analysis in which we ran a linear regression in each time frame (1600 frames in total, corresponding to 26.6 seconds). This procedure resulted in one β timeseries (i.e., the statistical load effect on pupil diameter) for each subject across the trial (Supplementary Figure S1A). After this, we performed a right tailed t test in each frame across subjects (n = 18 in each frame) to find the periods of time where the β values where higher than 0. Finally, we corrected by false discovery rate (FDR; Benjamini & Yekutieli, 2001) for multiple testing, which resulted in a period of time in which the load effect was higher than 0 (light gray area in Figure 1A). The mean β values during this section was calculated in each subject and defined as “βpupil.” Finally, following the same pipeline, we calculated the effect of attentional load on the brain-related signals (i.e., BOLD, participation coefficient [PC], and modularity [Q]). The effect of load on BOLD was calculated running a separate linear model in each subject and region within each TR (18 subjects; 375 regions; 10 TRs; 4 load conditions), resulting in a matrix of β values of 18 × 375 × 10.

To evaluate whether the modularity of the network we observed was higher than chance, we generated 100 random networks in each hierarchical level (300 random networks in total), with a preserved degree distribution (using the MATLAB randmio_und_signed function from the Brain Connectivity Toolbox; Rubinov & Sporns, 2010). We calculated the modularity value of each random network and used the resultant values to populate a null distribution (Figure 2D).

Expression of brain genetic atlas varies smoothly across the surface and thus is associated with nontrivial spatial autocorrelation that in turn violates the assumption of independence between samples (Burt et al., 2020; Markello & Misic, 2021; Vos de Wael et al., 2020). To account for the spatial autocorrelation in these brain maps, we used spatial autocorrelation null maps as implemented in Brain Surrogate Maps with Autocorrelated Spatial Heterogeneity (BrainSMASH) Python toolbox (Burt et al., 2020). A geodesic distance matrix of the atlas parcels using the surface of the Gordon atlas was obtained to build the surrogates using BrainSMASH functions. We generated 5,000 null maps that were used to generate null distribution of the different statistics corrected by spatial autocorrelation.

We thank P. Billeke for his thoughtful comments on our manuscript.

The anonymized preprocessed fMRI and pupillometry data can be found at https://figshare.com/articles/dataset/MOT_data_mat/13244504 (Wainstein et al., 2020). The ADRA2A expression atlas can be downloaded from https://www.meduniwien.ac.at/neuroimaging/mRNA.html. All analysis of the fMRI and pupil diameter data were performed on MATLAB 2020a. The surrogate maps of the ADRA2A atlas were generated on Python. Documented code for reproducing the analyses is provided in https://github.com/gabwainstein/MOT (Wainstein, 2021). Supporting information for this article is available at https://doi.org/10.1162/netn_a_00205.

Gabriel Wainstein: Conceptualization; Formal analysis; Investigation; Methodology; Visualization; Writing – original draft; Writing – review & editing. Daniel Rojas-Líbano: Visualization; Writing – original draft; Writing – review & editing. Vicente Medel: Formal analysis; Writing – review & editing. Dag Alnæs: Conceptualization; Data curation; Funding acquisition; Investigation; Methodology; Writing – review & editing. Knut K. Kolskår: Conceptualization; Investigation; Validation; Writing – review & editing. Tor Endestad: Conceptualization; Investigation; Methodology; Writing – review & editing. Bruno Laeng: Conceptualization; Investigation; Writing – review & editing. Tomas Ossandon: Supervision; Writing – review & editing. Nicolás Crossley: Validation; Writing – review & editing. Elie Matar: Conceptualization; Writing – review & editing. James M. Shine: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Project administration; Resources; Supervision; Validation; Writing – original draft; Writing – review & editing.

James M. Shine, the University of Sydney Robinson Fellowship. James M. Shine, National Health and Medical Research Council (https://dx.doi.org/10.13039/501100000925), Award ID: GNT1156536. Gabriel Wainstein, Becas Chile.

Network:

Theoretical representation of a system. Each part is represented as a node and the connection between nodes as edges.

Integration:

Number of connections of a region or set of regions outside its own module.

Locus coeruleus:

Principal noradrenergic neuromodulatory nuclei located in the brain stem. Projects towards the central nervous system.

Noradrenaline:

A monoaminergic neurotransmitter that modulates the neuronal activity of the target populations in the nervous system.

Drift rate:

A marker of the speed of the accumulation of decision evidence during the decision-making process.

Module:

Group of nodes that have more connection strength between them than to the rest of the system.

Modularity:

Quality function that represents the amount of connection of nodes within each module to nodes in the same module.

Load effect:

Statistical linear effect of the task load (i.e., difficulty) on a given dependent variable.

α2a:

Noradrenergic G-protein coupled receptor. It is the main adrenergic presynaptic and postsynaptic receptor subtype in the human brain.

ADRA2A atlas:

RNA expression atlas of the α2a. It is used as an indirect marker of the density of the α2a receptor.

https://github.com/gabwainstein/MOT

https://doi.org/10.6084/m9.figshare.13244504.v1