A temporal sequence of thalamic activity unfolds at transitions in behavioral arousal state

Awakening from sleep reflects a profound transformation in neural activity and behavior. The thalamus is a key controller of arousal state, but whether its diverse nuclei exhibit coordinated or distinct activity at transitions in behavioral arousal state is unknown. Using fast fMRI at ultra-high field (7 Tesla), we measured sub-second activity across thalamocortical networks and within nine thalamic nuclei to delineate these dynamics during spontaneous transitions in behavioral arousal state. We discovered a stereotyped sequence of activity across thalamic nuclei and cingulate cortex that preceded behavioral arousal after a period of inactivity, followed by widespread deactivation. These thalamic dynamics were linked to whether participants subsequently fell back into unresponsiveness, with unified thalamic activation reflecting maintenance of behavior. These results provide an outline of the complex interactions across thalamocortical circuits that orchestrate behavioral arousal state transitions, and additionally, demonstrate that fast fMRI can resolve sub-second subcortical dynamics in the human brain.


INTRODUCTION
Striking behavioral changes accompany arousal from sleep, as we transition from unresponsiveness to engaging with our sensory environment. These behavioral changes reflect a fundamental shift in neural dynamics throughout the cortex and subcortex. Extensive studies of arousal state transitions have identified characteristic electroencephalography (EEG) signatures of arousal [1][2][3][4] . Arousal is marked by changes in brain rhythms 2,5-10 and functional connectivity [11][12][13][14] throughout thalamocortical networks. Although the differences in neural dynamics between stable states of sleep and wakefulness have been well characterized, how thalamocortical networks implement transitions between behavioral states is not well understood. Furthermore, while arousals are often characterized by their EEG signatures, behavior is not always linked to a specific EEG signature [15][16][17] , which leaves open the question of which neural dynamics mediate behavioral state transitions. In particular, the specific sequence of neural activity that unfolds at the moment of behavioral state transitions has not yet been mapped throughout large-scale thalamocortical networks.
The thalamus is made up of many diverse nuclei with unique structural and functional properties [51][52][53] , and whether their individual activity is coordinated or distinct at state transitions has not yet been established. Invasive studies have provided important insights into the causal role of several nuclei of the thalamus, showing that individual nuclei can modulate arousal. However, the joint dynamics across multiple nuclei are not well understood, and addressing this question requires recording from many nuclei simultaneously. Traditional noninvasive methods such as conventional fMRI can provide whole-brain simultaneous imaging, but previously lacked the temporal resolution necessary to capture fast dynamics. A key open question is whether behavioral state transitions are due to synchronous changes throughout thalamocortical networks, or a specific cascade of neural events.
To address this question, we used high temporal resolution (fast) fMRI to delineate the thalamocortical dynamics underlying behavioral arousals. Advances in pulse sequence technologies now enable acquiring whole-brain fMRI data at fast rates [54][55][56] (e.g., repetition time (TR) under 400 ms), and investigations of the hemodynamic response have demonstrated that fast fMRI data can resolve neural dynamics on timescales of hundreds of milliseconds, in cortex and thalamus [57][58][59][60][61][62][63] . We exploited these fast fMRI techniques to track rapid thalamic dynamics at the moment of changes in behavioral arousal state. First, we used fast fMRI at 3T to investigate the moment of behavioral arousal, defined as the time when behavior returns after a period of inactivity, and found thalamic activation prior to arousals. Then, to delineate activity within distinct thalamic nuclei, we collected a second dataset using fast fMRI at 7T to image nine individual thalamic nuclei with sub-second resolution. We identified a temporal sequence of activity that unfolds across thalamic nuclei just prior to the moment of behavioral arousal, and further determined that this sequence precedes a brain-wide reorganization of thalamocortical dynamics that reflects subsequent behavioral state. These results precisely identify a cascade of neural events that unfold as the brain switches between behavioral arousal states.

Identifying moments of behavioral arousal state transition across sleep and wakefulness
We aimed to identify the thalamocortical dynamics underlying behavioral arousal state transitions by imaging participants spontaneously transitioning between sleep and wakefulness during nighttime fMRI scans. Arousal from sleep occurs spontaneously throughout the night 64 , and we observed occasional spontaneous awakenings in participants sleeping in the scanner. To track behavioral arousal state without a stimulus that would influence arousal directly, all subjects performed a self-paced behavioral task in which they pressed a button with every breath 17,65 . Behavioral arousal was defined as the first button press after at least 20 s (seconds) of unresponsiveness. Because movement sometimes accompanies arousal from sleep, arousals with high motion (>0.3 mm) were excluded to minimize fMRI artifacts.  (8)(9)(10)(11)(12)(13). Top: the alpha power increase coincides with the return of behavior. Power is temporally smoothed in 2-s sliding windows. Bottom: the spectrogram demonstrates a return of occipital alpha directly following behavioral arousal. b) Top: the mean alpha power significantly increases at behavioral arousal, consistent with a transition to the awake state (6 subjects, 66 arousals). Power is temporally smoothed in 2-s sliding windows. Shading is standard error. Bottom: spectrogram shows specific increase in alpha power during behavioral arousal.
z With this approach, we identified 66 behavioral arousals from 6 subjects (mean=11, std=11.78, min=1, max=29) with simultaneous EEG and fast fMRI at 3T. Traditional sleep scoring from EEG recordings revealed that most of these arousals occurred from NREM sleep: 15 s before the arousal event, 20% were awake, 3% were in N1 sleep, 51% were in N2 sleep, and 26% were in N3 sleep (Fig. S1a). The rigid 30-s intervals used in sleep-scoring do not always reflect the transient dynamics of arousal state transitions, and since some of the segments classified as awake likely corresponded to N1 sleep, the proportion of behavioral arousals representing transitions out of sleep may have been higher than the observed 80%. However, this result also demonstrated that behavioral arousal state transitions could occur at low rates during EEG-defined wakefulness.
We calculated the alpha (8-13 Hz) power in the occipital electrode to evaluate electrophysiological dynamics coupled to behavioral arousal, because occipital alpha defines wakefulness in polysomnography 2 . We found that behavioral arousal was accompanied by a rise in occipital alpha power ( Fig. 1; p<0.05; paired t-test), consistent with its link to awakening from sleep. In keeping with prior observations that behavioral and electrophysiological measures provide different definitions of arousal 16,17 , these results confirmed that our definition of behavioral arousal captured a state transition linked tobut not identical toawakening from EEGdefined sleep.

