Previous research has shown that the autonomic nervous system provides essential constraints over ongoing cognitive function. However, there is currently a relative lack of direct empirical evidence for how this interaction manifests in the brain at the macroscale level. Here, we examine the role of ascending arousal and attentional load on large-scale network dynamics by combining pupillometry, functional MRI, and graph theoretical analysis to analyze data from a visual motion-tracking task with a parametric load manipulation. We found that attentional load effects were observable in measures of pupil diameter and in a set of brain regions that parametrically modulated their BOLD activity and mesoscale network-level integration. In addition, the regional patterns of network reconfiguration were correlated with the spatial distribution of the α2a adrenergic receptor. Our results further solidify the relationship between ascending noradrenergic activity, large-scale network integration, and cognitive task performance.

In our daily lives, it is usual to encounter highly demanding cognitive tasks. They have been traditionally regarded as challenges that are solved mainly through cerebral activity, specifically via information-processing steps carried by neurons in the cerebral cortex. Activity in cortical networks thus constitutes a key factor for improving our understanding of cognitive processes. However, recent evidence has shown that evolutionary older players in the central nervous system, such as brain stem’s ascending modulatory systems, might play an equally important role in diverse cognitive mechanisms. Our article examines the role of the ascending arousal system on large-scale network dynamics by combining pupillometry, functional MRI, and graph theoretical analysis.

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.

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.

Figure 1.

Effect of task difficulty on pupil diameter. (A) Group average (z-score) pupil diameter time series for each load condition. Colors represent passive viewing (PV) in blue, and Loads 2 to 5 in green, orange, red, and black, respectively. The shaded area represents the standard error of the mean. We observed an average increase in pupil diameter, during tracking, with each load condition. The light gray area represents time points with significant parametric effect (βpupil > 0; FDR corrected at p < 0.01). Dotted lines represent the onset of each trial event (shown in the top part of the figure). The red dotted line (time = 0) is the tracking onset period when the dots began to move. (B) Drift rate in each load condition. Each dot is the drift rate for each subject and load (mean βDrift = −0.03, t(17) = −7.43, p = 9.7 × 10−7). (C) Pearson correlation between the pupil parametric effect of load (βpupil) with the average drift rate across subjects (rdrift = 0.8, p = 1.0 × 10−4). The x-axis is the mean beta estimate of the pupillary load effect of the significative time window (βpupil), and the y-axis represents the mean drift rate across loads.

Figure 1.

Effect of task difficulty on pupil diameter. (A) Group average (z-score) pupil diameter time series for each load condition. Colors represent passive viewing (PV) in blue, and Loads 2 to 5 in green, orange, red, and black, respectively. The shaded area represents the standard error of the mean. We observed an average increase in pupil diameter, during tracking, with each load condition. The light gray area represents time points with significant parametric effect (βpupil > 0; FDR corrected at p < 0.01). Dotted lines represent the onset of each trial event (shown in the top part of the figure). The red dotted line (time = 0) is the tracking onset period when the dots began to move. (B) Drift rate in each load condition. Each dot is the drift rate for each subject and load (mean βDrift = −0.03, t(17) = −7.43, p = 9.7 × 10−7). (C) Pearson correlation between the pupil parametric effect of load (βpupil) with the average drift rate across subjects (rdrift = 0.8, p = 1.0 × 10−4). The x-axis is the mean beta estimate of the pupillary load effect of the significative time window (βpupil), and the y-axis represents the mean drift rate across loads.

Close modal

### The Relationship Between Sympathetic Tone and Attentional Processing

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(17) = −5.19, p = 7.2 × 10−5; Supplementary Figure S1B) while the reaction time (RT) increased with attentional load (mean βRT = 0.06, t(17) = 5.10, p = 8.8 × 10−5). 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(17) = −2.70, p = 0.015) and the drift rate (mean βDrift = −0.03, t(17) = −7.43, p = 9.7 × 10−7; Figure 1B), and an increase in the non-decision time (mean βnd = 0.07, t(17) = 5.32, p = 5.5 × 10−5) 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, pFDR < 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 rdrift = 0.8, p = 1.0 × 10−4; Figure 1C; racc = 0.68, p = 1.5 × 10−3, Supplementary Figure S1C; rBound = 0.71, p = 9 × 10−4). The same relationships were not observed with non-decision time (Pearson’s rnd = −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−4; t(286) = 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−5; Drift: Pearson’s r = 0.73, p = 5 × 10−5), 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.

### Network Integration Increases as a Function of Attentional Load

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).

Figure 2.

Hierarchical functional topology analysis of the brain during tracking across all loads. We observed two large-scale modules and two mesoscale modules within each larger module (Module 1 [M1, red/blue] and Module 2 [M2, green/yellow], respectively): M1 corresponded to predominantly attentional and somatomotor network, and M2 to frontoparietal network (FPN) and default mode network (DMN), among others (B and E). (A) Forced directed plot representation of the average cluster across subjects. Edges stronger than 0.3 are shown. Each color represents a unique submodule. (B) A circle plot representing the resting-state regions that were included within each submodule, with networks with >30% of regions in each submodule shown in the plot. The diameter of the circles corresponds to the percentage of network regions that participated in that cluster. Connection width relates to average positive connection strength (functional connectivity); however, only connections with r > 0.1 are shown. (C) Connectivity matrix (Pearson’s r) between all pairs of regions ordered by module assignments—note the strong anticorrelation between the red and green/yellow submodules. (D) Correlation between parametric load effect on large-scale modularity (βQ value) and drift rate (Pearson’s r = 0.53; p = 0.022). (E) Hierarchical analysis representation: QL, QM1, and QM2 represent the modularity value for each level (QL large-scale, and QM1–M2 mesoscale level), and ** represents the probability of finding this value when running a null model (p = 0 for all three modularity values). The brain maps correspond to the cortical regions associated with each submodule. (F) BOLD mean effect for each subcluster. Each line represents the group average, and shaded areas are the standard error of the mean. X-axis is repetition time (TR) centered around tracking onset (TR = 0). DAN, dorsal attention; VN, visual; FPN, frontoparietal; SN, salience; CO, cingulo-opercular; VAN, ventral attention; SMm, somatomotor mouth; SMh, somatomotor hand; RSpN, retrosplenial; FTP, frontotemporal; DMN, default mode; AN, auditory; CPN, cinguloparietal; SubC, subcortex; Cer, cerebellar.

