Introduction

The human visual system faces a fundamental challenge: the visual inputs that bombard the retina are often fundamentally ambiguous. The inherent ambiguity of perceptual input is overcome by incorporating prior information about the causal structure of the world into sensory inferences. Consider the problem of watching a blurry item fly across the sky (“Is it a bird? Is it a plane?”) – this sensory information is often occluded and noisy, yet somehow our nervous system can adaptatively and dynamically refresh its activity to mediate the accurate perception (“It’s Superman!”). In cases of bistable perception, such as the Necker cube illusion, the inputs themselves offer two equally valid, yet mutually exclusive interpretations that observers periodically switch between. In both cases, perceptual awareness of the images (outside the moments of switching) is remarkably stable, suggesting that the nervous system can rapidly (and flexibly) identify the best ‘match’ between the visual data and stable (likely pre-learned) stimulus categories. The mapping between transitions in brain-states and transitions in perception remains an open question.

fMRI studies investigating perceptual updating and switches in bistable perception have typically identified a distributed set of regions within the cerebral cortex, including the anterior cingulate, anterior insula, dorsolateral prefrontal and superior parietal cortex^1,2. These cortical regions are presumed to play an attentional role with changes in attention driving switches in perceptual contents by selectively boosting activity within the relevant circuits^3,4. This interpretation is complemented by behavioural evidence showing that attention plays a prominent role in determining the contents of perception in bistable perception tasks where competition is not resolved at low-levels of the visual hierarchy^5,6. Despite these clear links, the ability to flexibly respond to ambiguous visual inputs according to changing task demands is a feature that is present across phylogeny⁷ and hence is present in a wide variety of animals that have poorly developed cerebral cortices⁸. This suggests that more highly conserved features of the nervous system may play a role in driving task-relevant switches in the contents of perception.

The ascending neuromodulatory system, and specifically the noradrenergic locus coeruleus (LC), is well-suited to mediate the large-scale, brain state switches required to flexibly alter perceptual contents^9,10. While the cell body of the LC is located in the brainstem, the nucleus sends projections throughout the central nervous system, wherein its axons release noradrenaline, which in turn modulate the excitability of targeted regions¹¹. In previous work, it has been argued that phasic release of noradrenaline from the LC acts as a “network reset” signal, which effectively disrupts ongoing processing, and hence allows animals to reconfigure their ongoing neural dynamics towards more salient (and hopefully, behaviourally-relevant) processes^12–14. This mechanism is of critical importance in ecological contexts in which an animal needs to be able to both focus on the current task in an exploitative mode (such as foraging), while being able to rapidly modify its internal, attentional and behavioural state when required (e.g., if resources are depleted, or in the presence of a predator). Further, in the context of bistable perception, pupil diameter (a non-specific and indirect readout of LC activity^9,15) is tightly linked to switches in the content of perception^16,17. Based on previous work^{12–14,18,19}, and the fact that the LC is known to be interconnected with many of the regions implicated in whole-brain imaging studies of perceptual uncertainty^20–22, we hypothesized that task-related perceptual switches should be intricately related to phasic bursts of LC activity, which then act as a ‘network reset’¹⁴, flattening the whole-brain energy landscape²³ and thus allowing cortical dynamics to evolve into a new state thereby changing the contents of perception.

To test this hypothesis, we analysed a cognitive task designed to investigate switches in perceptual categorisation²⁴. We first demonstrate that pupil diameter – a non-specific, indirect readout of phasic noradrenergic activity^15,20,25,26 – peaks at the point of the perceptual switch. We then trained a recurrent neural network (RNN) to perform an analogous change detection task, after which we altered the gain of activation function of the model – a demonstrated mechanism for the action of noradrenaline on the cerebral cortex^11,27,28 – and showed that switch times occurred earlier as a function of heightened gain. We then used principal component analysis (PCA) and energy landscape analysis²³ in order to translate fMRI data collected during the task into a low dimensional read-out of the network dynamics. We found that: i) increases in gain are associated with increases in the velocity of the low dimensional trajectory around the switch point, and ii) the energy landscape derived from the dynamics of the RNN flattens as a function of gain. We then re-analysed existing BOLD data from the same perceptual task^1,24 to confirm the prediction of brain state velocity increase and energy landscape flattening²³ at the point of perceptual transitions. Overall, our results confirm our hypothesis that phasic bursts of neuromodulatory activity act as a “network reset”^12,14, which via the downstream impact of neuromodulatory neurotransmitters such as noradrenaline on brain network dynamics, leads to a switch in perceptual categorisation.

Results

Evoked pupil dilations coincide with the resolution of perceptual ambiguity

To assess the role of the ascending arousal activity during task performance, we analysed a dataset of 35 participants who performed the ambiguous figures task while simultaneous pupil diameter was recorded (SR Research, 1000 Hz). Briefly, the task consisted of a set of continuously transforming images that transition from an initial object (e.g., a shark) into a second object (e.g., a plane), while preserving basic psychophysical attributes, such as luminance (Fig. 1A). Crucially, even though the task stimuli change incrementally and linearly, with maximal ambiguity at the mid-point of each trial (the peak of the dotted line curve in Fig. 1A), conscious awareness of a change in stimuli is known to ‘pop out’, often at different times on each trial²⁴. When these perceptual switches occurred, subjects were instructed to change the button they were pressing, thus indicating a change in perceptual interpretation while controlling for motor preparation and execution across stimuli. Participants viewed 20 unique sets of images, each of which morphed from a starting image into a second image through 15 equally spaced intermittent stages (Fig. 1A). For each participant, we identified the first and last time they viewed a sequence of images, as well as the three images leading up to and following an identified perceptual switch, irrespective of the categories associated with each specific object switch. The rest of the analyses in this manuscript are organised around this perceptual transition.