Thalamus activates before behavioral arousal
Our fast fMRI acquisition (TR=367 ms at 3T) enabled us to delineate temporally precise signals within the thalamus and cortex during behavioral arousal. To that end, we extracted the mean arousal-locked fMRI time series in each region. Activity in both the thalamus and cortex increased slightly in the ten seconds before arousal, but then their activity diverged. We observed a striking increase in thalamic activity at behavioral arousal, whereas the cortical signal decreased (Fig. 2a), consistent with prior studies 49,50,66 . Notably, the latency and waveform of the arousal-locked fMRI signatures differed across the thalamus and the cortex. The thalamus activated before behavioral arousal, with a latency (defined as 10% of maximum signal) of −1.65 s before arousal (95% confidence interval, CI: (−3.57 s, 0 s)), while the cortex deactivated afterward with a latency of 4.77 s after arousal (95% CI: (3.76 s, 5.60 s)). These results demonstrated that the increase in thalamic activity began seconds earlier than the return of behavior and diverged from the pattern seen in the cortex.
While most behavioral arousals were clearly from deeper stages of sleep (N2 and N3), our data did include some behavioral arousals from light N1 sleep and wakefulness (Fig. S1a). Therefore, we investigated whether this same thalamocortical pattern was present when analyzing only arousals from the deeper stages of sleep, yielding 50 arousals from 5 subjects. We found that the arousal-locked signatures in the thalamus and cortex were broadly conserved (Fig. S1b). The thalamus activated −2.8 s before arousal (10% latency; 95% CI: (−3.50 s, 0 s)), and the cortex deactivated 4.92 s after arousal (10% latency=95% CI: (3.67 s, 5.32 s)). Thus, we concluded that the moment of behavioral arousal was preceded by a large increase in thalamic activity and a subsequent decrease in cortical activity.
To investigate whether these arousal-locked patterns were present throughout the cortex or within specific areas, we used the Desikan-Killiany atlas 67 to extract 30 cortical regions of interest (ROIs). The vast majority of the cortex did not show a significant increase at arousal (26/30, two-sided t-test Bonferroni corrected); however, there were a few exceptions (Fig. 2b, S2), including the caudal anterior cingulate (cACC) and posterior cingulate cortex (PCC). The cACC and PCC activity increased before arousal (latency=−4.53 s, 95% CI: (−6.14 s, −0.83 s), and latency= −4.91 s, CI: (−6.06 s, −2.80 s), respectively) and showed higher correlations with the thalamic arousal-locked fMRI signal than other cortical regions (Fig. 2c). The cingulate cortex thus had unique arousal-locked activity within the cortex, resembling the thalamus in its profile of an early, large activation that preceded the moment of behavioral arousal.

Arousal-locked dynamics across individual thalamic nuclei and cortical regions detected at 7T
Our results confirmed a striking increase in thalamic activity at the moment of behavioral arousal and further identified the temporal dynamics of this activity, showing that thalamic signals precede behavioral state changes and cortical suppression by several seconds. However, the thalamus is a heterogeneous structure made up of many functionally distinct nuclei, and their respective roles in arousal dynamics is not well understood. The temporal signal-to-noise ratio of individual thalamic nuclei is low at 3T. Therefore, we repeated our study using 7T fMRI to enhance subcortical signals and enable imaging of individual thalamic nuclei. Using the same nighttime scanning procedure and behavioral task, we recorded 97 behavioral arousals from 13 subjects (mean=7.61 per subject, std=5.78, min=1, max=22). To minimize signal blurring across nuclei, we analyzed only the thalamic nuclei which had a 90% probability or higher of filling at least one functional voxel in all of our subjects 68 , yielding 9 nuclei of interest (Fig. 3a). In one subject, we again confirmed an electrophysiological arousal state transition during behavioral arousal by simultaneously recording EEG signals 69 to measure occipital alpha. As at 3T, we found that occipital alpha power rose during behavioral arousal (Fig. 3b (p=0.01 paired t-test), suggesting our behavioral definition captured similar arousal changes in the 7T setting. We then repeated our analysis of the thalamocortical activity locked to arousal, and replicated our observation . Thalamus activates seconds before behavioral arousal, while most of cortex deactivates after arousal. a) Arousal-locked blood-oxygenation-level-dependent (BOLD) signals in the whole cortex and thalamus. Thalamus (red) begins to activate (latency denoted by red arrow) seconds before behavioral arousal (vertical dashed line), while cortex (blue) deactivates afterwards (blue arrow). Shading represents standard error. b) Most of the cortex deactivates during arousal, but a subset of regions, including the caudal anterior cingulate cortex (cACC, purple) activate before arousal similarly to the thalamus. The rostral middle frontal (rmF, pink) and superior temporal (sT, dark blue) are shown as representative regions from frontal and temporal lobes. Results from all cortical ROIs in can be seen in Fig. S2. c) Correlation between each individual cortical region and the whole cortex or thalamus. While cortical activity is largely distinct from the thalamic rise, a subset of cortical regions are correlated with thalamus, notably the caudal anterior cingulate and the posterior cingulate. that the thalamus activated before arousal (10% latency=−2.60 s, 95% CI (−6.02 s, −1.08 s)) and the cortex deactivated afterward (10% latency=4.38 s, 95% CI:(3.67 s, 5.25 s); Fig.3c, S3). Thus, these thalamic and cortical activity patterns were remarkably preserved across the two experiments, with a consistent rise in thalamic activity beginning before behavioral arousal transitions.
While our initial results identified a major divergence in overall cortical and thalamic activity at behavioral arousal, we next aimed to determine whether these dynamics were consistent within individual thalamic nuclei and cortical regions or unfolded heterogeneously across subregions. We found that all thalamic nuclei activated during behavioral arousal, and inter-nucleus correlation was higher during behavioral arousal than during wakefulness (p<0.001; r=0.31 in wakefulness and r=0.46 during arousal). To identify common modes of activity across regions in a data-driven way, we performed a principal components analysis to extract the principal components of all thalamic and cortical regions during arousal. We found that most (98.25%) variance could be explained with just two components (Fig. 3d). The first component was driven primarily by cortical regions and exhibited a large decrease after arousal, and the second component was driven primarily by thalamic regions, with activity increasing prior to arousal (Fig. 3e). These results demonstrated that thalamic and cortical regions largely segregate into opposite dynamics during behavioral arousal. 7T imaging shows two primary activity modes at arousal, largely segregated into thalamic activation and cortical deactivation. a) Schematic of thalamic nuclei of interest that were resolved with this imaging protocol. b) In one subject scanned with simultaneous EEG, mean alpha significantly increases during behavioral arousal (dashed line), consistent with Experiment 1. c) Experiment 2 replicates the result that the thalamus activates (latency indicated by red arrow) prior to behavioral arousal, while cortical deactivation (blue arrow) follows (n=13 subjects, 99 arousals). d) A principal component analysis of all 9 thalamic regions and 30 cortical regions reveals two primary modes of activity in the thalamus and cortex: one which increases before arousal, and one which decreases after arousal. e) The first principal component is more heavily influenced by cortical regions, and the second is influenced primarily by thalamic regions, demonstrating a segregation of thalamic and cortical activity during behavioral arousal.