Figure 2.

Hierarchical functional topology analysis of the brain during tracking across all loads. We observed two large-scale modules and two mesoscale modules within each larger module (Module 1 [M1, red/blue] and Module 2 [M2, green/yellow], respectively): M1 corresponded to predominantly attentional and somatomotor network, and M2 to frontoparietal network (FPN) and default mode network (DMN), among others (B and E). (A) Forced directed plot representation of the average cluster across subjects. Edges stronger than 0.3 are shown. Each color represents a unique submodule. (B) A circle plot representing the resting-state regions that were included within each submodule, with networks with >30% of regions in each submodule shown in the plot. The diameter of the circles corresponds to the percentage of network regions that participated in that cluster. Connection width relates to average positive connection strength (functional connectivity); however, only connections with r > 0.1 are shown. (C) Connectivity matrix (Pearson’s r) between all pairs of regions ordered by module assignments—note the strong anticorrelation between the red and green/yellow submodules. (D) Correlation between parametric load effect on large-scale modularity (βQ value) and drift rate (Pearson’s r = 0.53; p = 0.022). (E) Hierarchical analysis representation: QL, QM1, and QM2 represent the modularity value for each level (QL large-scale, and QM1–M2 mesoscale level), and ** represents the probability of finding this value when running a null model (p = 0 for all three modularity values). The brain maps correspond to the cortical regions associated with each submodule. (F) BOLD mean effect for each subcluster. Each line represents the group average, and shaded areas are the standard error of the mean. X-axis is repetition time (TR) centered around tracking onset (TR = 0). DAN, dorsal attention; VN, visual; FPN, frontoparietal; SN, salience; CO, cingulo-opercular; VAN, ventral attention; SMm, somatomotor mouth; SMh, somatomotor hand; RSpN, retrosplenial; FTP, frontotemporal; DMN, default mode; AN, auditory; CPN, cinguloparietal; SubC, subcortex; Cer, cerebellar.

Close modal

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; QM1 = 0.137, p = 0; QM2 = 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−3, t(17) = 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).

Figure 3.

Relationships between load effect on participation, pupil load effect, and drift rate. (A) Average participation coefficient (PC) for each load, for the red module, during tracking. Each color represents the corresponding tracking load (from 2 to 5). Gray lines correspond to each subject. (B–C) A regression parameter (βPC) was calculated for each subject and then correlated to βpupil (B; r = 0.62; p = 0.006) and drift rate (C; r = 0.66; p = 2.4 × 10−3). Each circle corresponds to the mean value per subject.

Figure 3.

Relationships between load effect on participation, pupil load effect, and drift rate. (A) Average participation coefficient (PC) for each load, for the red module, during tracking. Each color represents the corresponding tracking load (from 2 to 5). Gray lines correspond to each subject. (B–C) A regression parameter (βPC) was calculated for each subject and then correlated to βpupil (B; r = 0.62; p = 0.006) and drift rate (C; r = 0.66; p = 2.4 × 10−3). Each circle corresponds to the mean value per subject.

Close modal

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 ± 3 × 10−3, t(70) = 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 (βQL, 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 βQL (r = 0.51, p = 0.034), and the partial correlation between drift rate and βQL 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.

### Network Mesoscale Integration and Adrenergic Receptor Density

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(194) = 3.82, p = 2 × 10−4, pSA = 0.02) and the differences between green and yellow submodules (t(177) = −4.47, p = 1.3 × 10−5, pSA = 0.004), while the other differences did not survive the spatial autocorrelation test (green-red: t(152) = 0.47; p = 0.635, pSA = 0.590; yellow-red: t(156) = −3.02, p = 0.003, pSA = 0.121; green-blue: t(173) = −0.68, p = 0.496, pSA = 0.324; red-blue: t(135) = −1.30, p = 0.195, pSA = 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.

Figure 4.

Receptor density analysis. (A) Spatial maps of α2a density (left), BOLD parametric effect (middle), and participation coefficient parametric effect (right). The ∼ symbol represents the linear model tested in the analysis. (B) Scatterplot depicting the relationship between βPupil and the random effect of α2a (RE α2a; r = 0.54, p = 0.02). (C) Scatterplot depicting the relationship between the random effect of α2a and drift rate (r = 0.70, p = 0.001); the colors of the dots represent the pSA value from the linear effect of α2a on βBOLD within each subject, and the marked circles correspond to subjects with pSA < 0.05. (D) Pearson correlation of the group average BOLD parametric effect (βBOLD) and participation coefficient (βPC; r = 0.26, p = 7 × 10−7). Colors represent each module assignment as in Figure 2.

Figure 4.

Receptor density analysis. (A) Spatial maps of α2a density (left), BOLD parametric effect (middle), and participation coefficient parametric effect (right). The ∼ symbol represents the linear model tested in the analysis. (B) Scatterplot depicting the relationship between βPupil and the random effect of α2a (RE α2a; r = 0.54, p = 0.02). (C) Scatterplot depicting the relationship between the random effect of α2a and drift rate (r = 0.70, p = 0.001); the colors of the dots represent the pSA value from the linear effect of α2a on βBOLD within each subject, and the marked circles correspond to subjects with pSA < 0.05. (D) Pearson correlation of the group average BOLD parametric effect (βBOLD) and participation coefficient (βPC; r = 0.26, p = 7 × 10−7). Colors represent each module assignment as in Figure 2.

Close modal

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(5992) = 2.29; p = 0.022; pSA = 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, pSA = 0.036; Figure 4B) and mean drift rate (Pearson’s r = 0.70, p = 0.001, pSA = 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 4BC, we observed a dependency between the pSA value, βpupil, and drift rate, in which the respective within-subject effects that survived the spatial autocorrelation correction are shown (pSA < 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(5992) = −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(5992) = 3.96; p = 7.55 × 10−5; 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(17) = −12.99, p = 2.93 × 10−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.

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.

### Participants

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).

### Behavior and EZ-Diffusion Model

The EZ-diffusion model was used to interpret the performance measures from the task (Ratcliff & Rouder, 1998; Wagenmakers et al., 2007). This model considers the mean RT of correct trials, the standard deviation of the reaction time (SD-RT) across correct trials, and mean accuracy across the task, and computes from these a value for drift rate (v, Equation 1), boundary separation (a, Equation 2) , and non-decision time (Equation 3)—the three main parameters for the drift diffusion model (Ratcliff et al., 2016; Ratcliff & Rouder, 1998).
$v=signP−12·0.1·logP1−P·P2·logP1−P−P·logP1−P+P−12VRT,4$
(1)
$a=0.01·logP1−Pv,$
(2)
$Ter=MRT−a2×v×1−e−100·v·a1+e−100·v·a,$
(3)
in which P is the average performance (range between 0 and 1); sign is an operator that will be −1 if P < 0.5 or +1 if P > 0.5; VRT is the standard deviation of reaction time (in seconds); and MRT is the mean reaction time (in seconds).

### Pupillometry

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.

### MRI Data

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).