Figure 1

Given the known (yet admittedly non-specific) relationship between LC activity and the first-order dynamics of pupil diameter²⁹, we were able to use the diameter of the pupil as an indirect means for testing the hypothesis that ascending neuromodulatory tone is associated with perceptual switching. The linear nature of the morphing procedure meant that luminance levels (which could otherwise bias pupil diameter^9,15,20) were kept constant across all trials. Additional, behavioural output was controlled as subjects needed to press one button on each image (depending on the perceptual content). Mapping all blink-corrected, filtered and normalized trials over time, we noticed two clear and robust patterns. Firstly, we observed a clear increase in the phasic pupillary response coincident with the frame when participants noticed a change in perceptual content (maximum pupil diameter in Fig. 1C). These results can be seen clearly when the pupil diameter at the perceptual switch point is viewed in comparison to the remaining stimuli. Secondly, when analysing the times before and after the change point, there was a significant increase of the mean pupil response, beginning 3 images before the change point (mean β = -0.22; t₍₃₂₎ = -8.02, p = 2.3 × 10⁻¹⁹), before returning to baseline levels. Together, these results clearly couple pupil dilation to the point of perceptual switching.

In subsequent analyses, we sought to elucidate the relationship between ascending arousal, quantified by pupil diameter, and the temporal dynamics of perceptual shifts on a trial-by-trial basis. Given that the LC plays a pivotal role in modulating both sensory processing and the perceptual switches (Fig. 1B), we hypothesized that the speed of a perceptual switch should be dependent on LC activity. Specifically, we predicted that trials with earlier perceptual changes (faster perceptual switches) would be associated with an increase in pupil diameter (and v.v.). To test this prediction, we built a linear model wherein the peak pupil diameter observed during the perceptual switch was designated as the independent variable, and the trial on which the perceptual shift was reported functioned as the regressor. To control for potential confounds, such as impulsive premature responses and to compensate for diminished statistical power in extreme) response epochs (both early and late) due to reduced sample sizes, our analytical framework was confined to responses within a boundary of one standard deviation from the median (8 ± 2; 87.1% of total trials). As expected, we observed an inverse relationship between the evoked pupil diameter and the trial marking the perceptual switch (β = -0.115, R² = 0.02, p = 1.5 x 10⁻³). Earlier responses showed a positive relationship with higher evoked pupil diameter during the switch epoch (red in Fig. 1D), whereas later responses were associated with a more constricted pupil (blue in Fig. 1D). In summary, these results provide indirect evidence for our hypothesis that ascending neuromodulation – such as LC activity – is associated with the resolution of perceptual uncertainty.

Computational evidence for neuromodulatory-mediated perceptual resolution in a recurrent neural network

Our initial results provided evidence implicating the LC in perceptual switching, however there is also somewhat conflicting evidence showing that simply changing stimulus categories can also induce similar pupillary dilations^16,30,31, suggesting the need for a more mechanistic test of the network reset hypothesis. Along with others^32–35, we have used a combination of computational modelling^27,36, neurobiological theory¹¹ and multi-model neuroimaging^19,27,31,32 to suggest that noradrenaline alters neural gain^28,39, which in turn affects inter-regional communication and flattens the energy landscape traversed by network dynamics allowing the state of the network to jump between perceptual attractors more easily. Whether these signatures of large-scale network reconfiguration are mechanistically related to network reset remains and important open question.

To test whether our hypothesised neuromodulatory mechanism could recapitulate the behaviour we observed in the ambiguous figures task, we trained 50 continuous time recurrent neural networks (RNN) that we constrained to respect Dale’s law (i.e., 80/20 split of purely excitatory/inhibitory units^40,41) to perform a perceptual change detection task analogous to the task performed by our participants (Fig. 2A). We also constrained the input and readout weights to be exclusively positive and only allowed the firing rate of excitatory units to contribute to the readout (see Materials and Methods). Each network was provided with a two-dimensional input , with each row representing the “sensory evidence” for each of the two stimulus categories (Fig. 2A). The task lasted for 50 time-steps: to mimic the linear transition between image categories in our task, each trial began with maximum evidence for one of the two categories and minimum evidence for the other (e.g., u₁ = 1, u₂ = 0), and then changed linearly over the course of each trial such that by the final time-step the evidence for each category had switched (e.g., u₁ = 0, u₂ = 1). At each time point, the network was trained to output a categorical response indicating which input dimension had a higher value.

Figure 2