An activity sequence across thalamic nuclei occurs prior to arousal
While fMRI has previously identified spatial maps of regions implicated in arousal state transitions 50 , we used the high temporal resolution of our imaging to test whether a temporal ordering of activity appeared across the thalamus during behavioral arousal. Specifically, does activity in one nucleus foreshadow the subsequent thalamic-wide activation, or does the thalamus activate as a whole during arousal transitions? To compare the temporal properties of individual thalamic nuclei, we calculated the mean arousal-locked fMRI signal across time (Fig. S4). To quantify relative timing across the thalamus, we implemented a cross-correlation analysis which calculated the lag between each thalamic nucleus and the mean signal from the whole thalamus. We found that the arousal-locked fMRI signal across thalamic nuclei exhibited highly heterogeneous timing properties, spanning a range of 2.1 s (Fig. 4a). The centromedian (CM) and ventral posterolateral (VPL) nuclei consistently led the rest of the thalamus, whereas the ventral lateral anterior (VLA), ventral anterior (VA), and anteroventral (AV) nuclei lagged (Fig. 4a). Thus, a subset of spatially distinct thalamic nuclei consistently showed early activity, with subsequent activity unfolding across the remaining nuclei prior to arousal (Fig. 4b).
We next investigated whether this thalamic heterogeneity reflected a temporal sequence specifically in the onset time of activity across nuclei, since onset time may be more reliable than peak amplitude time in determining lags between regions by using fMRI 62 . We fit a double-Gaussian model to estimate a smooth hemodynamic response function for each nucleus, which allowed us to define onset time (Fig. 4c, S5). We found that CM indeed activated first, 0.68 s before the next fastest nucleus, the VPL, and 4.08 s before the last nucleus, VA (Fig. 4c, d). This result confirmed a specific temporal sequence of arousal-locked thalamic activity underlies The relative timing of thalamic activation differs significantly across nuclei (mean lag shown by the solid vertical line; shaded region is 95% confidence interval, colored from red to purple based on order in sequence). Zero represents the timing of the whole thalamus (dashed line). The centromedian (CM, red) and ventral posterolateral (VPL, orange) thalamic nuclei lead the rest of the thalamus during behavioral arousal, while the ventral anterior (VA, light blue), anteroventral (AV, dark blue), and ventral lateral anterior (VLA, purple) lag behind. b) Image of mean lag times per nucleus in one subject. Color represents precise lag. c) The activation onset times (arrows) and the hemodynamic response (solid lines) show a sequence in activation across nuclei. Color represents order in lag sequence. d) Thalamic nuclei activation onset is sequenced similarly to the lag values. the transition between behavioral arousal states, with CM activity appearing first, VPL next, and these earliest nuclei activating seconds prior to the rest of the thalamus.