### fMRI Data Preprocessing

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.

### Brain Parcellation

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.

### Time-Resolved Functional Connectivity and Network Analysis

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.

### Community Structure

The Louvain modularity algorithm from the Brain Connectivity Toolbox (Rubinov & Sporns, 2010) was used in combination with the JC to estimate both time-averaged and time-resolved community structure. The Louvain algorithm iteratively maximizes the modularity statistic, Q, for different community assignments until the maximum possible score of Q has been obtained (Equation 4).
$QT=1v+∑ijwij+−eij+δMiMj−1v++v−∑ijwij−−eij−δMiMj.$
(4)
Equation 4: Louvain modularity algorithm, where v is the total weight of the network (sum of all negative and positive connections); wij is the weighted and signed connection between regions i and j; eij is the strength of a connection divided by the total weight of the network; and δMiMj is set to 1 when regions are in the same community and 0 otherwise. The + and − superscripts denote all positive and negative connections, respectively.

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 (QT) within each temporal window for each participant in the study.

### Regional Integration

Based on the group consensus community assignments, we estimated between-module connectivity using the participation coefficient, BT, which quantifies the extent to which a region connects across all modules (i.e., between-module strength; Equation 5). In our experiment, we used two separate community assignments, one for each of the modularity levels. In this manner we measure (a) how the first hierarchical-level (i.e., large-scale) topology changed during tracking across the complete brain; and (b) how the topology of the submodules changed across the task. These values were calculated in each time point using the time-resolved adjacency matrix across each load condition.
$BiT=1−∑s=1nMκisTκiT2.$
(5)
Equation 5: Participation coefficient BiT, where κisT is the strength of the positive connections of region i to regions in module s at time T, and κiT is the sum of strengths of all positive connections of region i at time T. The participation coefficient of a region is therefore close to 1 if its connections are uniformly distributed among all the modules and 0 if all of its links are within its own module.

### Neurotransmitter Receptor Mapping

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.

### The Relationship Between Sympathetic Tone and Attentional Processing

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 the statistical effect of pupil diameter on accuracy, we performed a logistic linear mixed-effects model. We used the mean pupil diameter of the significant time period (Figure 1A) of the high load trials (Loads 4 and 5) and the accuracy (i.e., correct or incorrect) as the predictor variable of each trial, grouping by subject as the random effect. The statistical model is described in the following equation:
$Accuracy∼Pupil+1+Pupil+1Subject.$
(6)

### Network Integration Increases as a Function of Attentional Load

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).

We analyzed the statistical effect of pupil diameter on the participation coefficient both within and between subjects by performing a linear mixed model using the time-varying PC of the red submodule (Figure 3A) of each load as a dependent variable (N = 72), and the respective pupil diameter as a regressor, with grouping by subject. The statistical model is described in the following equation:
$PC∼Pupil+1+Pupil+1Subject.$
(7)

### Network Mesoscale Integration and Adrenergic Receptor Density

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 measure the statistical difference in the receptor density between submodules by a two-tailed t test between each pair of modules. The same procedure was performed using the surrogate maps to generate a null distribution of t statistics. To evaluate the effect of the density of each adrenergic receptor on the neural activity in the attentional task, we built a linear mixed model aimed at predicting regional differences in BOLD activity and participation coefficient. We created a model using the receptor density atlas of α2a receptor to predict parametric BOLD activity (i.e., linear increase of BOLD activity with task load) during tracking (Equation 8). To evaluate the relationship between BOLD activity, adrenergic receptor expression, and changes in participation coefficient as a function of attentional load, we tested two models: one using the adrenergic receptor density as independent factor (Equation 9), and another using the parametric BOLD effect as an independent factor (Equation 10). Additionally, we assessed the across-subject variability using the subjects’ ID as a grouping variable in order to evaluate the random effects on the independent factor. We corrected the spatial autocorrelation by running the same model using 5,000 surrogate maps. Then we used the fixed-effect null distribution to calculate the pSA (i.e., the probability of finding the fixed effect within the 95th percentile of the null distribution). The deterministic part of the model is expressed in the following equations (Wilkinson & Rogers, 1973):
$βBOLD∼α2a+1+α2a+1Subject,$
(8)
$βPC∼α2a+1+α2a+1Subject,$
(9)
$βPC∼BOLD+1+BOLD+1Subject,$
(10)
where PC is the parametric effect of mesoscale participation coefficient (i.e., βPC), BOLD is the parametric effect of load on BOLD activity during tracking for each region, and α2a are the regional densities of the respective adrenergic receptor atlas. We then correlated the random-effects parameters to pupil diameter responses and behavior and then compared these with the Pearson’s correlation of the null distribution using the random effect of the surrogate maps. Finally, we performed a linear model within each subject with α2a as a regressor and βBOLD as dependent variable. Again, the statistical effect (i.e., β value) was compared against the null distribution when performing the regression using the surrogate maps (Figure 4BC).

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.