Following training, all networks achieved near perfect behavioural accuracy (0.99 ± 0.01 Because of the constraint that the input and output weights are strictly positive, we could use their (normalised) value as a measure of stimulus selectivity. Inspection of the firing rates sorted by input weights revealed that the networks had learned to complete the task by segregating the excitatory units into two stimulus-selective clusters, each of which was excited by a single dimension of the input stimulus, which in turn projected to independent dimensions of the network readout (i.e., for all networks input and readout weights were highly correlated; mean = 0.9849, std = .0076; Fig. 2B). The inhibitory units also showed strong stimulus selectivity in their input weights but because of recurrent projections with both excitatory populations their dynamics were only weakly stimulus-selective (Fig. 2C). Because network output was exclusively a function of the excitatory firing rate after training the switch in output from category 1 to category 2 occurred approximately when rate of one of the competing excitatory clusters exceeded the other at time point 25 (i.e., when u₂ > u₁; and v.v.). It is interesting to note that the computational architecture we arrived at by gradient based learning is strikingly similar to classic attractor network models of perceptual decision making^42,43 that have two stimulus selective excitatory populations with recurrent connections to a shared inhibitory pool (see insets of Fig. 2D-F).

With an operational model of the task in hand, we next sought to test our hypothesis about the role of neural gain in perceptual switches. In previous work, we (and others) have argued that the impact of neuromodulators (such as NA) on population-level activity can be effectively approximated by steepening (or flattening) the sigmoid activation function, thus mimicking the effect NA has on neuronal excitability by liberating intracellular calcium stores and/or opening (or closing) voltage-gated ions channels^11,44. To test this mechanism in our model, we manipulated the gain of the sigmoid activation function for both excitatory and inhibitory units by simulating an additional 500 trials in each network between gain values of 0.5 to 1.5 (in steps of 0.1). As predicted, increased gain (red; corresponding to heightened adrenergic tone) lead to earlier ‘perceptual switches’ in the network output whereas low gain caused later switches (Fig. 2E).

Inspection of the firing rates showed that the speeding effect of gain was mediated by a shift in the timepoint at which each cluster of stimulus-selective excitatory units was out competed by the opposing excitatory cluster (Fig. S2). To interrogate whether the earlier switching was caused by increased excitatory or increased inhibitory gain (or both), we re-ran the simulations with the gain manipulation applied exclusively to either the excitatory or inhibitory units in isolation (Fig. 2E-F). Intriguingly, manipulating excitatory gain in isolation slowed switching under reduced gain, but did not speed up switching under heightened gain (Fig. 2E). In contrast, manipulating inhibitory gain in isolation led to both slowed switches under reduced gain and faster switches under heightened gain (Fig. 2F) demonstrating that, in our model, the speeding effect of gain on perceptual switch time was largely mediated by enhanced inhibitory signalling.

Large-scale neural predictions of recurrent neural network model

Having confirmed the behavioural component of our gain modulation hypothesis in our model, we next sought to identify population-level neural signatures of noradrenergic tone. In particular, we sought to test our hypotheses that the speeding effect of gain on perceptual switches is mediated by a flattening of the energy landscape traversed by the network dynamics. Crucially, this shift in perspective from the firing rates of individual units to the dynamics of the RNN as a whole also allowed us to generate a series of predictions that we could later test in functional neuroimaging data.

In recent work^23,45, we have shown that peaks in blood flow within the LC precede large changes in brain state dynamics. Viewed through the lens of dynamical systems theory⁴⁶ in which the brain is treated as a dynamical system whose state space (i.e., an instantaneous snap-shot of the activity of all regions in the system) evolves over time according to the presence (or absence) of attractors, the effect of the LC can be conceptualised as akin to lowering the energy barrier required to escape an attractor defining the current state so as to reach a novel location in state-space²³ (see Fig. S3 for an intuitive example of the energy landscape analysis). In light of the task used in this experiment, stable perception of the first image in a sequence (e.g., the plane in Fig. 2A) is represented as a stable fixed-point attractor, to which all visual inputs somewhat similar to the plane would evolve, as if brain states were a ball drawn by gravity down into a valley. Note that, in this scenario, there also exists an attractor representing a shark to which the brain state could evolve, but the trajectory is blocked by an energy barrier (i.e., it is easier for a plane-like input to be interpreted as a plane). As the incoming image begins to look less like a plane, we hypothesize that this growing perceptual mismatch recruits increased neuromodulatory tone (likely mediated via the LC^47,48). The effect of this increased tone is to ‘flatten’ the energy barrier between the two attractors, and thus to make it more likely that the brain state could fall into the second attractor – i.e., to ‘switch’ the interpretation from a plane to a shark.

To test whether gain was playing an analogous role in the dynamics of the RNN, we leveraged the relationship between probability and energy in statistical mechanics to construct a measure of the energy landscape traversed by the dynamics of each network. Specifically, we calculated the mean-squared displacement of the activity of each RNN at each time point t₀ relative to a reference point t in the trial. We then used Gaussian kernel estimation to construct a distribution over MSD for each timepoint t, and then used this to estimate the energy^23,49 . Following the physical analogy, we conceptualise gain as akin to temperature – systems with a high temperature, have a higher average energy and a more uniform probability distribution. Correspondingly, low temperature systems, are associated with lower average energy and thus a more peaked probability distribution.

We ran this analysis over all trials at each gain value. To obtain a summary measure of the flatness of each landscape we integrated over time to visualise energy exclusively as a function of MSD. In line with our hypothesis the energy required for large movements in state space (i.e., large MSD values) decreased as a function of gain (Fig. 3B) similar to the manner in which increasing temperature increases average kinetic energy (i.e., particle motion). These results reinforce our previous work and clearly demonstrates that the implementation of neuromodulatory-mediated dynamics in the RNN acted in a similar fashion to previously observed patterns in resting-state fMRI²⁶. In addition, our results confirm that the putative impact of the release of noradrenaline from the locus coeruleus can change the manner in which brain states evolve over time by ‘turning up the temperature’ which ‘flattens’ the brains energy landscape, which in turn helps to facilitate the navigation of otherwise difficult state-to-state transitions²⁶.

Figure 3

Having confirmed our hypothesis that the energy landscape flattens under increased gain, we turned from the energy-landscape traversed by the population dynamics to the populations dynamics themselves. Specifically, we performed a Principal Component Analysis (PCA) on the concatenated activity of the network at gain = 1. The set of PCs was low-dimensional, with 75.98 ± 5.79% of the variance explained by the first principal component (PC₁). Based on this information, we projected the trial averaged activity at each gain value and timepoint onto the first PC. The resultant low dimensional trajectories all showed a change in direction around the timepoint of the switch in network output from category 1 to category 2 (and v.v.; Fig. 3B inset), reminiscent of a system jumping between attractors which, occurred earlier (and later) as a function of heightened (and lowered) gain (Fig. 3C). This switch not only occurred sooner as a function of heightened gain (Fig. 3B inset), but it also occurred at a higher neural “speed” with the velocity of the trajectory peaking sharply at the point of the switch under high gain, whereas the transition between states was relatively gradual under low gain (Fig. 3B).

The low-dimensional signature of ambiguity resolution and perceptual change

Having confirmed our hypothesis about the speeding effect of gain in our RNN model we next sought to test the predictions in the human brain – i.e., whether the increase in neural speed and the flattening of the energy landscape observed in the RNN were also present in functional neuroimaging data. To this end, we re-analysed an existing BOLD dataset collected while participants performed a similar version of the ambiguous figures task to identify the low-dimensional patterns that occurs during the perceptual change.

We were, however, left with a dilemma: RNNs provide a proof-in-principle of how computations can be instantiated in networks, however there are key differences between artificial neural networks and the human brain that require careful consideration^50,51. While both RNNs and the brain are thought to compute through dynamics⁵², the human brain is comprised of highly specialised neural circuits that have been shaped over evolutionary time to perform a range of highly idiosyncratic functions that matter for adaptive behaviour⁵², but aren’t necessarily related to task-switching. So where in the brain should we look for the same low-dimensional signatures we observed in the RNN as a function of gain? Rather than select a particular region a priori, we instead opted for a data-driven approach – principal components analysis (PCA) – which summarizes regional timeseries concatenated across all subjects and trials into a set of low-dimensional patterns that can then be interrogated in a similar fashion to the activity of the RNN (see Methods for details). Consistent with previous work⁵³, a small number of principal components (PCs) mapped onto distributed regions across the brain (Fig. 4A) and explained a substantial proportion of the variance observed in the task (PC_1-3 explained 32% of the total variance).

Figure 4

To isolate the low-dimensional component that best reflected the task (Fig. 1A), we performed a principal component regression⁵⁴ that modelled the switch point of each trial using the loadings of the top 3 PCs calculated from fMRI data. PC₁ was not selectively aligned with switches, both PC₂ and PC₃ showed a pronounced, isolated peak around the switch point across trials (Fig 4B), with PC₂ showing the most robust task-related engagement. To ensure that these results could not be explained by the spatial autocorrelation inherent within the PC maps, we created a null distribution of regression coefficients calculated using the same statistical model but with block-resampling applied to the switch times in the design matrix. The dotted grey line in Fig. 3C denotes the 95^th percentile of the null distribution, and clearly shows that the engagement of both PC₂ and PC₃ during the switch point was greater than to be expected by chance. Furthermore, to validate that the perceptual switch was predominantly represented by PC₂ and PC₃ (Fig. 4D), we conducted a regression with these two PCs as predictors and the evoked activity derived from the original BOLD time-series as the dependent variable. The resulting variance accounted for was 88% (R² = 0.88, β = 0.99, p = 9.2x10^-178).

To determine whether PC₂ or PC₃ was a better index of perceptual switching, we then correlated the spatial loadings of PC₂ and PC₃ with the spatial map associated with the term “switching” from a meta-analysis performed on the neurosynth database⁵⁵. We observed a significant positive correlation between the map for “switching” and both PC₂ (r = 0.453, p = 3.041x10^-18), and PC₃ (r = 0.115, p = 0.037), however the correlation for PC₃ was much lower than PC₃, suggesting that PC₂ was a better match for “switching”. The spatial map of PC₂ was also positively correlated with other terms putatively associated with the ambiguous figures task (notably, “effort”, “load” and “attention”; all r > 0.2; and not with “episodic”, which was included as a negative control), a partial correlation analysis revealed that PC₂ was selectively associated with “switching” and “attention” (Fig. 4C). Given the multifaceted nature of the ambiguous figures task, the convergence between brain maps for “switching”, “attention”, and “effort” was to be expected and we therefore did not try to dissociate them in further analysis.

Before turning to the predictions of the RNN, we first sought to further validate the face validity of focusing on a limited number of principle components. In previous work, we have linked the impacts of NA on systems-level neural dynamics to alterations in network topology^23,56,57, with NA increasing large-scale network integration. Given that PCA naturally captures patterns of covariance between regions, we expected to see that the observed time signatures of PC engagement at the switch-point should coincide with similar measures of network integration. To test this hypothesis, we clustered the time-averaged functional connectivity matrix using a hierarchical modular decomposition approach (see Methods) – doing so revealed three main clusters (Fig. 4E). For each participant, we used this matrix and the three clusters to estimate the amount of integration (using the participation coefficient) and segregation (using the module-degree Z-score, see Methods) of each region. We then correlated a joint histogram of these measures with the sum of subject-specific regression coefficients for PC₂ and PC₃ and observed a robust correlation with integration (Fig. 4F; p < 0.05 following permutation testing). These results clearly demonstrate the highly convergent nature of PCA and our previous network-based approaches.

Confirmation of model predictions in whole-brain BOLD data

Based on the patterns observed in the RNN (i.e., those in Fig. 2-3), we hypothesized that the energy landscape topography would decrease, and the velocity of the low dimensional brain patterns would peak at the switch point. Given the prominent role in switching we focused our analysis on the PC₂ time series. To estimate the energy landscape, we first estimated the mean displacement of PC₂ by averaging the β value around the switch-point and then divided this term by the logarithm of the inverse probability of the loading of PC₂, which was also inferred from the GLM. Using this approach, we observed that PC₂ was maximally displaced at the perceptual change (Fig. S4C), suggesting that the brain state showed a substantial shift from baseline during the perceptual change. Energy (log[1/p_switch]; see Methods) showed a U-shaped pattern around the perceptual change point – i.e., with a minimum value in the perceptual change along with the first and last images (Fig. 5B-D). To relate this measure to the energy landscape framework, and to control by the specific displacement occurring at each image, we then calculated the ratio between energy and the mean displacement (i.e., energy landscape ’depth’; Fig. 5A). As predicted, the brain-state reduced the amount of energy per displacement towards its minimum around the perceptual change (Fig. 5B). We interpret this set of results as the system flattening the energy landscape, reducing the energy (i.e., higher system changes become more common) required for large displacement values, effectively generating a ‘network reset’^12,14,58 of the brain state, which ultimately facilitated an updating of the content of perception.

Figure 5

To analyse speed-evoked changes in brain trajectories, we used a GLM to analyse each PC time series as a function of each perceptual switch. Our design matrix included the first and last images seen in each set, along with the three images leading up to the switch, the switch trial itself and the three images following the switch (see Methods for details). This approach thus allowed us to track the low-dimensional signature of the brain through the processing and resolution of perceptual ambiguity. As predicted (Fig. 5C), we found evidence that PC₂ showed a peak in velocity at the change point (Fig. 5D) providing confirmatory evidence that the low-dimensional brain state dynamics observed in whole-brain fMRI were highly similar to those observed in the trained RNN.

Discussion

Here, we studied the relationships between the ascending arousal system, low-dimensional trajectories and energy landscape dynamics during a perceptual switch task. Our results provide evidence that the ascending arousal system is involved in the modulation of dynamic brain state topography during switches in perceptual interpretation. We found that pupil diameter tracked with the task stimuli and was directly related to the resolution of the ambiguity of the image being presented (Fig. 1). Next, we confirmed that this process could be replicated using a simple RNN Model (Fig. 2), in which the functional behaviour of the network could be directly modulated by neuromodulatory-mediated alteration of activation function (Fig. 3). Finally, we confirmed that low-dimensional dynamics around the perceptual change could be inferred from whole brain neuroimaging data (Fig. 4) and that the low-dimensional signatures of the impacts of NA-mediated gain changes on RNN behaviour could also be inferred from whole-brain functional neuroimaging data (Fig. 5). Together, these results provide evidence of how the ascending arousal system facilitates perceptual dynamics during ambiguity by flattening the energy landscape thereby reducing the energy needed to reset the system topology in an adaptive and task dependent manner.

The relationship between perception and pupil diameter is consistent with the role of the noradrenergic system in cognition and attention^39,44. For instance, the LC dynamically changes its activity according to external and cognitive demands imposed on the system^{29,39,44,59,60}. Our results demonstrate a more precise role for LC-mediated alterations in neural gain in the transition between the two images. Even though the image changes occurred linearly across each trial, such that the middle trial contained a perfect blend of both images, the perceptual ambiguity is hypothesized to place demands on the system that we propose are resolved through the recruitment of the LC and noradrenaline, which then subsequently increase brain-wide communication by increasing the gain and signal transmission in targeted brain regions^11,27,60,61. As the pupil diameter increased as a function of perceptual ambiguity, which rose sharply in the few images prior to a reported perceptual change (Figure 1D), and pupil diameter is an indirect marker of the noradrenergic system^9,15, we thus provide early evidence for the role of the ascending neuromodulatory activity to facilitate perceptual changes and decision making under perceptual ambiguity^17,21,22. Taken together, we have provided evidence that the ascending arousal system, measured indirectly through pupillometry, can modulate the conscious perception of ambiguous images by facilitating the arbitration of uncertain information decision-making processes.

A core neuroanatomical property of the LC noradrenergic system is that a relatively small number of neurons (~fifty thousand in an adult human) projects to almost all of the regions in the brain^20,62. This organisation implies that the LC acts as a low dimensional modulator of the much more high-dimensional cerebral cortex. Subtle changes in the activity of LC can have significant effects on how different brain regions communicate^{27,44,59,63–65}. Based on this, we hypothesized that the ascending neuromodulatory system needs to increase its activity in order to facilitate network communication and brain state transitions when the sensory information is limited or presents ambiguity. In line with this hypothesis, we showed that increasing the gain of the sigmoid transfer function – which at a neuronal level is mediated by NA increasing intracellular calcium and opening (or closing) voltage-gated ions channels^11,44 – in our RNN model of perceptual switching had a speeding effect on the timing of perceptual switches, an effect that we found was largely mediated by enhanced inhibitory signalling. At the large-scale level, increasing gain led to a flattening of the energy landscape traversed by the dynamics of the RNN increasing the likelihood of large shifts in state space analogous to increasing the ‘temperature’ of a physical system^36,66. Indeed, analogous to an increase in kinetic energy we observed an increase in the velocity around the switch point as a function of gain.

In line with the predictions of the RNN in our analysis of the BOLD data, we showed that the velocity of the low dimensional brain state trajectory most associated with perceptual switching increased significantly during the point of reported perceptual change in comparison (Fig. 4B), which we interpret as the brain moving from one attractor to another (Fig. 5A). Importantly, we showed that around the perceptual switch, the energy needed for each unit of change in brain state (i.e., displacement) is smaller than at other points in the task (Fig. 5A-B). Under the energy landscape framework^23,38, this shows that the landscape is flattened, and the energy required to transition between states is reduced. Together with the pupillary findings (Fig. 1), the computation model (Fig. 2-3), and replication of former results^23,38 (Fig. 4E-F), we propose that the ascending neuromodulatory system is responsible for the large-scale flattening of the energy landscape which facilitates the change in reported perceptual changes.

This work is not without limitations. First, the pupil diameter dataset and the fMRI analysis came from different participants, such that the link between the pupil diameter and the fMRI results is inherently indirect. Moreover, the specificity of the pupil diameter as a marker of the LC activity is under active debate¹⁵. For instance, there is evidence supporting a role of the superior colliculus, the dorsal raphe nucleus and central cholinergic system in driving pupil dilations^29,67,68. There is uncertainty regarding whether these other nuclei are directly related to pupil dilation, or only indirectly via their connections with other neural regions and nuclei. Despite this, we believe that our pupillometry dataset captures the ascending neuromodulatory role on perceptual ambiguity as there is strong evidence showing that pupil diameter is a reliable marker of noradrenergic activity during evoked cognitive tasks^{19,63,69–71}. Additionally, the sample size of our fMRI study makes it difficult to generalize our results. In spite of this, we believe that the converging evidence from the pupillometry dataset, the fMRI dataset and the computational model, supports the role of the ascending neuromodulation on the resolution of perceptual ambiguity hypothesized here. Future work is needed both in humans, with higher sample sizes utilizing fMRI and eye-tracking recordings, as well as animal studies, to directly modulate and record the LC activity in a perceptual ambiguity task of the kind used here.

Conclusion

In summary, we provide computational and empirical evidence for the association between neuromodulation, pupil dilation, energy landscape flattening and the resolution of perceptual ambiguity. Our results strengthen our understanding of the complex, distributed neurobiological processes that in combination facilitate our ability to resolve perceptual ambiguity. Specifically, we suggest that the gain-modulated dynamics of the brain system use the widespread projections of the noradrenergic arousal system to mediate the reconfiguration of the systems-level organisation of cortical network architecture^62,63 to complete challenging perceptual tasks. Rather than working alone, we further hypothesize that the noradrenergic system cooperates and competes with other distinct arms of the ascending arousal system to together facilitate complex, adaptive network dynamics in the brain.

Methods

Overview

There were two independent groups analysed in this study: 35 subjects performed a perceptual decision-making perceptual task while pupil diameter was recorded; and a separate group of 17 subjects performed a version of the task adapted for the MRI scanner.

Perceptual Task

Twenty picture sets were used in which line drawings of common objects morphed over 15 iterations into a different object (Fig. 1A). Picture sets were selected from a larger set validated in an earlier study²⁴. Picture sets for the current study were selected with the criterion that all sets were perceived categorically in the normative study (i.e., that the majority of participants in the normative study categorized each picture they saw as either the first object or second object in the set¹). Selecting only the categorically perceived image sets guaranteed that pictures in the middle of the morphing sequence were not simply ‘noisier’ than pictures at the beginning or end. In other words, the ambiguous images were still easily categorised by participants as either object 1 or object 2. All images were a standard size (316x316 pixels) and were displayed on a white background. In addition, participants were presented with two kinds of control picture sets to ensure that they were responding to changes in the pictures in the set rather than simply to the position in the set (e.g., always switching after the 8th picture). In these control picture sets, a salient deviating picture was presented either after three pictures or after thirteen pictures resulting in an early or late abrupt shift. Those sets served as controls and were not analysed further.

The picture morphing task consisted of five experimental runs. We randomized the order in which the picture sets were presented in each run and kept this randomized order consistent across participants. Picture morphing in each picture set occurred over fifteen discrete steps, each corresponding with the acquisition of a whole-brain image. Each picture within a set was presented for two seconds. Pictures were randomly intermixed with eight inter-stimulus-intervals (2, 4, 6 or 8 seconds) during which participants saw a fixation cross. Participants provided their responses in the scanner using two buttons on a four button Cedrus fibre optic system. In a two-alternative forced-choice task, participants were asked to press the first button when they ‘saw the first object' and the second button when they ‘saw the second object' – this ensured that there was not a motor confound present on only the switch trials. All participants were ignorant as to the identity of the second object in each picture set. At the end of each set of 15 images the word END was presented for 2 s to indicate that the next picture set would begin shortly.

Participants

A total of seventeen (6 male) neurologically healthy participants with normal or corrected to normal vision took part in this study (mean age 27.65 ± 8.01). Fifteen were right-hand dominant. A separate cohort of 35 participants performed the task while simultaneous pupil diameter was recorded using an eye tracker device (SR Research, 1000 Hz). None of the participants had a history of brain injury. Participants received $30 for their participation. All participants provided informed consent prior to participation. The research protocol was approved by the Office of Research Ethics at the University of Waterloo and the Tri-Hospital Research Ethics Board of the Region of Waterloo in Ontario, Canada.

Pupillometry

Fluctuations in pupil diameter of the left eye were collected using an Eyelink 1000 (SR Research Ltd., Mississauga, Ontario, Canada), with a 1 kHz sampling frequency. Blinks, artifacts, and outliers were removed and linearly interpolated⁷². High-frequency noise was smoothed using a second-order 2.5-Hz low-pass Butterworth filter. To obtain the pupil diameter average profile, data from each participant were normalized across each trial (corresponding to the 15 consecutive image set). This allowed us to correct for low-frequency baseline changes without eliminating the load effect and baseline differences due to load manipulations^73,74.

Recurrent Neural Network Modelling

We used PyTorch⁷⁵ to implement and train 50 continuous-time recurrent neural networks that we constrained to respect Dale’s law (N_E+I = 40, 80% excitatory N_E = 32, and 20% inhibitory N_I = 8) using the procedure set out in Ref #⁴¹. The dynamics of each network evolved according to the following system of ordinary differential equations:

Where x ∈ ℝ^N×1 represents the sub-threshold activation of each unit, u ∈ ℝ^2×1 the external input into the network, W^rec ∈ ℝ^40×40 the recurrent weights, Wⁱⁿ ∈ ℝ^40×2 the input weights, and τ the time constant which we set to 100ms. In addition to task input each unit in the network received independent Gaussian noise Ν with 0 mean and standard deviation (with σ_rec set to 0.1). The subthreshold activation x was converted into a vector of firing rates by applying a sigmoid function where g ∈ ℝ^N×1 is a vector containing the gain control parameter of each unit’s activation function that was multiplied element wise with x. Network outputs z ∈ ℝ^2×1 were given by a linear readout of the excitatory population’s firing rate z = W^outr_E. Where . The network’s choice at each time point was the maximum of the two-dimensional output z.

We imposed Dale’s law on the recurrent weights of the network by parametrising the weight matrix with a mask W^mas ∈ ℝ^40×40 which contained zeros in the leading diagonal (removing self-connections), +1 in all non-diagonal entries of the first 32 rows/columns and −1 in the remaining 8 rows/columns. We obtained the constrained recurrent weight matrix by multiplying the absolute value of the trained weights element wise with the mask thereby imposing an 80/20 E/I ratio. Similarly, in line with the finding that long range connection in mammalian cortex are purely excitatory, we constrained the projection of the input to the network and the readout projection to be strictly positive by taking the absolute value of the trained input and output weights .

Following standard practice⁴¹, we simulated the network by discretising the system of ODEs using a forward Euler integration scheme with Δt set to 20 (denoting ).

Each network was trained by optimising , to minimise a cross entropy loss function through 2000 iterations of back propagation through time⁷⁶ (BPTT) with ADAM⁷⁷. Batches consisted of single trials which for our simple task (described below) was sufficient for each network to converge on near perfect behavioural accuracy. All training was performed with the gain control parameter set to 1.

The task consisted of a simple change detection paradigm analogous to the task performed by our human participants. Specifically, at each time point the network was fed a two-dimensional input with each row representing the “sensory evidence” for each of the two stimulus categories. The task lasted for 50 time-steps beginning with maximum evidence for one of the two categories and over the course of each trial changed linearly such that at time-step 26 the sensory evidence for each stimulus category changed to slightly favour the second stimulus category and by the final time-step consisted of maximum evidence for the second stimulus category . We trained the network to output a response for stimulus category 1 whenever u₁ > 0.5, and u₂ < 0.5, and category 2 whenever u₁ < 0.5, and u₂ > 0.5.

Following training we simulated a further 500 trials over the gain parameter range [0.5, 1.5] in steps of 0.1. To study how the population dynamics of the trained networks changed as a function of gain we performed a Principal Component Analysis (PCA) on the concatenated activity of the network at gain = 1. The set of principal components was highly low-dimensional, with 75.98 ± 5.79% of the variance explained by the first principal component (PC₁) and 13.75 ± 4.13 % explained by the second. We then projected the trial averaged activity at each gain value at each timepoint onto the top PC. The resultant low dimensional trajectories all showed a change in direction around the timepoint of the switch in network input from category 1 to category 2. We then calculated the gradient of the trajectory around the switch point observing a consistent increase in velocity as a function of gain (Fig. 2D).

Leveraging previous work from our group²³ we constructed a measure of the energy landscape traversed by each network through an analogy to the relationship between probability and energy in statistical mechanics⁴⁹ which is given by:

Where p_i denotes the probability of each state, E_i the energy of each state, β thermodynamic beta, and z the canonical partition function. Solving for E_i we obtain:

We then calculated the mean-squared displacement (MSD) of the activity of the RNN at each time point t₀ relative to a reference point t in the trial.

The MSD is a measure of the average deviation in the activity, x_t = [x_1,t x_2,t, …, x_r,t] of each unit n, with respect to the activity at reference point t0. We then constructed a distribution over MSD values at each timepoint using Gaussian kernel estimation , where . As Σ_i p_i = 1 by construction z = 1, which after setting β = 1 gives:

Following the physical analogy, we think of the MSD of the RNN activity as akin to temperature with large displacements, like systems with a high temperature, having a higher average energy and thus a more uniform probability distribution. Correspondingly, low temperature systems, are associated with lower average energy and thus a more peaked probability distribution. We ran this analysis over all trials at each gain value and then correlated the resulting energy, E, with gain (Fig. 2D).

MRI Data

Functional data were acquired using gradient echo-planar T2*-weighted images collected on a 1.5 T Phillips scanner located at Grand River Hospital in Waterloo, Ontario (TR = 2000 ms; TE = 40 ms; slice thickness = 5 mm with no gap; 26 slices/volume; FOV = 220x220 mm²; voxel size = 2.75x2.75x5mm³; flip angle = 90°). Each experimental run consisted of 26 slices per volume and 285 volumes. At the beginning of each run, a whole brain T1-weighted anatomical image was collected for each participant (TR = 7.4 ms; TE = 3.4 ms; voxel size = 1x1x1mm³; FOV = 240x240 mm²; 150 slices with no gap; flip angle = 8°). The experimental protocol was programmed using E-Prime experimental presentation software (v1.1 SP3; Psychology Software Tools, Pittsburgh, PA). Stimuli were presented on an Avotec Silent Vision™ fibre-optic presentation system using binocular projection glasses (Model SV-7021). The onset of each trial was synchronized with the onset of data collection for the appropriate functional volume using trigger pulses from the scanner.

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, including regression of estimated motion parameters), the unwrapped EPI images were then smoothed at 6-mm FWHM, and nonlinearly co-registered 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⁷⁸, 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 five runs (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⁷⁹. Finally, in keeping with previous time-resolved connectivity experiments⁸⁰, 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⁸¹; (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/); (c) 28 cerebellar regions from the SUIT atlas⁸² for each participant in the study.

Neuroimaging Analysis

In order to analyse task evoked activity related to stimulus presentations, we first performed a principal component analysis (PCA)⁵³ on the pre-processed BOLD time-series (per subject/session), to extract orthogonal low-dimensional time-series. The top 3 PCs explained ~30.6% of the variance. The time-series of these PCs was entered into a general linear model, in which we modelled the following 9 event-types across an entire session, centred around the perceptual switch point, which changed on a trial-by-trial basis: the first two images (modelled as a single regressor), the seven images surrounding each perceptual change (i.e., the switch trial and the three images surrounding the change point, modelled as seven separate regressors) and the last two images (modelled as a single regressor). Each of the event onset times was also convolved with a canonical hemodynamic response function. This left us with nine unique β values per principal component, which we could use to determine how each PC differentially engaged as a function of the task. To test the hypothesis that the rate of change of PC engagement peaked at the perceptual change point, we calculated the difference between the β value for each of the top 3 PCs for each of the 9 event-types, and then plotted the resultant series in order to identify whether a peak occurred at the perceptual switch point (i.e., the middle β value in the series). A block-resampling null (n = 5,000 permutations) was used as a permutation test (p < 0.05).

Spatial maps associated with the terms: “switching”, “effort”, “attention”, “perception” and “load” were downloaded from the neurosynth repository⁵⁵ and mapped into our parcellation space by calculating the mean value within each independent parcel. These values were then correlated with the spatial loading of each of the top 3 PCs. A separate partial correlation analysis was conducted in which the same correlation was estimated after controlling for each of the other spatial maps.

Topological analyses

A hierarchical modularity approach was used to collapse the mean time-averaged correlation matrix across participants into a set of four spatially non-overlapping modules. Briefly, this involved running the Louvain modularity algorithm, which iteratively maximizes the modularity statistic, Q, for different community assignments until the maximum possible score of Q has been obtained.

The community assignment for each region was then estimated 500 times across a range of γ values (0.5-2.0, in steps of 0.1). In order to identify multi-level structure in our data, we repeated the modularity analysis for each of the modules identified in the first step⁸³. Finally, a consensus partition was identified using a fine-tuning algorithm from the Brain Connectivity Toolbox (http://www.brain-connectivity-toolbox.net/). We subsequently used this final module assignment to estimate the cartographic profile of the each participant’s time-averaged adjacency matrix⁵⁷. Specifically, we estimated Integration using the participation coefficient, which quantifies the extent to which a region connects across all modules (i.e., between-module strength⁸⁴), and Segregation using the module-degree Z-score. These measures were entered into a joint-histogram (101x101 unique bins, equally-spaced between 0-1 [for Integration] and -1 and 1 [for Segregation]). The value within each bin of this joint histogram was then correlated with the combined regression weights of PC₂ and PC₃ for each subject. A permutation test that scrambled the order of participants was used to assess statistical significance (p < 0.05).

Brain-state displacement and the energy landscape

To quantify the change in the evoked BOLD activity following each stimulus we calculated the main BOLD displacement (MBD). The MBD is a measure of the absolute evoked deviation in BOLD activity. The evoked activity is measured through a general linear model using a canonical hemodynamic response function convolved on a design matrix. We are interested in the probability, p(MBD, re), that we will observe a given displacement in BOLD at a given regressor re. The probability is calculated through the null model of the general lineal model (the probability that the observed evoked value of the corresponding region is different from 0). As described above we then calculated the energy for each displacement value as . Finally, to measure the surprise per displacement, we divided the absolute β for PC2 from E_MBD,re for each regressor re (Figure 3).

Supplemental Materials

Figure S1

Figure S2

Figure S3

Figure S4