Hemodynamic lags and systemic physiological changes cannot explain thalamic latency differences
Since blood oxygenation level-dependent (BOLD) fMRI signals arise from hemodynamic responses to neural activity, and this hemodynamic response can be heterogeneous across the brain 60,70-73 , we next examined whether differences in vascular reactivity existed across thalamic nuclei that could contribute to these temporal sequences. As a control experiment, we measured fMRI signals caused by the widespread vascular response to a brief breathhold in order to assess hemodynamic latency differences across nuclei [74][75][76][77][78][78][79][80] . We found that the breathhold task produced robust fMRI signals across thalamic nuclei (Fig. S6a). The amplitude of breathlocked responses across thalamic nuclei was strongly correlated with the previously observed arousal-locked signal amplitudes (r=0.88; Fig. S6b), demonstrating that the breathhold task provided a reliable control experiment to investigate hemodynamic properties of thalamic nuclei. Therefore, we repeated our temporal sequence analysis in these control data. We found that the thalamus led the cortex by only −0.12 s (95% CI: (−0.43 s, −0.062 s)) during the breathhold task (Fig. 5a), far less than the 5.3 s latency difference observed at arousal, and consistent with prior studies showing slightly faster hemodynamics in thalamus than in cortex 60,81,82 . Furthermore, while minor differences in the breathhold fMRI response were detected across thalamic nuclei, the lag differences between nuclei could not account for the differences observed during arousal (Fig. 5b). One nucleus, the AV, lagged significantly behind the whole thalamus by 0.99 s (CI=(0.12 s,1.3 s)), suggesting that the late arousal-locked activation of AV was partly due to a slower hemodynamic response rather than delayed neural activation. The mediodorsal (MD) nucleus had slightly slower, the VPL and VLA nuclei exhibited slightly faster (less than 0.25 s in each nucleus), and the CM, pulvinar (PUL), lateral geniculate (LGN), ventral lateral posterior (VLP), and VA nuclei showed similar vascular reactivity to the whole thalamus. Altogether, these differences across nuclei were small, and could not explain the multi-second thalamic sequence observed at behavioral arousals. Furthermore, we corrected our arousal sequence for each nucleus' hemodynamic latency and still found equivalent groupings of early and late nuclei (Fig. S6c). These results supported the conclusion that a temporal sequence of neural activity across thalamic nuclei occurs prior to the moment of arousal, with the earliest activation first in CM and then in VPL.
A second modulator of fMRI signals is systemic physiology: systemic physiology changes across arousal states, and in particular, changes in heart rate and respiration accompany arousal [83][84][85] . Such shifts in systemic physiology can also contribute to the fMRI signal 74,[86][87][88][89] . To investigate how these systemic physiological dynamics were linked to behavioral arousal, we calculated the average respiratory and fingertip pulse-oximeter amplitude. We indeed found that changes in the respiratory (Fig. 5c) and pulse-oximeter (Fig. 5d) signal amplitudes were locked to behavioral arousal, consistent with known physiological modulation linked to arousal. However, both the respiratory and pulse-oximeter signals changed after behavioral arousal. These systemic physiological effects therefore could not produce hemodynamic confounds to explain the thalamic activity beginning seconds earlier. Altogether, these results demonstrated that shifts in systemic physiology accompany the transition in behavioral state, consistent with the expected physiological changes during shifts in arousal, but that they occurred after the arousal-locked thalamic activity sequence began. We thus concluded that these thalamic temporal dynamics represented a neural activity sequence that began to unfold across thalamic nuclei before the moment of behavioral arousal.

Thalamic signals, but not cortical signals, differ in sustained and transient behavioral arousals
Arousal can vary in length-the subject may stay awake and continue to perform the task, or may fall directly back into non-behaving sleep after a brief arousal. Our study thus provided an opportunity to ask whether these thalamic dynamics reflected simply the moment of transition between states or also the maintenance of a subsequent higher arousal state. We separated the behavioral arousals into two types (Fig. 6a): sustained arousals (at least five responses following arousal) and transient arousals (two or fewer responses), resulting in 45 sustained arousals and 34 transient arousals. Surprisingly, even in brief arousals consisting of two or fewer button presses, we found strong activity increases in the thalamus (Fig. 6b). However, the thalamic activity plateaued at a level higher than its pre-arousal baseline during sustained arousals, whereas after transient arousal, the thalamus deactivated and returned to baseline (Fig. 6b). Interestingly, the mean cortical signal did not have a notable difference between fMRI signatures in transient versus sustained arousals (Fig. 6c). Therefore, a stable plateau of thalamic activity reflected differences in subsequent maintenance of arousal.
We next evaluated dynamics across individual thalamic nuclei in sustained and transient arousals to test whether they were linked to the subsequent maintenance of arousal. We found that most nuclei exhibited higher post-arousal plateau activity in sustained arousals, similar to the whole thalamus (Fig. S7). The LGN demonstrated the smallest difference between post-arousal plateaus in sustained and transient arousals, whereas all other nuclei showed striking significant differences (Fig. S7). Next, to determine if the temporal sequence differed across arousal types, we repeated our cross-correlation lag analysis. We found that in both sustained and transient arousals, there was an activity sequence across thalamic nuclei, beginning with the CM Figure 5. Hemodynamic latencies cannot explain the arousal-locked temporal dynamics in the thalamus. a) fMRI signals increase after breathhold release (purple dashed line). The time of the breathhold is shaded in light purple. The thalamus leads cortex by −0.12 s. b) Small hemodynamic lags between thalamic nuclei exist, but are not large enough to account for the arousal-locked sequence. The color and ordering of thalamic nuclei represent the lag sequence during behavioral arousal. Mean lag (solid line) and 95% confidence interval (shading) calculated by same method as Fig. 4d. Purple dashed line represents a zero lag relative to the whole thalamus. c) Respiration amplitude increases at behavioral arousal (dashed line). d) Amplitude of pulse-oximeter signal decreases after behavioral arousal. Shading is standard error. (Fig. S8). However, thalamic nuclei activated more closely together in time during sustained arousals with a range of 1.17 s (Fig. S8a). During transient arousals, the activity sequence was more distributed with a range of 2.46 s (Fig. S8b). Therefore, during sustained switches in behavioral state, the temporal activity sequence across nuclei concluded in coordinated activity throughout most of the thalamus, whereas transient arousal was accompanied by a more sluggish sequence across nuclei. Thus, network activity across thalamic nuclei reflected whether the transition in behavioral arousal state was subsequently maintained.