•

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.

•

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.

•

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

Akiki
,
T. J.
, &
Abdallah
,
C. G.
(
2019
).
Determining the hierarchical architecture of the human brain using subject-level clustering of functional networks
.
Scientific Reports
,
9
(
1
),
1
15
. ,
[PubMed]
Alnæs
,
D.
,
Kaufmann
,
T.
,
Richard
,
G.
,
Duff
,
E. P.
,
Sneve
,
M. H.
,
,
T.
,
Nordvik
,
J. E.
,
Andreassen
,
O. A.
,
Smith
,
S. M.
, &
Westlye
,
L. T.
(
2015
).
Attentional load modulates large-scale functional brain connectivity beyond the core attention networks
.
NeuroImage
,
109
,
260
272
. ,
[PubMed]
Alnæs
,
D.
,
Sneve
,
M. H.
,
Espeseth
,
T.
,
,
T.
,
van de Pavert
,
S. H. P.
, &
Laeng
,
B.
(
2014
).
Pupil size signals mental effort deployed during multiple object tracking and predicts brain activity in the dorsal attention network and the locus coeruleus
.
Journal of Vision
,
14
(
2014
),
1
20
. ,
[PubMed]
Arnsten
,
A. F. T.
, &
Haven
,
N.
(
2013
).
The neurobiology of thought: The groundbreaking discoveries of Patricia Goldman-Rakic 1937–2003
.
Cerebral Cortex
,
23
(
10
),
2269
2281
. ,
[PubMed]
Arnsten
,
A. F. T.
,
Wang
,
M. J.
, &
Paspalas
,
C. D.
(
2012
).
Neuromodulation of thought: Flexibilities and vulnerabilities in prefrontal cortical network synapses
.
Neuron
,
76
(
1
),
223
239
. ,
[PubMed]
Aston-Jones
,
G.
, &
Cohen
,
J. D.
(
2005
).
An integrative theory of locus coeruleus-norepinephrine function: Adaptive gain and optimal performance
.
Annual Review of Neuroscience
,
28
(
1
),
403
450
. ,
[PubMed]
Aston-Jones
,
G.
, &
Waterhouse
,
B.
(
2016
).
Locus coeruleus: From global projection system to adaptive regulation of behavior
.
Brain Research
,
1645
,
75
78
. ,
[PubMed]
Avery
,
M. C.
, &
Krichmar
,
J. L.
(
2017
).
Neuromodulatory systems and their interactions: A review of models, theories, and experiments
.
Frontiers in Neural Circuits
,
11
(
December
),
1
18
. ,
[PubMed]
Bassett
,
D. S.
,
Greenfield
,
D. L.
,
Meyer-Lindenberg
,
A.
,
Weinberger
,
D. R.
,
Moore
,
S. W.
, &
Bullmore
,
E. T.
(
2010
).
Efficient physical embedding of topologically complex information processing networks in brains and computer circuits
.
PLoS Computational Biology
,
6
(
4
),
e1000748
. ,
[PubMed]
,
Y.
,
Restom
,
K.
,
Liau
,
J.
, &
Liu
,
T. T.
(
2007
).
A component based noise correction method (CompCor) for BOLD and perfusion based fMRI
.
NeuroImage
,
37
(
1
),
90
101
. ,
[PubMed]
Bekar
,
L. K.
,
Wei
,
H. S.
, &
Nedergaard
,
M.
(
2012
).
The locus coeruleus-norepinephrine network optimizes coupling of cerebral blood volume with oxygen demand
.
Journal of Cerebral Blood Flow and Metabolism
,
32
(
12
),
2135
2145
. ,
[PubMed]
Benjamini
,
Y.
, &
Yekutieli
,
D.
(
2001
).
The control of the false discovery rate in multiple testing under dependency
.
Annals of Statistics
,
29
(
4
),
1165
1188
.
Bouret
,
S.
, &
Sara
,
S. J.
(
2005
).
Network reset: A simplified overarching theory of locus coeruleus noradrenaline function
.
Trends in Neurosciences
,
28
(
11
),
574
582
. ,
[PubMed]
Briand
,
L. A.
,
Gritton
,
H.
,
Howe
,
W. M.
,
Young
,
D. A.
, &
Sarter
,
M.
(
2007
).
Modulators in concert for cognition: Modulator interactions in the prefrontal cortex
.
Progress in Neurobiology
,
83
(
2
),
69
91
. ,
[PubMed]
Burt
,
J. B.
,
Helmer
,
M.
,
Shinn
,
M.
,
Anticevic
,
A.
, &
Murray
,
J. D.
(
2020
).
Generative modeling of brain maps with spatial autocorrelation
.
NeuroImage
,
220
(
February
),
117038
. ,
[PubMed]
Campos-Arteaga
,
G.
,
Forcato
,
C.
,
Wainstein
,
G.
,
Lagos
,
R.
,
Palacios-García
,
I.
,
Artigas
,
C.
,
Morales
,
R.
,
Pedreira
,
M. E. E.
, &
Rodríguez
,
E.
(
2020
).
Differential neurophysiological correlates of retrieval of consolidated and reconsolidated memories in humans: An ERP and pupillometry study
.
Neurobiology of Learning and Memory
,
107279
. ,
[PubMed]
Cazettes
,
F.
,
Reato
,
D.
,
Morais
,
J. P.
,
Renart
,
A.
, &
Mainen
,
Z. F.
(
2020
).
Phasic activation of dorsal raphe serotonergic neurons increases pupil size
.
Current Biology
,
192
197
. ,
[PubMed]
de Gee
,
J. W.
,
Colizoli
,
O.
,
Kloosterman
,
N. A.
,
Knapen
,
T.
,
Nieuwenhuis
,
S.
, &
Donner
,
T. H.
(
2017
).
Dynamic modulation of decision biases by brainstem arousal systems
.
eLife
,
6
,
e23232
. ,
[PubMed]
de Gee
,
J. W.
,
Knapen
,
T.
, &
Donner
,
T. H.
(
2014
).
Decision-related pupil dilation reflects upcoming choice and individual bias
.
Proceedings of the National Academy of Sciences
,
111
(
5
),
E618
E625
. ,
[PubMed]
Diedrichsen
,
J.
,
Balsters
,
J. H.
,
Flavell
,
J.
,
Cussans
,
E.
, &
Ramnani
,
N.
(
2009
).
A probabilistic MR atlas of the human cerebellum
.
NeuroImage
,
46
(
1
),
39
46
. ,
[PubMed]
Donner
,
T.
,
Kettermann
,
A.
,
Diesch
,
E.
,
Ostendorf
,
F.
,
Villringer
,
A.
, &
Brandt
,
S. A.
(
2000
).
Involvement of the human frontal eye field and multiple parietal areas in covert visual selection during conjunction search
.
European Journal of Neuroscience
,
12
(
9
),
3407
3414
. ,
[PubMed]
Douglas
,
R. J.
, &
Martin
,
K. A. C.
(
2004
).
Neuronal circuits of the neocortex
.
Annual Review of Neuroscience
,
27
(
1
),
419
451
. ,
[PubMed]
Fornito
,
A.
,
Arnatkevičiūtė
,
A.
, &
Fulcher
,
B. D.
(
2019
).
Bridging the gap between connectome and transcriptome
.
Trends in Cognitive Sciences
,
23
(
1
),
34
50
. ,
[PubMed]
Gilzenrat
,
M. S.
,
Nieuwenhuis
,
S.
,
Jepma
,
M.
, &
Cohen
,
J. D.
(
2010
).
Pupil diameter tracks changes in control state predicted by the adaptive gain theory of locus coeruleus function
.
Cognitive, Affective and Behavioral Neuroscience
,
10
(
2
),
252
269
. ,
[PubMed]
Gordon
,
E. M.
,
Laumann
,
T. O.
,
,
B.
,
Huckins
,
J. F.
,
Kelley
,
W. M.
, &
Petersen
,
S. E.
(
2016
).
Generation and evaluation of a cortical area parcellation from resting-state correlations
.
Cerebral Cortex
,
26
(
1
),
288
303
. ,
[PubMed]
Gryglewski
,
G.
,
Seiger
,
R.
,
James
,
G. M.
,
Godbersen
,
G. M.
,
Komorowski
,
A.
,
Unterholzner
,
J.
,
Michenthaler
,
P.
,
Hahn
,
A.
,
,
W.
,
Mitterhauser
,
M.
,
Kasper
,
S.
, &
Lanzenberger
,
R.
(
2018
).
Spatial analysis and high resolution mapping of the human whole-brain transcriptome for integrative analysis in neuroimaging
.
NeuroImage
,
176
(
January
),
259
267
. ,
[PubMed]
Gu
,
S.
,
Pasqualetti
,
F.
,
Cieslak
,
M.
,
Telesford
,
Q. K.
,
Yu
,
A. B.
,
Kahn
,
A. E.
,
Medaglia
,
J. D.
,
Vettel
,
J. M.
,
Miller
,
M. B.
,
Grafton
,
S. T.
, &
Bassett
,
D. S.
(
2015
).
Controllability of structural brain networks
.
Nature Communications
,
6
,
1
10
. ,
[PubMed]
Hawrylycz
,
M. J.
,
Lein
,
E. S.
,
Guillozet-bongaarts
,
A. L.
,
Shen
,
E. H.
,
Ng
,
L.
,
Miller
,
J. A.
,
van de Lagemaat
,
L. N.
,
Smith
,
K. A.
,
Ebbert
,
A.
,
Riley
,
Z. L.
,
Abajian
,
C.
,
Beckmann
,
C. F.
,
Bernard
,
A.
,
Bertagnolli
,
D.
,
Boe
,
A. F.
,
Cartagena
,
P. M.
,
Chakravarty
,
M. M.
,
Chapin
,
M.
,
Chong
,
J.
, …
Zielke
,
H. R.
(
2012
).
An anatomically comprehensive atlas of the adult human brain transcriptome
.
Nature
,
489
(
7416
),
391
399
. ,
[PubMed]
Joshi
,
A.
,
,
V.
,
Vemana
,
V.
,
McGinnity
,
T. M.
,
,
G.
, &
Wong-Lin
,
K. F.
(
2017
).
An integrated modelling framework for neural circuits with multiple neuromodulators
.
Journal of the Royal Society Interface
,
14
(
126
). ,
[PubMed]
Joshi
,
S.
, &
Gold
,
J. I.
(
2020
).
Pupil size as a window on neural substrates of cognition
.
Trends in Cognitive Sciences
,
24
(
6
),
466
480
. ,
[PubMed]
Joshi
,
S.
,
Li
,
Y.
,
Kalwani
,
R. M.
, &
Gold
,
J. I.
(
2016
).
Relationships between pupil diameter and neuronal activity in the locus coeruleus, colliculi, and cingulate cortex
.
Neuron
,
89
(
1
),
221
234
. ,
[PubMed]
Kahneman
,
D.
, &
Beatty
,
J.
(
1966
).
Pupil diameter and load on memory
.
Science
,
154
(
3756
),
1583
1585
. ,
[PubMed]
Kaymak
,
H.
,
Fricke
,
A.
,
Mauritz
,
Y.
,
Löwinger
,
A.
,
Klabe
,
K.
,
Breyer
,
D.
,
Lagenbucher
,
A.
,
Seitz
,
B.
, &
Schaeffel
,
F.
(
2018
).
Short-term effects of low-concentration atropine eye drops on pupil size and accommodation in young adult subjects
.
Graefe’s Archive for Clinical and Experimental Ophthalmology
,
2211
2217
. ,
[PubMed]
Lecrux
,
C.
, &
Hamel
,
E.
(
2016
).
Neuronal networks and mediators of cortical neurovascular coupling responses in normal and altered brain states
.
Philosophical Transactions of the Royal Society B: Biological Sciences
,
371
(
1705
),
20150350
. ,
[PubMed]
Liu
,
Y.
,
Rodenkirch
,
C.
,
Moskowitz
,
N.
,
Schriver
,
B.
, &
Wang
,
Q.
(
2017
).
Dynamic lateralization of pupil dilation evoked by locus coeruleus activation results from sympathetic, not parasympathetic, contributions
.
Cell Reports
,
20
(
13
),
3099
3112
. ,
[PubMed]
Logothetis
,
N. K.
(
2003
).
The underpinnings of the BOLD functional magnetic resonance imaging signal
.
Journal of Neuroscience
,
23
(
10
),
3963
3971
. ,
[PubMed]
Markello
,
R. D.
, &
Misic
,
B.
(
2021
).
Comparing spatial null models for brain maps
.
NeuroImage
,
236
(
April
),
118052
. ,
[PubMed]
Mather
,
M.
,
Clewett
,
D.
,
Sakaki
,
M.
, &
Harley
,
C. W.
(
2016
).
Norepinephrine ignites local hotspots of neuronal excitation: How arousal amplifies selectivity in perception and memory
.
Behavioral and Brain Sciences
,
39
. ,
[PubMed]
McGinley
,
M. J.
,
Vinck
,
M.
,
Reimer
,
J.
,
Batista-Brito
,
R.
,
Zagha
,
E.
,
,
C. R.
,
Tolias
,
A. S.
,
Cardin
,
J. A.
, &
McCormick
,
D. A.
(
2015
).
Waking state: Rapid variations modulate neural and behavioral responses
.
Neuron
,
87
(
6
),
1143
1161
. ,
[PubMed]
Meunier
,
D.
,
Lambiotte
,
R.
, &
Bullmore
,
E. T.
(
2010
).
Modular and hierarchically modular organization of brain networks
.
Frontiers in Neuroscience
,
4
(
200
),
1
11
. ,
[PubMed]
Meunier
,
D.
,
Lambiotte
,
R.
,
Fornito
,
A.
,
Ersche
,
K. D.
, &
Bullmore
,
E. T.
(
2009
).
Hierarchical modularity in human brain functional networks
.
Frontiers in Neuroinformatics
,
3
(
37
),
1
12
. ,
[PubMed]
Mohr
,
H.
,
Wolfensteller
,
U.
,
Betzel
,
R. F.
,
Mišić
,
B.
,
Sporns
,
O.
,
Richiardi
,
J.
,
Ruge
,
H.
,
Mohr
,
H.
,
Wolfensteller
,
U.
,
Betzel
,
R. F.
, &
Mis
,
B.
(
2016
).
Integration and segregation of large-scale brain networks during short-term task automatization
.
Nature Communications
,
7
(
1
),
1
12
. ,
[PubMed]
Mulder
,
G.
(
2012
).
The concept and measurement of mental effort
.
Energetics and Human Information Processing
,
175
198
.
Munn
,
B.
,
Müller
,
E. J.
,
Wainstein
,
G.
, &
Shine
,
J. M.
(
2021
).
The ascending arousal system shapes low-dimensional brain dynamics to mediate awareness of changes in intrinsic cognitive states
.
bioRxiv
.
Murphy
,
P. R.
,
Boonstra
,
E.
, &
Nieuwenhuis
,
S.
(
2016
).
Global gain modulation generates time-dependent urgency during perceptual choice in humans
.
Nature Communications
,
7
(
1
),
1
15
. ,
[PubMed]
Murphy
,
P. R.
,
Vandekerckhove
,
J.
, &
Nieuwenhuis
,
S.
(
2014
).
Pupil-linked arousal determines variability in perceptual decision making
.
PLoS Computational Biology
,
10
(
9
). ,
[PubMed]
Palomero-Gallagher
,
N.
, &
Zilles
,
K.
(
2019
).
Cortical layers: Cyto-, myelo-, receptor- and synaptic architecture in human cortical areas
.
NeuroImage
,
197
,
716
741
. ,
[PubMed]
Parvizi
,
J.
, &
Damasio
,
A.
(
2001
).
Consciousness and the brainstem
.
Cognition
,
79
(
1–2
),
135
160
. ,
[PubMed]
Perez
,
D. M.
(
2020
).
A1-adrenergic receptors in neurotransmission, synaptic plasticity, and cognition
.
Frontiers in Pharmacology
,
11
(
September
),
1
22
. ,
[PubMed]
Power
,
J. D.
,
Mitra
,
A.
,
Laumann
,
T. O.
,
Snyder
,
A. Z.
,
Schlaggar
,
B. L.
, &
Petersen
,
S. E.
(
2014
).
Methods to detect, characterize, and remove motion artifact in resting state fMRI
.
NeuroImage
,
84
,
320
341
. ,
[PubMed]
Ratcliff
,
R.
, &
Rouder
,
J.
(
1998
).
Modeling response times for two-choice decisions
.
Psychological Science
,
9
(
5
),
347
356
.
Ratcliff
,
R.
,
Smith
,
P. L.
,
Brown
,
S. D.
, &
Mckoon
,
G.
(
2016
).
Diffusion decision model: Current issues and history
.
Trends in Cognitive Sciences
,
20
(
4
),
260
281
. ,
[PubMed]
Ratcliff
,
R.
,
Thompson
,
C. A.
, &
McKoon
,
G.
(
2015
).
Modeling individual differences in response time and accuracy in numeracy
.
Cognition
,
137
,
115
136
. ,
[PubMed]
Reimer
,
J.
,
Froudarakis
,
E.
,
,
C. R.
,
Yatsenko
,
D.
,
Denfield
,
G. H.
, &
Tolias
,
A. S.
(
2014
).
Pupil fluctuations track fast switching of cortical states during quiet wakefulness
.
Neuron
,
84
(
2
),
355
362
. ,
[PubMed]
Reimer
,
J.
,
McGinley
,
M. J.
,
Liu
,
Y.
,
Rodenkirch
,
C.
,
Wang
,
Q.
,
McCormick
,
D. A.
, &
Tolias
,
A. S.
(
2016
).
Pupil fluctuations track rapid changes in adrenergic and cholinergic activity in cortex
.
Nature Communications
,
7
(
May
),
1
7
. ,
[PubMed]
Richiardi
,
J.
,
Altmann
,
A.
,
Milazzo
,
A.-C.
,
Chang
,
C.
,
Chakravarty
,
M. M.
,
Banaschewski
,
T.
,
Barker
,
G. J.
,
Bokde
,
A. L. W.
,
Bromberg
,
U.
,
Büchel
,
C.
,
Conrod
,
P.
,
Fauth-Bühler
,
M.
,
Flor
,
H.
,
Frouin
,
V.
,
Gallinat
,
J.
,
Garavan
,
H.
,
Gowland
,
P.
,
Heinz
,
A.
,
Lemaître
,
H.
, …
IMAGEN Consortium
. (
2015
).
Correlated gene expression supports synchronous activity in brain networks
.
Science
,
348
(
6240
),
1241
1244
. ,
[PubMed]
Robbins
,
T. W.
, &
Arnsten
,
A. F. T.
(
2009
).
The neuropsychopharmacology of fronto-executive function: Monoaminergic modulation
.
Annual Review of Neuroscience
,
32
,
267
289
. ,
[PubMed]
Rojas-Líbano
,
D.
,
Wainstein
,
G.
,
Carrasco
,
X.
,
Aboitiz
,
F.
,
Crossley
,
N.
, &
Ossandón
,
T.
(
2019
).
A pupil size, eye-tracking and neuropsychological dataset from ADHD children during a cognitive task
.
Scientific Data
,
6
(
1
),
25
. ,
[PubMed]
Rubinov
,
M.
, &
Sporns
,
O.
(
2010
).
Complex network measures of brain connectivity: Uses and interpretations
.
NeuroImage
,
52
(
3
),
1059
1069
. ,
[PubMed]
Samuels
,
E.
, &
,
E.
(
2008
).
Functional neuroanatomy of the noradrenergic locus coeruleus: Its roles in the regulation of arousal and autonomic function part II: Physiological and pharmacological manipulations and pathological alterations of locus coeruleus activity in humans
.
Current Neuropharmacology
,
6
(
3
),
254
285
. ,
[PubMed]
Santana
,
N.
, &
Artigas
,
F.
(
2017
).
Laminar and cellular distribution of monoamine receptors in rat medial prefrontal cortex
.
Frontiers in Neuroanatomy
,
11
(
87
),
1
13
. ,
[PubMed]
Sara
,
S. J.
(
2009
).
The locus coeruleus and noradrenergic modulation of cognition
.
Nature Reviews Neuroscience
,
10
(
3
),
211
223
. ,
[PubMed]
Shenhav
,
A.
,
Musslick
,
S.
,
Lieder
,
F.
,
Kool
,
W.
,
Griffiths
,
T. L.
,
Cohen
,
J. D.
, &
Botvinick
,
M. M.
(
2017
).
Toward a rational and mechanistic account of mental effort
.
Annual Review of Neuroscience
,
40
(
1
),
99
124
. ,
[PubMed]
Shine
,
J. M.
(
2019
).
Neuromodulatory influences on integration and segregation in the brain
.
Trends in Cognitive Sciences
,
23
(
7
),
572
583
. ,
[PubMed]
Shine
,
J. M.
(
2021
).
The thalamus integrates the macrosystems of the brain to facilitate complex, adaptive brain network dynamics
.
Progress in Neurobiology
,
199
,
101951
. ,
[PubMed]
Shine
,
J. M.
,
Bissett
,
P. G.
,
Bell
,
P. T.
,
Koyejo
,
O.
,
Balsters
,
J. H.
,
Gorgolewski
,
K. J.
,
Moodie
,
C. A.
, &
Poldrack
,
R. A.
(
2016
).
The dynamics of functional brain networks: Integrated network states during cognitive task performance
.
Neuron
,
92
(
2
),
544
554
. ,
[PubMed]
Shine
,
J. M.
,
Breakspear
,
M.
,
Bell
,
P. T.
,
Ehgoetz Martens
,
K.
,
Shine
,
R.
,
Koyejo
,
O.
,
Sporns
,
O.
, &
Poldrack
,
R. A.
(
2019
).
Human cognition involves the dynamic integration of neural activity and neuromodulatory systems
.
Nature Neuroscience
,
22
(
2
),
289
296
. ,
[PubMed]
Shine
,
J. M.
,
Hearne
,
L. J.
,
Breakspear
,
M.
,
Hwang
,
K.
,
Müller
,
E. J.
,
Sporns
,
O.
,
Poldrack
,
R. A.
,
Mattingley
,
J. B.
, &
Cocchi
,
L.
(
2019
).
The low-dimensional neural architecture of cognitive complexity is related to activity in medial thalamic nuclei
.
Neuron
,
104
(
5
),
849
855
. ,
[PubMed]
Shine
,
J. M.
,
Müller
,
E. J.
,
Munn
,
B.
,
Cabral
,
J.
,
Moran
,
R. J.
, &
Breakspear
,
M.
(
2021
).
Computational models link cellular mechanisms of neuromodulation to large-scale neural dynamics
.
Nature Neuroscience
,
24
(
6
),
765
776
. ,
[PubMed]
Shine
,
J. M.
,
van den Brink
,
R. L.
,
Hernaus
,
D.
,
Nieuwenhuis
,
S.
, &
Poldrack
,
R. A.
(
2018
).
Catecholaminergic manipulation alters dynamic network topology across cognitive states
.
Network Neuroscience
,
2
(
3
),
381
396
. ,
[PubMed]
Smiley
,
J. F.
,
Subramanian
,
M.
, &
Mesulam
,
M. M.
(
1999
).
Monoaminergic-cholinergic interactions in the primate basal forebrain
.
Neuroscience
,
93
(
3
),
817
829
. ,
[PubMed]
Thompson
,
W. H.
,
Richter
,
C. G.
,
Plavén-Sigray
,
P.
, &
Fransson
,
P.
(
2018
).
Simulations to benchmark time-varying connectivity methods for fMRI
.
PLoS Computational Biology
,
14
(
5
),
1
23
. ,
[PubMed]
van den Brink
,
R. L.
,
Pfeffer
,
T.
,
Warren
,
C. M.
,
Murphy
,
P. R.
,
Tona
,
K.-D.
,
van der Wee
,
N. J. A.
,
Giltay
,
E.
,
van Noorden
,
M. S.
,
Rombouts
,
S. A. R. B.
, &
Donner
,
T. H.
(
2016
).
Catecholaminergic neuromodulation shapes intrinsic MRI functional connectivity in the human brain
.
Journal of Neuroscience
,
36
(
30
),
7865
7876
. ,
[PubMed]
Varazzani
,
C.
,
San-Galli
,
A.
,
Gilardeau
,
S.
, &
Bouret
,
S.
(
2015
).
Noradrenaline and dopamine neurons in the reward/effort trade-off: A direct electrophysiological comparison in behaving monkeys
.
Journal of Neuroscience
,
35
(
20
),
7866
7877
. ,
[PubMed]
Varela
,
F.
,
Lachaux
,
J.-P.
,
Rodriguez
,
E.
, &
Martinerie
,
J.
(
2001
).
The brainweb: Phase synchronization and large-scale integration
.
Nature Reviews Neuroscience
,
2
(
4
),
229
239
. ,
[PubMed]
Vos de Wael
,
R.
,
Benkarim
,
O.
,
Paquola
,
C.
,
Lariviere
,
S.
,
Royer
,
J.
,
Tavakol
,
S.
,
Xu
,
T.
,
Hong
,
S. J.
,
Langs
,
G.
,
Valk
,
S.
,
Misic
,
B.
,
Milham
,
M.
,
Margulies
,
D.
,
Smallwood
,
J.
, &
Bernhardt
,
B. C.
(
2020
).
BrainSpace: A toolbox for the analysis of macroscale gradients in neuroimaging and connectomics datasets
.
Communications Biology
,
3
(
1
). ,
[PubMed]
Wagenmakers
,
E. J.
,
Van Der Maas
,
H. L. J.
, &
Grasman
,
R. P. P. P.
(
2007
).
An EZ-diffusion model for response time and accuracy
.
Psychonomic Bulletin and Review
,
14
(
1
),
3
22
. ,
[PubMed]
Wainstein
,
G
. (
2021
).
MOT, GitHub
, https://github.com/gabwainstein/MOT
Wainstein
,
G.
,
Rojas-Líbano
,
D.
,
Alnæs
,
D.
,
Kolskår
,
K.
,
,
T.
,
Laeng
,
B.
,
Ossandon
,
T.
,
Crossley
,
N.
,
Matar
,
E.
, &
Shine
,
J. M.
(
2020
).
MOT_data.mat (Version 1), Figshare
.
Wainstein
,
G.
,
Rojas-Líbano
,
D.
,
Crossley
,
N. A.
,
Carrasco
,
X.
,
Aboitiz
,
F.
, &
Ossandón
,
T.
(
2017
).
Pupil size tracks attentional performance in attention-deficit/hyperactivity disorder
.
Scientific Reports
,
7
(
1
),
1
9
. ,
[PubMed]
Wang
,
M.
,
Ramos
,
B. P.
,
Paspalas
,
C. D.
,
Shu
,
Y.
,
Simen
,
A.
,
Duque
,
A.
,
Vijayraghavan
,
S.
,
Brennan
,
A.
,
Dudley
,
A.
, &
Nou
,
E.
(
2007
).
A2A-adrenoceptors strengthen working memory networks by inhibiting cAMP-HCN channel signaling in prefrontal cortex
.
Cell
,
129
(
2
),
397
410
. ,
[PubMed]
Wang
,
X. J.
(
2020
).
Macroscopic gradients of synaptic excitation and inhibition in the neocortex
.
Nature Reviews Neuroscience
,
21
(
3
),
169
178
. ,
[PubMed]
Wilkinson
,
G. N.
, &
Rogers
,
C. E.
(
1973
).
Symbolic description of factorial models for analysis of variance
.
Journal of Applied Statistics
,
22
(
3
),
392
399
.
Yüzgeç
,
Ö.
,
Prsa
,
M.
,
Zimmermann
,
R.
, &
Huber
,
D.
(
2018
).
Pupil size coupling to cortical states protects the stability of deep sleep via parasympathetic modulation
.
Current Biology
,
28
(
3
),
392
400
. ,
[PubMed]
Zerbi
,
V.
,
Floriou-Servou
,
A.
,
Markicevic
,
M.
,
Vermeiren
,
Y.
,
Sturman
,
O.
,
Privitera
,
M.
,
Von Ziegler
,
L.
,
Ferrari
,
K. D.
,
Weber
,
B.
,
De Deyn
,
P. P.
,
Wenderoth
,
N.
, &
Bohacek
,
J.
(
2019
).
Rapid reconfiguration of the functional connectome after chemogenetic locus coeruleus activation
.
Neuron
,
103
(
4
),
702
718
.
Zilles
,
K.
, &
Palomero-Gallagher
,
N.
(
2017
).
Multiple transmitter receptors in regions and layers of the human cerebral cortex
.
Frontiers in Neuroanatomy
,
11
,
1
26
. ,
[PubMed]

## Author notes

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

Handling Editor: Christopher Honey

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. For a full description of the license, please visit https://creativecommons.org/licenses/by/4.0/legalcode.