DISCUSSION
We conclude that behavioral arousal is preceded by a temporal sequence of activity across thalamic nuclei, followed by widespread cortical deactivation. A consistent temporal pattern emerged across individual thalamic nuclei at the moment of arousal, and this thalamic activity was linked to subsequent maintenance of behavioral arousal. Our results are consistent with prior animal [18][19][20][21][22] and human fMRI studies 49,50,66 showing that thalamic activity is coupled to arousal state. Our work further demonstrates that thalamocortical networks exhibit highly structured spatio-temporal dynamics that evolve throughout this transformation between unresponsiveness and active behavior, with distinct profiles for distinct thalamic nuclei.
We consistently found that an intralaminar thalamic nucleus-specifically CM-activated first within this thalamic temporal sequence. This early activity in CM is consistent with studies in mice that recorded from intralaminar nuclei during arousals. During NREM to wake transitions in mice, the central medial nucleus (CeM), an intralaminar thalamic nucleus, activates earlier than the ventrobasal nucleus, a sensory thalamic nucleus, and optogenetically driving CeM induces awakening 34 . CeM efferents terminate on cingulate and insular cortices, and optogenetic activation of these axon terminals induces awakening 34 . The temporal sequence we observed is consistent with these dynamics, with the CM, an intralaminar nucleus that receives input from the reticular activating system as well as sensorimotor information 90,91 , activating earliest. Furthermore, electrical stimulation of the central thalamus causes behavioral arousal transitions in anesthetized macaques 32,39 , and deep brain stimulation of the central thalamus modulates behavioral responses in minimally conscious state patients 92 . Our results are thus consistent with invasive studies showing that intralaminar thalamic nuclei can drive changes in arousal state and further demonstrate that in spontaneous arousals, activity originates in an intralaminar nucleus and subsequently spreads to other thalamic nuclei.
Our results suggest potential mechanisms that could generate such a thalamic temporal sequence. The early activation of the CM and VPL nuclei indicates that there may be a sequential chain, with activity in CM spreading to neighboring nuclei either directly or via thalamocortical loops. Arousal-driving thalamic nuclei may shift baseline excitation in the other nuclei of the thalamus, setting the stage for subsequent thalamic activation and arousal transition. In this case, the lagged dynamics could reflect the dispersion of activity following the lead of intralaminar nuclei such as CM. However, driving the activation of any subcomponent of the chain would not necessarily trigger the rest of the sequence, since optogenetic stimulation of lower-order nuclei does not cause arousal transitions 34,39 . If only a subset of nuclei, such as the intralaminar nuclei, can effectively raise baseline excitation in neighboring nuclei and induce broader thalamic synchronization, this may explain why stimulating subcomponents of the chain does not cause arousal transitions.
We found that while activity in most of cortex declined at arousal, the ACC and PCC notably diverged and instead mirrored the thalamus, with a large increase in activity before arousal. Prior studies have shown that these regions are linked to state-switching and attentional modulation within the awake state: the ACC is involved in sustained arousal in complex, novel environments 102 and attentional-state switching 103 , and the dorsal ACC is active during moments of surprise 104 . The PCC supports internally directed cognition 105 , regulation of the focus of attention 106,107 , and is a component of the default mode network, which is more active in the task-negative state 108,109 . Our results align with this role for cingulate cortex in attentional state switching, and further suggest that these cingulate regions may also have a broad, non-specific role in task engagement, as they are linked to effective thalamic modulation of behavioral arousal state, and the transition from complete unresponsiveness to active behavior.
Our finding that most cortical regions deactivated after the moment of behavioral arousal raises several possible interpretations, as both neural mechanisms as well as indirect neuromodulation of vasculature could potentially contribute to this decrease in cortical BOLD signal 93 . Neural inhibition could cause a decrease in BOLD 94 ; however, the spike rate of neurons in the cortex is typically higher during wakefulness than during sleep 95,96 , suggesting the decline in cortical signal may have different origins. One possibility is that the shift from coherent to desynchronized cortical activity that occurs with arousal may drive this decrease, since changes in neural synchrony modulate the BOLD signal 97 . Additionally, the decrease in cortical BOLD could also reflect vascular factors. In some circumstances, it is possible for cerebral blood flow to decrease despite increased neural activity 98 . This can be caused by a redistribution of oxygenated blood to more active areas, such as deep brain regions 99 , or endogenous alterations in vasoconstriction due to neuromodulators such as acetylcholine and noradrenaline 100,101 . Notably, this cortical decrease may also be partly a result of systemic vascular changes, as arousal is also followed by systemic vasoconstriction 86,88,89 . A combination of hemodynamic and neural mechanisms may thus contribute to the decrease in cortical BOLD signal after arousal.
Since vascular dynamics can influence the fMRI signal and are strongly modulated by arousal state, we performed control analyses to confirm that the activity sequence observed across thalamic nuclei was not driven by hemodynamic confounds. Our breathhold experiment tested for local differences in cerebrovascular reactivity, as while the hemodynamic response to a myogenic breathhold task may not exactly replicate the neurogenic vascular response 77 , it can be used to correct for hemodynamic lags across the brain 75 . Correcting for hemodynamic lags demonstrated that most thalamic nuclei had only very small differences in vascular reactivity, with the exception of a notably slow response in AV. Exactly why AV has a slower hemodynamic lag is not clear, but could reflect differences in its vascular architecture 71 . In addition, we found that systemic physiological changes in heartbeat and respiration occur after arousal, and therefore could not explain the early thalamic activation prior to arousal. Our observations of temporally patterned thalamocortical dynamics were therefore robust even after accounting for vascular physiology.
Our observation of an overall thalamic activation and cortical deactivation at arousal replicates the main findings of prior fMRI studies of awakening from sleep 50 , but with important differences. In this prior work, the thalamic and cortical dynamics were reported to be synchronous, whereas our study found that the thalamus activates before the cortex deactivates, replicated in our two separate experiments. Our detection of precise temporal sequences can be primarily attributed to the fast temporal sampling of our acquisition, which enhances sensitivity to timing differences 110,111 . However, a second potential contributor to this difference is that prior work aimed to reduce physiological noise by regressing the cerebrospinal fluid (CSF) signal, a common approach in fMRI. However, CSF signals are strongly coupled to neural activity and the fMRI signal during sleep 65 . Since CSF signals can be highly collinear with neural signals during sleep, regressing out the CSF could remove neurally relevant dynamics and distort the timing of the underlying signals 112 . This issue reflects a broader challenge in neuroimaging studies of sleep and arousal: the neural dynamics that drive arousal state changes are often collinear with the systemic physiological changes they induce, and common noise removal practices can have unintended effects.
Future work could investigate how these thalamic sequences are related to distinct types of arousal transitions. Here, we studied spontaneous behavioral arousal, but there are many ways to awaken from sleep and different sleep stages from which arousal can occur. An intracranial study in humans found that arousallocked thalamic spectral signatures are highly stereotyped across different sleep stages and in stimulus-driven arousals 113 . However, the precise timing of these activity patterns has not yet been measured, and the order of activation across thalamic nuclei may differ in arousals driven by sensory stimuli. For example, one might expect the visual nuclei of the thalamus (LGN and PUL) to activate earlier during arousal driven by visual stimulation, such as the lights turning on in a darkened room. Additionally, arousal transitions can result in different subsequent awake states, such as waking up afraid, which may also correspond to distinct thalamocortical dynamics. Future fast fMRI studies investigating stimulus-driven arousals could potentially uncover different sequences of activation associated with such arousals. In addition, arousal itself can be defined through multiple metrics. Importantly, here we focused on behavioral arousals, but thalamocortical network activity may differ during electrophysiological arousals, such as those defined by the return of occipital alpha in the EEG. In particular, we would predict that the LGN may deactivate during the return of occipital alpha rather than activate as we report in behavioral arousals, due to the established anticorrelation between LGN activity and alpha rhythms 114 . Future studies could use this accelerated functional imaging approach to delineate the thalamic activity linked to electrophysiological, as compared to behavioral, arousal transitions.
In conclusion, we identified a sequence of activity that unfolds across thalamic nuclei during behavioral arousal state transitions. Our results reinforce the role of the thalamus in arousal state transitions and identify distinct functional profiles for the diverse nuclei within the thalamus, uncovering precise temporal patterns that are coupled to both the moment of transition and the stability of the subsequent arousal state. Finally, these results suggest a broad potential for fast, ultra-high field fMRI to identify the temporal dynamics of subcortical activity in the human brain. As subcortical structures are challenging regions to image but play fundamental roles in cognition and awareness, this highlights a valuable approach for studying their functional roles at sub-second timescales.

Acknowledgments
We are grateful to Sydney Bailes, Dr. Thomas Witzel, and Dr. Donald Straney for their contributions to imaging and technical help. This study was funded by National Institutes of Health Grants R00-MH111748, R01-EB019437, P41-EB300006, S10-OD010759, the NARSAD Young Investigator Award from the Brain and Behavior Research Foundation, the Searle Scholars Program, the One Mind Rising Star Award, the Alfred P. Sloan Fellowship, the Pew Biomedical Scholar Award, and a Hariri Institute for Computing Fellowship to B.S.

Participants
The current study includes three datasets. All experimental procedures were approved by the Massachusetts General Hospital Institutional Review Board. Participants in all datasets were screened not to have any neurological, psychiatric, or sleep disorders and not currently be taking psychiatric or sleep medications. The first dataset (Experiment 1) used 3T fMRI data from a previously published study of cerebrospinal fluid dynamics during sleep 65 . Our current analysis of these data included the subset of subjects who were instructed to perform a behavioral task during the sleep imaging (n=6; one male and five female), with a mean age of 24.6 years (range: [23][24][25][26]. The second dataset (Experiment 2) was a newly acquired imaging experiment using the same behavioral task, performed with fast fMRI at 7T. Written informed consent was obtained from 20 healthy adults (14 female and six male; mean age: 24.9, age range: [22][23][24][25][26][27][28][29][30][31][32][33]. Seven of these subjects were scanned with simultaneous EEG. Participants were excluded from all analyses if they lacked any behavioral arousals (n=5) or excessive motion artifacts (n=2). In total, 13 subjects were included in the analysis (nine females and four males, mean age= 25.2, age range: 22-33), four of whom had simultaneous EEG. Subjects who did not present occipital alpha were excluded from the EEG analysis portion only, which resulted in one subject with simultaneous 7T fMRI and EEG being included in our EEG analysis.
For the third dataset (Experiment 3), eight subjects participated in the 7T fMRI breathhold control experiment (three females and five males, mean age=25.5, range: 20-34), with the same exclusion criteria as Experiment 2. Three of these individuals were also participants in Experiment 2. No participants were excluded.

Sleep and behavioral protocols -Experiments 1 and 2
Participants were asked to restrict their sleep to four hours the night prior to the sleep scan (Experiments 1 and 2) and avoid caffeine 24 hours prior to the study. Subjects self-reported the number of hours they slept. Imaging sessions began between 10:30 pm and midnight, the night after sleep restriction.
Participants in Experiments 1 and 2 were instructed to keep their eyes closed for the scan duration. Participants were instructed to press a button on every breath in and out as long as they were awake and that it was permitted to fall asleep during the scan. Behavioral arousals were defined as the first response to this behavioral task after at least 20 s of unresponsiveness. Arousals that included excessive motion (>0.3 mm) were excluded from the analysis.

BOLD fMRI acquisition
Participants in Experiment 1 lay supine in an MRI scanner (Siemens 3T Prisma, 64-channel head-andneck coil). First, we performed a T1-weighted anatomic multi-echo MPRAGE scan with 1-mm 3 isotropic voxels to provide an anatomical reference 115 . Participants then underwent BOLD-weighted echo-planar-imaging (EPI) scans for fMRI data acquisition covering most of the brain, containing up to 1500 volumes and lasting up to 90 minutes per run. Single-shot gradient-echo SMS-EPI 56 data were acquired with a voxel size of 2.5 mm 3 isotropic, MultiBand factor=8, matrix=92×92, shift factor=4, TR=367 ms, TE=30 ms, flip angle=32-37°, echo spacing=0.53 ms, and no in-plane acceleration.
Participants in Experiment 2 lay supine in an MRI scanner (Siemens 7T Magnetom, custom-built 32channel head coil). First, a T1-weighted anatomic multi-echo MPRAGE scan with 0.75 mm 3 isotropic voxels was performed to provide an anatomical reference 115,116 . Participants then underwent one to three EPI scans for fMRI data acquisition, containing up to 8000 volumes and lasting up to 33 minutes. 40 oblique slices were acquired with a voxel size of 2.5 mm 3 isotropic, TE=24 ms, TR=247 ms, MultiBand factor=8, shift factor=4 matrix=84×84, flip angle=30°, echo spacing=0.53 ms, and no in-plane acceleration.
Participants in Experiment 3 underwent the same imaging protocol as Experiment 2, but with functional runs split into 8-minute periods. Each participant underwent four to eight functional runs.
Fast fMRI is susceptible to physiological noise 87,119 , including cardiac rhythms and respiration, which appear in the data as quasiperiodic BOLD oscillations. We used dynamic regression based on the concept of RETROICOR 120 to remove signals driven by the heartbeat and respiration from the data while allowing the peak frequency of these physiological rhythms to vary over time. The cardiac signal was bandpass filtered between 0.2 and 10 Hz. Peaks of the cardiac signal were identified using the automated peak detection technique in the Chronux toolbox (http://chronux.org/) 121 , and the interpeak intervals were transformed into phases. The respiratory signal was bandpass filtered between 0.16-0.4 Hz using a finite impulse response filter, and the instantaneous phase was computed as the angle of the Hilbert transform. This phase information was transformed into sine functions, and beta values were estimated in a window of 1000 s sliding every 400 s voxelwise using a general linear model. These values were then interpolated across each time-point and used to remove each voxel's first and second harmonic frequencies of the cardiac and respiratory signals.
Spatial smoothing was not applied to avoid blurring and to maximize spatial accuracy within the thalamus. All fMRI signal analyses were performed in the original spatial frame of the fMRI acquisition space for each individual without transforming to a common average. Registration between the functional and anatomical images was completed using Freesurfer boundary-based registration 122 . Cortical ROIs were extracted using the Desikan-Killiany atlas 67 to identify functional voxels that were at least 70% filled by each region. The thalamic segmentation was done using the individual-level probabilistic atlas in the Freesurfer developmental version which provides voxelwise segmentation probabilities in individual anatomical space, that have been previously validated against ex vivo histology 68 . Nuclei were defined by selecting functional voxels that had at least a 90% chance of falling within a given thalamic nucleus to minimize partial-volume effects. All voxels that fell in any of the four pulvinar sub-regions were combined into a single pulvinar nucleus ROI, and similarly, and any voxels than fell in either of the two mediodorsal sub-regions were combined into the mediodorsal nucleus ROI. Cortical and thalamic regions which were not captured in every subject were excluded from the analysis. The mean timecourse in each cortical and thalamic ROI was then extracted using FSL version 6.

Physiological recordings and analysis
In subjects with simultaneous EEG, physiological signals were simultaneously recorded using a Physio16 device (Electrical Geodesics, Inc., Eugene, OR USA). ECG was measured through two disposable electrodes placed on the chest diagonally across the heart, with an MR-compatible lead (InVivo Corp, Philips). Respiration was measured through a piezo-electrical belt (UFI systems, Morro Bay, CA, USA) around the chest.
In the remaining subjects in Experiments 2 and 3, pulse was measured with a pulse-oximeter around the left thumb (UFI, Morro Bay, CA, USA). Respiration was measured through a piezo-electrical respiration transducer belt around the chest (UFI, Morro Bay, CA, USA).
To investigate how respiration and heart rate change during arousal, we calculated the mean amplitude of the respiratory and cardiac signals. The respiratory signal was bandpass filtered between 0.16-0.4Hz, and the pulse-oximeter signal was bandpass filtered between 0.7-1.8 Hz. Then, the instantaneous amplitude was calculated using the absolute value of the Hilbert transform and was averaged across all arousals.

EEG recordings, preprocessing, and analysis
In Experiment 1, EEG was measured by an MRI-compatible 256-channel geodesic net and a NA410 amplifier (Electrical Geodesics, Inc., Eugene, OR USA) at a sampling rate of 1000 Hz. Each subject wore a reference layer cap composed of an isolating vinyl layer and conductive satin layer on the head, with grommets inserted to allow electrodes to pass through and make contact with the scalp 123 , while other electrodes remained isolated from the scalp and recorded the noise, resulting in a total of 30-36 EEG electrodes per subject. The MR scanner's helium pump was temporarily shut off during EEG acquisition to reduce vibrational artifacts in the EEG signal. In a subset of the analyzed subjects (n=4) in Experiment 2, EEG was measured by a specially designed 256-channel polymer thick film MR-compatible EEG net, called the InkNet, for EEG collection at 7T 69 . A custom reference cap was placed under the InkNet to measure ballistocardiogram (BCG) noise 65,123 . EEG acquisition was synchronized to the scanner 10 MHz master clock to reduce the aliasing of high-frequency gradient and RFinduced artifacts. A baseline of eyes open versus eyes closed, eye movement, and jaw clenching were recorded outside of the scanner.
Gradient artifacts were cleaned from the EEG data using average artifact subtraction 124 . EEG data were filtered using a linear-phase filter between 0-50 Hz and then down-sampled to 200 Hz. The EEG channels and BCG channels were re-referenced to their respective average. All channels on the cheeks were excluded from this average. The BCG artifact was then removed by adaptively regressing the influence of the BCG channels average signal from each EEG channel over a sliding window of 30 s [123][124][125] .
EEG was used to determine the sleep stage prior to arousal and if these behavioral arousals correlated with a shift in electrophysiological arousal. The EEG data in Experiment 1 were sleep-scored in 30-second intervals using the standard vigilance state scoring criteria defined by the AASM 2 . The sleep stage in the 30second window 15 s prior to arousal was recorded.
The spectrogram of the cleaned EEG data was computed in order to identify shifts in electrophysiological arousal. The alpha power (7-12 Hz) was calculated in the electrode closest to the Oz position. Power in the alpha frequency range was calculated using Chronux 126 (tapers=3, sliding window length=2, sliding step= 1]). To test for a significant rise in alpha power during behavioral arousal, we compared the average occipital alpha power in the 10 s before behavioral arousal versus the 10 s after. We averaged the occipital alpha power in these time windows for each arousal and calculated if there was a significant change in alpha power using a two-sided, paired t-test. A p-value of less than 0.05 was considered significant.

Arousal-locked signal analysis
After extracting the ROI timeseries, we interpolated the fMRI signal to four times its original sampling rate (61.75 ms), and identified the sample closest to the arousal. We extracted the fMRI signal 10 s before and 20 s after this sample. Arousal windows that had a framewise displacement at any point over 0.3 mm were excluded from analyses. The resulting signals with low motion were averaged to calculate the mean fMRI response during arousal for each ROI. These means were centered at the average value during the first 3 s of the time-series (−10 s, −7 s).
To determine the latency of the thalamus and cortex, we identified the time at which the fMRI signal reached 10% of its maximum absolute amplitude. The baseline of the mean fMRI signal was calculated as the mean BOLD signal in the −10 to −7 seconds before arousal. We used a hierarchical bootstrap analysis to calculate the 95% confidence interval of the thalamic and cortical latencies. First, we resampled from the subjects with replacement, and then we resampled the arousals with replacement to re-calculate the mean fMRI timecourse and latency of each region. We did this 1000 times. The 95% confidence interval was calculated by using the 2.5 and 97.5 percentiles.
To test for a significant increase in any cortical region's fMRI signal during behavioral arousal, we used a time-binned, two-sided t-test. We averaged each ROIs fMRI signal across 1-second bins. Each bin was tested using a two-sided t-test to determine if the fMRI signal was significantly above zero. The p-value was Bonferroni corrected for the number of time bins and cortical regions tested, and time-bins with corrected p-values less than an alpha value of 0.05 were considered significant.
A principal component analysis was used to determine the major activation patterns that underlie arousal state transitions in all of the cortical and thalamic ROIs. The mean timeseries of each ROI 10 seconds before and 20 seconds after arousal was centered to the first 3 s of the timeseries (-10,-7), and singular value decomposition was used.
The estimation of lags between ROIs during arousal was possible due to the high temporal resolution of our fMRI data. A cross-correlation analysis was used to compute the time lag that produced the maximal correlation between the mean signals of each thalamic nucleus and the whole thalamus. A hierarchical bootstrap analysis was used to resample with replacement from the arousals, generate a resampled mean timeseries, and calculate 95% confidence intervals of the lags between regions, resampling 1000 times. First, we resampled the subjects with replacement and then resampled arousals with replacement from individual subjects. The mean fMRI signal during arousal for each ROI was computed, and the cross-correlation analysis was repeated. The upper and lower bounds of the 95% confidence interval of the lag between each nucleus and the whole thalamus were computed by calculating the 2.5 and 97.5 percentiles.
Next, we determined the onset time of each thalamic ROI. Since these small, deep regions can be noisy, we fit a linear combination of two Gaussian curves to each ROI using a simple search method 127 to find the parameters with the smallest root mean squared error. We defined onset time as the time that the model fit reached 10% of its maximum amplitude. In order to allow for temporal flexibility, we included the temporal derivative of each Gaussian in the model fit.
In order to determine if thalamic and cortical activity deviated between arousals in which the subject stayed awake after the initial button-press or fell directly back asleep, we segregated the behavioral arousals into arousals with five or more button-presses (sustained arousals) and two or fewer button-presses (transient arousals). Significant differences between ROI's fMRI signals in sustained and transient arousals were identified by computing a time-binned, paired t-test. The fMRI signal was averaged across 1-s bins, and the p-value was Bonferroni corrected for the number of time bins and ROIs tested. We then repeated the cross-correlation, and bootstrap analyses in individual thalamic nuclei as described above.

Breathhold cerebrovascular reactivity measurements -Experiment 3
We used a breathhold task to measure fMRI signals caused by a vascular response to test whether the fMRI sequences observed during arousal state changes are due to hemodynamic latency differences 74,75,86 . When subjects hold their breath, the brain compensates for the reduction in oxygen by increasing blood flow to the brain, and when the breathhold is released, there is a momentary increase in fMRI signal due to the increased vasodilation. The dynamics during this vascular response were compared to the dynamics during arousal to test whether results were confounded by inherent differences in vascular reactivity.
A projector displayed instructions for subjects to breathe freely for 27 s, followed by three paced breaths of 3 s in and out, then to hold their breath for 15 s. This cycle was repeated eight times per run. 4-8 runs were collected per subject. Breathhold releases were confirmed using simultaneously recorded respiratory signals. These scans underwent the same preprocessing steps as the sleep scans: motion correction, slice time correction, and physiological noise removal. Breathholds with excessive motion (>0.3mm) were excluded. Timeseries were extracted 20 s before and 20 s after the subject released their breathhold and began to breathe freely in order to capture both the beginning and end of the breathhold. ROI mean signals were averaged over this period, and the cross-correlation between regions was calculated. The hierarchical bootstrap analysis was repeated to generate 95% confidence intervals. These results were compared to the mean fMRI signal of ROIs during arousal state transitions. To correct for the hemodynamic latencies based on the breathhold response, we subtracted the average breathhold lag from the behavioral arousal lag of each nucleus.