Brain-wide resting-state fMRI network dynamics elicited by activation of single thalamic input

Domain: neuroscience
Published: 2025-12-17
Editorial outcome: major_revision

Abstract
Understanding how the brain rapidly processes neural information through large-scale network dynamics is fundamental to neuroscience. Resting-state fMRI (rsfMRI) captures networks by probing spontaneous neural activity fluctuations, yet whether these networks can accommodate fast, sensitive information processing remains unknown. Here, we reveal brain-wide rsfMRI network dynamics (i.e., states and their transitions) that enable fast information processing through rapid activity-triggered network reconfiguration in male rats. Following activation of a single optogenetically evoked somatosensory thalamic input, states and transitions associated with internally-oriented information processing are quickly suppressed and diverted through basal forebrain and hypothalamus to strengthen externally-oriented (i.e., sensory) processes. Simultaneous rsfMRI and electrophysiology measurements in somatosensory and cingulate cortices demonstrate that intra-regional wide-band and inter-regional narrow-band (i.e., slow oscillations) neural activity synchronizations directly drive rsfMRI network dynamics and reconfiguration to support fast processing of single thalamic input. Our study uncovers fast interactions between neural activity and rsfMRI networks in processing single neural input.

Introduction
The brain is a complex system composed of diverse interacting neural populations, circuits, and networks that are essential for everyday functions and behaviors1. In recent decades, resting-state functional MRI (rsfMRI) has emerged as a robust, non-invasive imaging technique for mapping brain-wide functional networks across species in both humans and rodents2–7. Notably, such rsfMRI measurable networks are highly correlated with diverse sensory, behavioral, and cognitive task performance8,9, and can be used to phenotype brain dysfunction10,11. Furthermore, rsfMRI networks are increasingly recognized for exhibiting dynamic processes driven by spontaneous neural activity fluctuations12–15, which are crucial for neural information flow and processing16–20. Dissecting these large-scale, time-varying networks during both resting and task execution can help unravel the brain’s functional architecture and its role in neural information processing.

Brain network communication primarily relies on the flow of activity among neuronal populations through numerous pathways21. This process is intrinsically dynamic because the brain must continuously adapt the behaviors to environmental changes, switching between brain states that process externally- and internally-oriented information22–24. Externally-oriented processing is driven by stimuli and task demands22, while internally-oriented processing relies on intrinsic cognitive processes such as memory retrieval23,24. Recent studies have explored network characteristics or configurations, such as network connectivity states and their transitions, which underlie the dynamics of large-scale rsfMRI networks25–27. These rsfMRI network states and transitions closely parallel the brain states and the switching within and between these brain states, which are characterized by patterns of neural activity associated with physiological or cognitive processes28,29, like wake-NREM sleep stages30 and shifts in cognitive task demands31,32. Studies have shown that rsfMRI networks can adjust their configurations to accommodate alternating task demands and behavioral contexts17,32–35. For instance, rsfMRI network reconfiguration supports the selective rerouting of visual information in sensory systems34 and the access of learned information from cognitive systems to aid decision-making33. However, in these sensory and cognitive tasks, external stimuli—such as an oriented bar in the visual field36, whisker deflection37, an odor38, or an auditory tone39 evoke a cluster of complex and long-lasting neural activities that intermingle and cascade with neural information flows from multiple parallel pathways, making it challenging to investigate whether and how rsfMRI network dynamics rapidly process specific neural activity. Additionally, in contrast to localized sensory-evoked or optogenetically-evoked hemodynamic responses (with onset delays as short as 400 milliseconds40,41), the change in rsfMRI network configurations is often perceived as slow. This is primarily due to the low-frequency fluctuations (0.005–0.1 Hz) of blood-oxygenation-level-dependent (BOLD) rsfMRI signals and the nature of the tasks or stimuli. In other words, while it is known that rsfMRI network configuration can probe the dynamics of brain-wide functional networks, it remains unclear and challenging to determine whether rsfMRI networks can support fast neural information flow and processing by sensitively accommodating relatively simple and short-lasting neural activity, as well as to understand the mechanisms underlying these processes.

In this study, we investigated how brain-wide rsfMRI networks enable fast neural information flow and processing. We examined the effects of activating single thalamic input on rsfMRI network dynamics using single-pulse optogenetic stimulation of the somatosensory thalamus. Such stimulation allowed us to bypass parallel sensory pathways typically engaged by external stimuli, and to avoid the complex cascades and excessive synchronization of neural populations caused by prolonged or repetitive stimulation. By using a single, brief neural population activity input, we were able to capture activity-triggered rapid reconfiguration of rsfMRI networks while minimizing disruption of ongoing spontaneous neural activities. Using unsupervised machine learning analysis on rsfMRI data obtained with and without optogenetic stimulation, we identified key elements of brain-wide rsfMRI network dynamics, including transient network states and their hierarchical transitions, with spatiotemporally specific evolution patterns both within and between states. We demonstrated that these network states and transitions could rapidly reconfigure on the order of a few seconds (~2 s to begin reconfiguring; ~10 s to normalize) in response to a single neural population activity input from the somatosensory thalamus, suppressing internally-oriented processing activities and redirecting them to externally-oriented sensorimotor information processing. Combining simultaneous electrophysiology and fMRI measurements, we further revealed that this reconfiguration was driven by distinct synchronizations of neural activity within and between networks, occurring across specific time and frequency windows.

Results
To investigate the rapid dynamics of brain-wide rsfMRI networks and their responses to activation of a single brief thalamic input, we conducted interleaved rsfMRI measurements both with and without sparse single-pulse (10 ms, 0.066% duty cycle) optogenetic stimulation of the ventral posteromedial thalamic nucleus (VPM) thalamocortical excitatory neurons42,43 in lightly anesthetized rats (~1% isoflurane; Fig. 1A, B). We identified typical rsfMRI networks using seed-based analysis5,43,44 and derived elements of brain-wide rsfMRI network dynamics (i.e., states and their transitions) through Gaussian Hidden Markov Modeling45 (HMM; Fig. 1C). To further analyze the spatio-temporal evolutions of rsfMRI network states, each state was decomposed into three substates based on its probability time course. Voxel-wise activation maps were generated from pre-processed BOLD rsfMRI time courses aligned with the corresponding substates. By examining the evoked neural activity propagation structure and neural activity variabilities upon single-pulse optogenetic stimulation, we first demonstrated that this stimulation paradigm served as an effective tool for studying single thalamic activity-triggered rsfMRI network reconfiguration. We then examined the reconfiguration of rsfMRI networks and the stability of their dynamics in response to single-pulse VPM stimulation. Guided by the imaging findings, we conducted simultaneous rsfMRI and electrophysiology experiments to investigate the neural bases underlying the observed rsfMRI network dynamics upon a single neural population activity input from the somatosensory thalamus (Fig. 1D, E).

RsfMRI network dynamics characterized by transient states and transition structures with distinct spatiotemporally evolving patterns

Firstly, we examined the basic properties of rsfMRI network dynamics, such as brain-wide network states and their transitions, in the absence of VPM optogenetic stimulation. Previous studies using Gaussian HMM typically presented region of interest (ROI)-wise rsfMRI network state activation maps based on principal component analysis (PCA) time courses30,45,46. This approach often obscures the spatiotemporal neural activity distributions across multiple voxels within each ROI due to averaging, and it also discards the variance inherent in each imaging voxel due to PCA. Here, we derived rsfMRI network state voxel-wise activation maps from BOLD time courses, enabling the visualization and examination of detailed spatial distributions of neural activities and the polarity of BOLD signals. We further decomposed these states to dissect their temporal properties of the transitions within and between states (Fig. 1C).

The typical homotopic sensory, default mode, salience, limbic, and striatal rsfMRI networks were first computed using seed-based analysis (ROI definition & functional connectivity maps: Supplementary Fig. S1A, B). We then characterized brain-wide rsfMRI network states and their transitions in a data-driven manner via an 18-state HMM model. Eighteen states were chosen based on model consistency47 (i.e., high mean pairwise similarity of state time courses across five iterations), clustering performance48 (i.e., high Calinski–Harabasz score), and model stability49 (i.e., low free energy; Supplementary Fig. S2). States 1–14 exhibited relatively short lifetimes (mean: 1–6 s; standard error of the mean, SEM: ±0.6 s) and high occurrence rates (>17 times/h), indicating their transient and repeatedly occurring properties (Fig. 2A). Notably, these HMM-derived rsfMRI network states resembled combinations of several seed-based rsfMRI networks (Supplementary Fig. S1). States 1, 2, 4, and 5 were dominated by positive (states 1, 2, 4) or negative (state 5) activities in the bilateral somatosensory cortices (primary somatosensory barrel field, S1BF, limb, S1Limb, upper lip area, S1ULp, and secondary somatosensory, S2, cortices), with additional activities in other networks. For instance, state 1 also involved negative activities in auditory (Aud), visual (VC), retrosplenial (RSC) and cingulate (Cg) cortices, while states 2 and 5 showed positive and negative activities in other sensorimotor cortices (Aud, VC, and motor cortex, MC), respectively. Beyond the somatosensory cortices, state 4 also exhibited positive activities in striatal regions (caudate putamen, CPu, and nucleus accumbens, NAc) and default-mode network (DMN)-related regions (Aud, insular area, Ins, Cg, and orbitofrontal cortex, OFC)50. Moreover, state 3 exhibited negative activities in DMN-related regions (RSC, Cg, and OFC), while state 6 exhibited negative striatal (CPu and NAc), hippocampus (HP), Ins, and OFC activities, accompanied by positive activities in the basal forebrain (BF) and hypothalamus (HTh). We observed that states 7–10 exhibited strong brain-wide positive (states 7 and 8) or negative (states 9 and 10) activities located in the DMN-related regions (Aud, Ins, temporal association cortex, TeA, Cg, prelimbic cortex, PrL, and OFC), sensorimotor cortices (primary somatosensory cortices, S1, S2, MC, and VC), and striatal regions (CPu, NAc, globus pallidus and ventral pallidum, GP & VP). Meanwhile, states 11 and 12 exhibited robust negative and positive HP activities, respectively, while states 13 and 14 primarily exhibited positive activities in DMN-related regions (Aud, Ins, TeA, RSC, Cg, PrL, and OFC). It is important to note that we observed activities in other sensory cortices (e.g., Aud and VC), besides the somatosensory cortices, albeit they were often accompanied by other networks, such as HP, RSC, and Cg. Several states that were not detected in most animals and exhibited excessively long lifetimes, high fractional occupancies, and/or low occurrence rates were excluded from our analysis. Specifically, states 15, 17, and 18 were only observed in a few animals and had relatively long lifetimes (>19 s) and high fractional occupancies (>28%), while states 16 and 17 occurred only 3.7±5.2 and 1.4±2.4 times/h, respectively.

We then examined neural information flow across the transient network states (i.e., states 1–14) within the state transition space and identified the clustering of and switching between sensorimotor network (SN)- and DMN-associated states (Fig. 2B, C). In the state transition space capturing the 90th percentile of transitions, states were organized based on the temporal characteristics of their transitions and activation maps. Hierarchical clustering analyses were conducted based on the cross-correlations among state time courses51 and the correlations of state fractional occupancies across animals45, respectively. We showed three distinct hierarchical organizations of rsfMRI network state transitions, wherein states within the same transition structures exhibited temporally correlated dynamics (Supplementary Fig. S3). Sensorimotor network-driven feedforward pathways (SN-FP) were predominantly initiated from state 1, which exhibited strong bilateral somatosensory network activities. We observed a sequential transition route from positive SN activities in state 1 to negative SN activities in state 5 (i.e., state 1 → 2 → 3 → 4 → 6 → 5). Meanwhile, the DMN-dominant recurrent loop (DMN-RL) was marked by a loop transition structure (i.e., state 6 → 7 → 8 → 9 → 10 → 6). DMN-RL featured states with distributed activities in key DMN-related regions (e.g., Aud, Ins, Cg, PrL, and OFC) that fluctuated between positive and negative activities. Notably, states 6 and 8, which exhibited robust activities in BF and HTh, were found to be positioned at the intersection of transitions between rsfMRI networks in DMN-RL and SN-FP. These two transition structures represented the predominant transitions in our rsfMRI network dynamics. We also observed background states that primarily exhibited outputs to other network states but no dominant inputs. These background states were marked by robust activities in HP (states 11 and 12) or DMN-related regions (states 13 and 14). Note that these states and their transitions were derived from the HMM model trained on rsfMRI data acquired both in the absence and presence of optogenetic stimulations. The findings did not differ from the HMM model trained on data without optogenetic stimulation (Supplementary Fig. S4), demonstrating that the predominant transition structures remained similar regardless of the rsfMRI data with or without single-pulse stimulation. For comparison, we also trained 14-state (Supplementary Fig. S5) and 22-state (Supplementary Fig. S6) models. Spatial features and characteristics of the states (e.g., short lifetimes, high occurrence rates), along with the predominant transition structures, remained consistent across models with varying numbers of states. However, the 14-state model tended to underfit the data by merging distinct spatial patterns and transitions (e.g., state 1 → 2 reflecting somatosensory neural activity propagation was not captured), whereas the 22-state model introduced redundant or infrequently occurring states, indicating potential overfitting.

We further characterized the temporal properties of transitions within and between rsfMRI network states (Fig. 2C). To do this, we decomposed each state into three substates: the entering, plateau, and exiting portions of a state’s time course (Fig. 1C). For within-state transitions, we found that the majority of SN-FP states (states 2, 3, and 5), one DMN-RL state (state 10), and most background states (states 11–13) were persistent, exhibiting little to no changes in the spatial distribution of neural activities and the polarity of BOLD signals across three substates. In contrast, the spatiotemporal evolution patterns for the remaining SN-FP states (states 1 and 4), DMN-RL states (states 6–9), and background state (state 14) showed significant changes in neural activity distribution and/or polarity of BOLD signals. For instance, in state 6 of DMN-RL, the negative striatal and cortical activities were strengthened with a larger spatial distribution, whereas the opposite was observed for state 8. Overall, states 1 and 4 of SN-FP and states 7 and 9 of DMN-RL exhibited increased strength and larger spatial distribution from the entering to the plateau substate. In background state 14, activities in DMN-related regions gradually peaked at the plateau but substantially decreased at the exiting substate, with a notable presence of negative striatal and cortical activities. These findings indicated that rsfMRI networks were capable of accommodating both stable and flexible information flow. For between-state transitions, we observed smooth or sharp transitions bridged by substates with or without similar regional activities, respectively. Smooth transitions mainly occurred in state 1 → 2 of SN-FP, and state 6 → 7, state 7 → 8, and state 9 → 10 of DMN-RL, while sharp ones were more common in other SN-FP transitions (e.g., state 2 → 3, state 3 → 4, state 5 → 4, state 13 → 5, state 11 → 2, state 12 → 4). Notably, states 6, 8, and 14 exhibited sharp within-state changes in BOLD activity polarities, which were consistent with the polarities of their dominant input and output states, suggesting their role in facilitating between-state transitions with opposite BOLD activity polarities.

Taken together, our analysis above reveals the transient nature of brain-wide rsfMRI network states and their hierarchical transition structures, characterized by distinct spatiotemporal evolution patterns across substates. Moreover, we underscore the potential role of neuromodulatory regional activities, such as BF and HTh activities in states 6 and 8, in switching between DMN- and SN-associated networks during the resting state.

Single-pulse-evoked VPM activity exhibiting hierarchical somatosensory thalamo-cortical neural activity propagation structures and minimal disruption on rsfMRI BOLD fluctuations

In parallel, we investigated whether sparse single-pulse optogenetic stimulation of the VPM in the somatosensory thalamus (Fig. 1A, C) could be used to interrogate single thalamic activity-triggered changes in rsfMRI networks. Precise targeting of VPM with an optical fiber was confirmed through high-resolution anatomical MRI images during experiment (Fig. 1B) and histological verification of brain slices after experiment (Fig. 1A). Histological slices confirmed that the injection target was in the VPM, where the majority of the expression of optogenetic construct (i.e., ChR2) was found (Supplementary Fig. S7). Expression was also found in parts of the adjacent thalamic nuclei, such as posterior nucleus (PO) and ventral posterolateral nucleus (VPL), due to diffusion of the injected viral vector solution from the middle of VPM42–44,52.

By averaging across numerous (i.e., 3192) optogenetic stimulation events, we demonstrated the expected evoked BOLD responses in somatosensory regions, including ipsilateral VPM, bilateral S1BF, S1Limb, S1ULp, and S2 (ROI definition & activation map: Fig. 3A, B; BOLD single profiles: Supplementary Fig. S8A). Consistent with a previous study53, we also observed BOLD activations in other sensorimotor regions such as visual (lateral geniculate nucleus, LGN, superior colliculus, SC, and visual cortex, V1 and V2), auditory (medial geniculate body, MGB, and Aud), and motor (MC) regions, primarily in the ipsilateral hemisphere. Importantly, we detected activations in higher-order cortices associated with cognition, such as the bilateral Ins and piriform cortex (Pir) and the ipsilateral parietal cortex (PtA), limbic regions, including the bilateral amygdala (Amg), entorhinal cortex (EC), and ipsilateral ventral and dorsal hippocampus (vHP and dHP), and striatal regions, such as ipsilateral CPu, NAc, and GP & VP. Despite ChR2 expression in PO and VPL (Fig. 1A and Supplementary Fig. S7), the single-pulse-evoked BOLD responses also confirmed the specificity of stimulating VPM, as the downstream targets differed from those observed when stimulating PO42 or VPL53. In general, we demonstrate that single-pulse-evoked neural activities propagate brain-wide, beyond the primary sensory regions traditionally detected in our previous fMRI studies52 and without massive averaging, but echo multisite electrophysiological recordings54. This is consistent with our knowledge that a single excitatory spike can trigger multiple postsynaptic spikes, enabling quick amplification of neural activation at downstream targets55.

We found that single-pulse-evoked neural activity propagated along central pathways as demonstrated by natural sensory53 or thalamic pulse-train stimulations42,56, reflecting the hierarchical neural information flow in the somatosensory pathways (Fig. 3C). Initially, neural activity arrived at the ipsilateral somatosensory cortex before engaging the contralateral counterpart. It then propagated from primary somatosensory regions, including S1BF, S1ULp, and S1Limb, to the adjacent higher-order somatosensory cortex, S2. Subsequently, we observed cross-sensory-modality activity propagation from somatosensory to motor, visual, and auditory cortices, followed by cross-system activity propagation from sensorimotor to limbic, striatal, and higher-order cortical regions. To visualize this process, we generated a movie demonstrating the single-pulse-evoked activity propagation (Supplementary Movie S1). Moreover, we found that voxels within the ipsilateral S1 and S2 (the earliest four activated ROIs) were instantaneously and widely activated upon VPM stimulation, whereas voxels in other regions showed more gradual and spatially restricted activation (Fig. 3C). Additionally, previous electrophysiological studies showed that sparse single-pulse stimulation promptly evoked thalamic and cortical activities while avoiding the prolonged and complex neural interactions induced by external sensory or pulse-train stimulations42,52,57. This highlighted the efficiency and simplicity of single-pulse-evoked activity in eliciting hierarchical neural activity propagation.

Variability is an inherent property of spontaneous neural activities and plays a crucial role in shaping evoked activities58 and subsequent behavioral functions59,60. However, under conditions of pulse-train stimulations, both local and distant neural activity can exhibit excessive synchronization, overriding this intrinsic variability59,61. We demonstrated that sparse single-pulse-evoked activity propagation maintained local-to-global neural response amplitude and latency variabilities (amplitude: Fig. 3D; latency: Fig. 3C), indicating minimal disruption to ongoing spontaneous neural activities. Specifically, variations in the local temporal signal-to-noise ratio (tSNR) and BOLD activation amplitudes at the primary projection target of the VPM (ipsilateral/right S1BF, S1BFR) followed a lognormal distribution, signifying local neural response amplitude variability upon single-pulse stimulation. Meanwhile, the strength of long-range activity propagation, measured by cumulative activated voxels excluding S1BFR, reflected global neural response amplitude variabilities. Lognormal distributions, which spontaneous neuronal firing rates typically follow62, imply multiplicative processes involved in both local and long-range propagation of neural activity. We also observed response latency variabilities upon stimulation within each region and across the whole brain, as indicated by the activated voxel distribution and voxel-wise activation latency.

Furthermore, we showed that brain network dynamics, particularly rsfMRI network states 1 and 2 and their transitions, mirrored the hierarchical somatosensory thalamo-cortical activity propagation structures (from S1 to remote brain areas). We found high spatiotemporal similarities between the early S1 activation patterns evoked by stimulation and the activation map of state 1 (highest correlation coefficient = 0.56), and between the activation patterns covering both S1 and subsequent remote regions and the activation map of state 2 (highest correlation coefficient = 0.46; Fig. 3E vs. Fig. 3F). Notably, the sequential transition from state 1 to state 2 (Fig. 3G) coincided with the peak time frames of these spatial similarities (Fig. 3F). Additionally, other higher-order region-associated states, i.e., states 3–14, also exhibited spatial similarities with the later evoked brain-wide activation patterns, albeit to a lesser extent (Supplementary Fig. S8B).

Our findings here demonstrate that single-pulse optogenetic VPM stimulation can serve as an effective tool for probing the response properties of rsfMRI network dynamics while minimally disrupting the intrinsic neural activities and networks. This is supported by the exhibition of hierarchical somatosensory neural activity propagation structures and the preservation of both local and brain-wide neural activity variabilities. Moreover, the flow of neural information in rsfMRI network dynamics (i.e., state transitions) appears to reflect these hierarchical activity propagation structures.

Rapid rsfMRI network reconfiguration while retaining overall stability of network dynamics upon a single thalamic input

To dissect the brain-wide resting functional network dynamics in response to the single thalamic input in VPM, we compared rsfMRI network dynamics in the absence (OG-) vs. presence (OG+) of sparse single-pulse stimulation.

Our results showed that a single somatosensory thalamic input rapidly altered the probability profiles of most transient rsfMRI network states, as evidenced by statistics and large effect sizes (FDR-corrected p value < 0.05 & Cohen’s d value > 0.8; states 1–10, 13, and 14; Fig. 4A, B). For the majority of these states (states 1, 2, 4–7, 9, and 13), the alterations occurred rapidly within ~2 s and returned to baseline levels ~10 s after stimulation. The probabilities were selectively increased (states 1–3 & 9–11) or decreased (states 4–8 & 12–14). Specifically, the probabilities of states with positive somatosensory cortical activities (states 1, 2, 9, 10) increased, while those with negative somatosensory cortical activities (states 5–8) decreased, as expected in the presence of an input from the somatosensory thalamus. Conversely, the probabilities of states with positive activities in DMN-related regions, either with (state 13) or without (state 14) HP activities, and with CPu and S2 activities (state 4) decreased, while the probability of a state with negative activities in DMN-related regions (state 3) increased. Further, the probability of the state with robust positive HP activities (state 12) decreased, with an increase in the state with negative HP activities (state 11). Notably, the latencies of peak probability changes in DMN-RL states followed their corresponding transition structures, with the earliest significant probability change observed in state 6, followed sequentially by states 7, 8, 9, and finally 10 (Fig. 4B). However, the overall state characteristics such as lifetime, occurrence rate, and fractional occupancy remained largely unaltered (OG- vs. OG+: Fig. 2A vs. Fig. 4C and Supplementary Fig. S9. Additionally, we trained HMM models with different dataset combinations to ensure our findings were not impacted by individual animal variability. Here, each dataset excluded one animal from the original, generating 17 and 11 different dataset combinations for OG- and OG+ animal groups, respectively. We confirmed that state time courses generated across different combinations were not skewed by any individual animal (i.e., Gaussian distributed with ~65% of datapoints within one standard deviation of the mean; Supplementary Fig. S10). We also noted that the single thalamic input did not disrupt the within- and between-state spatiotemporal evolution patterns, although, as expected, positive ipsilateral sensorimotor activities increased, accompanied by a parallel decrease in negative bilateral activities (OG-: Fig. 2C; OG+: Fig. 4D).

Further, we demonstrated that a single thalamic input reconfigured the rsfMRI network state transitions while largely maintaining the overall transition pathways and durations. The dominant state transition structures (i.e., pathways shown in the state transition space) were mostly preserved upon stimulation. Only transition probabilities from state 8 to SN-FP initiating state 1 and state 13 → 2 were significantly increased, whereas those of state 11 → 2 and state 2 → 3 were significantly decreased (Fig. 5A, B). Single thalamic input reduced the transition probabilities to DMN-RL and their occurrence rates, while increasing those to SN-FP, without affecting the durations of their respective transition structures (probability: Fig. 5A; occurrence rate & duration: Fig. 5C). We further examined whether such reconfiguration of state transitions was affected by the specific rsfMRI network state that is active at the time of single-pulse stimulation. Overall, we found that transitions toward states 1 and 2 were increased to initiate the SN-FP transition structure, regardless of the ongoing rsfMRI network dynamics at the time of stimulation (Supplementary Fig. S11). Further, when stimulation occurred during an active DMN-RL state, subsequent transitions to maintain the loop structure were suppressed (i.e., state 6 → 7, state 7 → 8, state 8 → 9, state 9 → 10, state 10 → 6; Supplementary Fig. S11B). Notably, states 6 and 8, characterized by robust BF and HTh activities, could transition into both SN-FP and DMN-RL, suggesting their bifurcation role in state transitions between SN-FP and DMN-RL (Fig. 5B). State 6-initiated DMN-RL transition structure exhibited the highest occurrence rate compared to DMN-RLs initiated from other states (states 7–10) at resting state, indicating the crucial role of state 6 in initiating the loop (Fig. 5C). However, in the presence of a brief thalamic input, the occurrence rate of the state 6-initiated DMN-RL structure was significantly suppressed (Fig. 5C). We observed only a moderate decrease (-21%, not significant; Fig. 5D) in the occurrence rate of state 6 → 7 transitions, implicating that other rsfMRI network state(s) might be involved in terminating this loop structure. State 8 showed a significant decrease (-49%) in transitioning to state 9 but a significant increase (+61%) in transitioning to SN-FP through state 1 (Fig. 5D). These findings indicated that a single-pulse somatosensory thalamic input likely disrupted the neural information flow in DMN-RL, diverting it to SN-FP through state 8.

Together, our results reveal that a single somatosensory thalamic input rapidly decreases internally-driven DMN-associated information processing and diverts resources to boost externally-oriented sensorimotor network activities, all while maintaining stable rsfMRI network dynamics such as state characteristics and overall transition structures. This reconfiguration likely involves neuromodulatory hubs, i.e., BF and HTh, which initiate and terminate internally-driven neural information flow and processing upon a thalamic input.

Local wide-band and long-range narrow-band neural activity synchronizations underlying rsfMRI network dynamics and their reconfiguration upon single thalamic input

Having established the key elements of brain-wide rsfMRI network dynamics and rapid network reconfiguration upon a single thalamic input, we further investigated their underlying neural bases. We employed rsfMRI along with simultaneous electrocorticography (ECoG) recordings in key cortical regions of the DMN (i.e., Cg) and SN (i.e., S1; Fig. 1D), providing high temporal and spectral resolution of network dynamics. Note that ECoG is a surface cortical recording and is largely insensitive to deep brain regions. Properties of rsfMRI network states and their transitions were consistent with those acquired without ECoG recordings (Supplementary Fig. S12), and the rsfMRI networks were similarly reconfigured by single-pulse stimulation (Supplementary Fig. S13). We also confirmed the consistency of brain states across animals based on stable physiological ranges (rectal temperature: 36.5–37.5°C, heart rate: 350–420 beats/min, breathing: ~60 breaths/min, oxygen saturation: >95%) and the occurrence rates of neural oscillations (Supplementary Fig. S14).

In the absence of optogenetic stimulations, we found that transient rsfMRI network states were accompanied by distinct patterns of intra- and inter-regional neural activity synchronization, primarily occurring before state onsets, as indicated in state-aligned local Cg and S1 power spectrograms and inter-regional Cg-S1 coherence spectrograms (Fig. 6A). In the local power spectrograms, DMN-RL states and their transition destination states in SN-FP displayed robust wide-band biphasic (states 7–9) or triphasic (states 4–6 & 10) patterns of neural activity synchronizations, whereas other SN-FP states exhibited weaker monophasic patterns (states 1–3; with S1 showing more robustness than Cg; Fig. 6A and Supplementary Fig. S15A, B). Slow (0.2–1 Hz), delta (1–4 Hz), spindle (7–15 Hz), and gamma (>25 Hz) neural oscillatory events in both Cg and S1 coincided with these neural activity patterns, confirming that these local synchronizations were indeed wide-band (OG-: Supplementary Figs. S16, S17). Meanwhile, long-range neural activity synchronizations were primarily pronounced at slow oscillations, as shown in the inter-regional Cg-S1 spectrograms. Most DMN-RL states (states 6–9) exhibited the most robust Cg-S1 coherence, followed by several SN-FP states (states 1, 2, 5), while the remaining states (states 3, 4, 10) exhibited relatively weaker coherence (Fig. 6A and Supplementary Fig. S15C). To explore the relationship between neural activity and rsfMRI network states, we conducted multivariate Granger causality analysis. We observed that local wide-band neural activity synchronizations at resting-state significantly predicted the state time courses (Cg: states 1–4, 6–10 & S1: states 2, 4, 6–8, 10), while long-range low-frequency (<15 Hz) synchronizations predicted the time courses of states 3 and 5 (Supplementary Fig. S18), indicating that neural activity directly contribute to rsfMRI network dynamics.

To investigate the neural underpinnings of the transitions between rsfMRI network states, we examined neural activity synchronizations in ECoG recordings during these transitions. We focused on states that were crucial for mediating transitions between SN-FP and DMN-RL, specifically states 6 and 8. We observed that the local neural activity synchronization patterns in Cg and S1 varied and shared certain similarities with those of their corresponding transition destinations (states 1, 4, 5, 7, 9; Fig. 6A and Supplementary Fig. S19A), supporting the switching between different transition structures. For example, state 6 → 7 and state 8 → 9 transitions within the DMN-RL structure exhibited larger wide-band biphasic or triphasic patterns compared to state 6 → 5 and state 8 → 4/1 transitions out of DMN-RL into the SN-FP structure. In parallel, Cg-S1 coherence was primarily concentrated at slow oscillations, with the exception of the state 8 → 1 (Supplementary Fig. S19B). These pronounced long-range neural activity synchronizations were distributed in distinct time windows mostly prior to the onset of states, depending on whether the states were transitioning to DMN-RL or SN-FP (to DMN-RL: state 6 → 7: -3 to 1 s, state 8 → 9: -4 to -2 s and ~10 s; to SN-FP: state 6 → 5: -4 to 1 s, state 8 → 4: -4 to -1 s).

Subsequently, we examined the underlying neural activity of rsfMRI network states in response to a single brief somatosensory thalamic input. We confirmed that single-pulse stimulation could modulate intra- and inter-regional neural activity synchronizations across all frequency bands (Supplementary Figs. S20, S21). Importantly, the thalamic input consistently modulated local neural activity synchronizations in Cg and S1 across a wide range of frequency bands, as well as their long-range synchronizations at low-frequency bands (especially at slow and delta oscillations). Specifically, in local spectrograms, multiple slow (Cg: states 1–4, 6–10 & S1: states 1–3, 5–10) and delta (Cg & S1: states 1–10) oscillation clusters in most SN-FP states and all DMN-RL states were modulated (OG+ vs. OG-: Fig. 6B and Supplementary Figs. S22, S23). Some states also exhibited modulated theta (4–7 Hz; Cg: states 3, 6, 7, 10 & S1: states 6, 8, 10), spindle (Cg: states 4, 6, 8, 10 & S1: states 3–10), beta (15–25 Hz; Cg: states 5, 7, 8 & S1: states 2, 5–8, 10), and gamma clusters (Cg: states 3, 7 & S1: states 1, 2, 8, 10). Notably, the thalamic input led to substantial attenuation of the triphasic pattern in state 6, making it more similar to that of state 5, while the biphasic pattern of state 8 developed into a triphasic pattern across slow and delta oscillations, resembling state 4. Moreover, in Cg-S1 coherence spectrograms, we observed that slow and delta oscillation clusters of most SN-FP states (slow: states 1, 2, 4, 5; delta: states 2, 4, 5) were altered, while only the coherence of slow oscillations was affected for DMN-RL states (states 6–10; Fig. 6B and Supplementary Fig. S24). Multivariate Granger Causality analysis indicated that both local and long-range low-frequency neural activity synchronizations, which were modulated by single-pulse stimulation, preceded and predicted the reconfiguration of state probability profiles, especially at slow oscillations (Cg: states 1, 7, 8 & S1: states 2, 5–7 & Cg-S1: 1, 2; Supplementary Fig. S25). These results demonstrate that local and long-range low-frequency neural activity synchronizations can directly shape the probability profiles of rsfMRI network states during neural information flow and processing.

Lastly, we examined the effects of a single somatosensory thalamic input on underlying neural activities during rsfMRI network transitions. In local power spectrograms, we identified clusters in the slow and delta frequency bands in transitions involving states 6 and 8 that were altered upon stimulation, along with some clusters in the theta, spindle, beta, and/or gamma frequency bands (theta: state 6 → 7/5, state 8 → 4; spindle: state 6 → 5 and state 8 → 4/1; beta: state 6 → 5 and state 8 → 9; gamma: state 6 → 7 and state 8 → 9/1; Supplementary Fig. S26). These altered neural activity patterns were distinct for transitions within the DMN-RL and SN-FP structures. In particular, state 6 → 7 and state 8 → 9 transitions in DMN-RL featured respective triphasic and biphasic patterns that were largely attenuated across wide frequency bands, while transitions to SN-FP, i.e., state 6 → 5 and state 8 → 4 or 1, only showed modulation at the slow and delta oscillations. Moreover, the brief thalamic input consistently altered long-range Cg-S1 neural activity synchronizations across multiple low-frequency clusters, particularly at slow and delta oscillations (Fig. 6C).

Our findings collectively demonstrate that local neural activity synchronizations across various frequency bands, along with parallel long-range synchronizations at slow oscillations, drive brain-wide rsfMRI network dynamics. Single-pulse somatosensory thalamic input primarily alters DMN- and SN-associated neural activities at low-frequency bands, leading to the reconfiguration of rsfMRI networks to support externally-oriented information processing.

Local gross neural activity synchronizations maintaining the overall stability of rsfMRI network dynamics upon a single thalamic input

We also investigated the neural underpinnings of stable rsfMRI network dynamics in the presence of thalamic inputs, such as unaltered state characteristics (i.e., lifetime, fractional occupancy, and occurrence rate) and overall state transition structures (i.e., spatiotemporal evolutions, dominant transition pathways, and duration of SN-FP and DMN-RL).

We observed largely preserved triphasic, biphasic, or monophasic patterns of gross (i.e., spectrally non-specific) neural activity synchronizations in local Cg (Fig. 7A) and S1 (Supplementary Fig. S27) recordings across transient rsfMRI network states. Quantifications of the ECoG power fluctuations (Supplementary Fig. S28A) at different frequency bands (slow, delta, spindle, and gamma oscillations) and their summations (gross) confirmed comparable fluctuation strengths between experimental conditions (OG- vs. OG+: Supplementary Fig. S28B, S28C). Importantly, we found significant correlations between the peak latencies (both negative and positive peaks) of gross power fluctuations of DMN-RL states and the typical transition sequence of the DMN-RL transition structure (Fig. 7B and Supplementary Fig. S28D). The onset of state 6 was closest to the peaks of power fluctuations, while the onset of state 10 was the furthest, indicating that neural activity synchronizations dictated the predominant transition sequence of DMN-RL (6 → 7 → 8 → 9 → 10). These latencies were also preserved in the presence of a single brief thalamic input, consistent with our observation that the DMN-RL transition structure remained intact. Cg and S1 neural oscillatory events were also unaltered upon single-pulse stimulation, coinciding with the local gross neural activity patterns (Supplementary Figs. S25, S16, S29).

We further examined whether local gross neural activity synchronizations remained unaltered during rsfMRI network transitions. Gross power fluctuations in Cg and S1 were larger in a DMN-RL transition (state 6 → 7) compared to a transition into SN-FP (state 6 → 5), while such fluctuations were comparable for transitions involving state 8 (state 8 → 1/4/9; Fig. 7C and Supplementary Fig. S30). The larger fluctuations observed in state 6 → 7 transitions were also found across multiple frequency bands (Cg: slow, delta, and spindle frequencies, Supplementary Fig. S31A; S1: all frequencies, Supplementary Fig. S31B). Quantifications further confirmed that the presence of a single thalamic input did not alter these power fluctuations. We also noted a pronounced decrease in Cg and S1 low-frequency neural oscillatory events (Cg: delta and spindle; S1: slow and delta; Supplementary Fig. S32) when state 6 transitioned to state 7 to maintain the DMN-RL structure, instead of transitioning into the SN-FP structure (i.e., state 6 → 5), as DMN-RL states were accompanied by substantial changes in neural oscillatory events (Supplementary Fig. S29). Granger causality analysis indicated that local gross ECoG power fluctuations in Cg and S1 largely predicted state time courses (Cg: states 4, 6–8; S1: states 2, 4, 6–10), especially states 6–8 (uncorrected p-value < 0.0001; Supplementary Fig. S33A). These observations were generally preserved upon single-pulse stimulation (Cg: states 6–8; S1: states 2, 6–8; Supplementary Fig. S33B).

These results demonstrate that the overall stability of rsfMRI network dynamics is supported by local gross neural activity synchronizations in both DMN and SN regardless of whether an external neural input is present.

Methods
Animal subjects
All animal experiments were approved by the University of Hong Kong’s Committee on the Use of Live Animals in Teaching and Research (CULATR). Animals were individually housed at 20–26°C, 30–70% humidity and a 12 h light/dark cycle with access to food and water ad libitum. Two groups of adult male outbred Sprague Dawley (SD) IGS rats (originally derived from Charles River Laboratories, USA) were used. Only male rats were used in this study to ensure consistency with previous optogenetic fMRI studies from our group42–44 and others53. Seventeen animals were used in group I for rsfMRI experiments in the absence (OG-) and presence (OG+) of optogenetic stimulation. Of these seventeen, rsfMRI data from six animals were only acquired in the absence of optogenetic stimulation. In group II, six were used for simultaneous rsfMRI and electrophysiological recording experiments, with electrodes placed above Cg and S1. Two animals were identified with poor electrophysiological recording quality in S1 and were removed from subsequent analyses.

Virus packaging
Recombinant AAV vectors were serotyped with AAV5 coated proteins and produced by the Vector Core of University of North Carolina. Viral titer in particles/mL was 2.1×1012 for AAV5-CaMKIIα-ChR2(H134R)-mCherry (Lot# AV4317M). Maps are available online from http://www.stanford.edu/group/dlab/optogenetics.

Stereotactic surgery for viral injection
All stereotactic surgeries for viral injection were performed with rats at 6–7 weeks of age. Rats were anesthetized with an intraperitoneal bolus injection of a ketamine (90 mg/kg) and xylazine (40 mg/kg) mixture. Buprenorphine (0.05 mg/kg, subcutaneous) was administered to minimize pain. Following a small craniotomy, viral injections were performed at two depths in the middle of VPM (-3.6 mm posterior to Bregma, +3.0 mm medial-lateral (ML) right hemisphere, -5.8 and -6.2 mm from the surface of dura). 3 μL of viral constructs were delivered (1.5 μL at each depth) through a 5-μL syringe and 33-gauge beveled needle at 150 nL/min. The injection needle was held in place for 10 min before being slowly retracted from the brain. Scalp incision was sutured, and animals were kept on a heating pad until recovery from anesthesia. Buprenorphine (0.05 mg/kg, subcutaneous) was administered post-injection twice daily for 72 h to minimize discomfort. Enrofloxacin was also administered orally for 72 h to minimize infection and inflammation post-surgery. Animals recovered for 4–5 weeks before conducting MRI and electrophysiological experiments.

Animal preparation for rsfMRI and rsfMRI-ECoG experiments
Surgery was performed under 1.2–2.0% isoflurane to implant an opaque custom-made plastic optical fiber cannula (d=450 μm) at VPM with subsequent high-resolution anatomical MRI verification, following the fiber implantation procedures established in our previous optogenetic fMRI studies5,42–44,52,130,131. The fiber-tip surface was beveled to minimize injury to brain tissue. Animals were only stabilized at 2% isoflurane during the first 15 min of surgery for craniotomy to minimize the duration of exposure to high-dosage anesthesia. Subsequently, animals were maintained under 1.2–1.5% isoflurane for the remaining 20–30 min surgical procedure, i.e., optical fiber implantation. Animals were then maintained under 1.0–1.2% isoflurane during MRI setup and anatomical scans. The fiber was fixed on the skull using dental cement. The fiber cannula was made opaque using heat-shrinkable sleeves and the eyes of the animals were blindfolded in all rsfMRI experiments to prevent light leakage that may cause undesired visual activation during optogenetic stimulation. The tip of opaque optical fiber was typically placed at the center of the VPM nucleus during stereotaxic implantation and verified by high-resolution anatomical MRI. The spatial spread from the fiber tip for 473 nm blue light is rather small (200 μm and 350 μm at 50% and 10% of initial light intensity, respectively) and has little or no spreading in the backward and lateral directions42,132, confining the optogenetic stimulation within the VPM nucleus (Fig. 1B). Heat spread around the tip of the optical fiber was modeled using the well-known Pennes bio-heat equation133, demonstrating that our single-pulse stimulation paradigm minimized the potential thermal effects, with maximal illumination-induced heating remaining below 0.06°C (Supplementary Fig. S34). For simultaneous rsfMRI and ECoG experiments, two insulated platinum-iridium wires were additionally placed above the surface of Cg (+1.5 mm anterior to Bregma, +0.2 mm ML) and S1 (-2.3 mm posterior to Bregma, +5 mm ML) of SD rats for recording electrophysiological signals (127 μm diameter; 90% platinum and 10% iridium composition; impedance ~120 k Ohm at 1 kHz).

After surgery, buprenorphine (0.05 mg/kg, subcutaneous) was administered post-implantation and one drop of 2% lidocaine was applied to the chords to minimize discomfort before endotracheal intubation and MRI experiments. The animals were mechanically ventilated at a rate of 60 breaths per minute with 1.0% isoflurane in room-temperature air using a ventilator (TOPO, Kent Scientific, Torrington, CT). During all fMRI experiments, animals were placed on a plastic holder with heads fixed by a tooth bar and ear bars. The physiology of animals was continuously monitored by sensors of an MRI-compatible system (SA Instruments, Stony Brook, NY) throughout the duration of the rsfMRI experiment. Rectal temperature was maintained at ~37.0°C using a water circulation system. Vital signs were within normal physiological ranges (rectal temperature: 36.5–37.5°C, heart rate: 350–420 beats/min, breathing: ~60 breaths/min, oxygen saturation: >95%) throughout the experiments.

MRI scanner-synchronized optogenetic stimulation and simultaneous rsfMRI-ECoG recordings
A multi-channel neurophysiology recording system (sampling frequency: 24 kHz; notch-filter: 50 Hz, 100 Hz, and 150 Hz; Synapse v86, Tucker Davis Technologies/TDT, Alachua, FL) was utilized to synchronize the MRI scanner trigger, lasers for optogenetic stimulation, and ECoG recordings (Fig. 1D). Computer, light delivery system, and electrophysiological recording amplifier were kept outside the magnet. The MRI-compatible headstage was used for recording electrophysiological signals. Long optical patch cables (5–10 m) delivered light into the bore of the scanner. Blue light was delivered for optogenetic stimulation using a 473 nm Diode-pumped solid-state laser measured before the start of scanning as 8 mW at the fiber-tip (diameter=450 μm, NA=0.5) corresponding to a light intensity of 40 mW/mm2.

To examine the effect of a single thalamic input from VPM on brain-wide rsfMRI network dynamics, rsfMRI scans without and with single-pulse stimulations were interleaved. Each rsfMRI experiment lasted 2~3 h in total. For each optogenetic session, 10 ms short single-pulse stimulation (light intensity=8 mW corresponding to 40 mW/mm2) was presented at VPM42,43,52 once every 15 s (fMRI no. of time points = 600). A total 14 fMRI sessions were performed on each animal.

MRI acquisitions
All MRI experiments were performed on a 7T MRI scanner (PharmaScan 70/16, Bruker Biospin GmbH, Ettlingen, Germany) using a transmit-only birdcage coil in combination with an actively decoupled receive-only surface coil. After placing the animal into the magnet, scout and anatomical RARE T2-weighted (T2W) images were first acquired for accurate positioning and reference with the field of view (FOV)=32×32 mm2, matrix=256×256, RARE factor=8, echo time (TE)=36 ms, and repetition time (TR)=4200 ms. Sixteen contiguous 1 mm slices were positioned in the transverse orientation according to the rat brain atlas to cover the majority of the brain. RsfMRI measurements were made using a single-shot gradient-echo echo-planar imaging (GE-EPI) sequence with FOV=32×32 mm2, matrix=64×64, flip angle=56°, and TE/TR=20/1000 ms. A total of 600 volumes were acquired for each rsfMRI scan.

RsfMRI data preprocessing
RsfMRI preprocessing was performed following the standard procedures established in our previous studies5,43,52 using MATLAB 2022a (MathWorks, Natick, MA) and SPM12 toolbox (Wellcome Department of Imaging Neuroscience, University College London, UK). For each fMRI session, all EPI images were first corrected for slice timing differences and then realigned to the mean images of the first fMRI scan. fMRI scans suffering from motion artifacts (>0.05 mm voxel shifts detected by realignment) and sudden physiological changes (i.e., abrupt changes in respiration pattern, heart rate, and oxygen saturation level) were discarded. After resampling EPI and T2 images to 128×128×16, EPI images of each animal were registered to their corresponding T2 images before they were co-registered to a representative T2W brain template using affine transformation and Gaussian smoothing to maximize normalized mutual information. The transformation matrix was then applied to co-register the rsfMRI images. Voxel-wise linear detrending with least-squares estimation was subsequently performed temporally to remove the baseline drift caused by physiological noises and system instability.

Coherence analysis and time delay quantification of rsfMRI data
For rsfMRI data with single-pulse stimulations, voxel-wise coherence analysis134 was performed on one long trial concatenated from data acquired from all animals (i.e., a total of 3192 optogenetic single-pulse stimulation events) to identify brain-wide activated voxels using the threshold p<0.01; hence, no group-level correction was applied. Coherence value is defined as the magnitude of the frequency of stimulation, f0, (|F(f0)|; in our experiment f0 is 1/15), divided by the sum-of-squares of all frequency components [sum(|F(f)|2)]1/2. BOLD signal profiles were also extracted from atlas-based ROIs for further evaluation and comparison of BOLD activations.

To quantify and compare the time delay of single-pulse-evoked activity propagation at the ROI level, each TR (1 s) of ROI-wise BOLD signal profiles was equally divided into ten time points (0.1 s resolution) to increase temporal resolution. A paired-sample t-test was performed between baseline (3 data points before stimulation) and each time point after stimulus to identify significantly activated brain regions when p<0.01. These brain regions were then sorted according to the order of being activated. To visualize and examine the spatiotemporal activity propagation at the voxel level, voxel-wise BOLD amplitude maps at each time point were similarly generated after equally dividing each TR into ten time points (0.1 s resolution). A paired-sample t-test was performed between the baseline and each subsequent data point after stimulus to determine whether BOLD amplitude was significantly above/below the baseline (0 to 3 s before stimulation). No correction method was applied due to our experimental design, which included a large number of stimulation events (3192) and large variability in BOLD responses (Fig. 3D). Significant data points were identified when p<0.01, while insignificant voxels were masked. Lastly, voxel-wise time-lapsed BOLD amplitude maps were then plotted at 0.1 s temporal resolution to visualize the voxels that were significantly activated/deactivated above/below baseline (Supplementary Movie S1). To further examine the spatial extent of activation within each brain region, the percentage of activated voxels within each atlas-defined ROI at each time point was calculated, where the higher percentage value meant larger activation spreads within brain regions upon single-pulse stimulation.

Seed-based analysis of rsfMRI data
For rsfMRI data without single-pulse stimulation, seed-based analysis (SBA) was employed to investigate rsfMRI functional connectivity. Firstly, a temporal band-pass filter (0.005–0.1 Hz) without spatial smoothing was applied to the pre-processed rsfMRI data, followed by the global signal regression to remove non-neuronal global variance. SBA was then conducted to quantify the interhemispheric rsfMRI connectivity of sensorimotor networks, default mode networks, salience networks, limbic regions, and striatal regions. Two atlas-defined 2×2-voxel regions were selected as the ipsilateral and contralateral seed, respectively. Reference signal profiles were calculated by regionally averaging signals from the voxels within each seed. Pearson’s correlation coefficients (CCs) were calculated between the reference signal profiles and the BOLD signals of every other voxel to generate correlation coefficient maps (CC maps) for each seed location. The group-averaged rsfMRI functional connectivity maps were generated by averaging individual CC maps across animals and applying a threshold (CC>0.14, corresponding to p<0.001) on the group-averaged CC maps.

Hidden Markov Modeling (HMM) of rsfMRI data
To examine the rsfMRI network dynamics and their rapid reconfiguration upon single-pulse thalamic stimulation in a data-driven fashion, we applied the HMM method (https://github.com/OHBA-analysis/HMM-MAR)30 to extract the transient network states and their transitions. The schematic pipeline of dynamic fMRI analysis is shown in Fig. 1C. After standard rsfMRI preprocessing mentioned above, 64 atlas-defined ROI time courses were first extracted and concatenated from pre-processed rsfMRI data without and with optogenetic stimulation (OG-: n=17, 90 trials in total; OG+: n=11, 102 trials in total). PCA was used to reduce dimensionality and increase the signal-to-noise ratio of fMRI data. Approximately 70% variance of ROI time courses was retained by PCA, and HMM was applied to the PCA time courses. Further, the number of states was determined (state number = 18) when the number of states provided the high model’s consistency, clustering performance, and stability (Supplementary Fig. S2). Specifically, the HMM models with different state numbers were repeated 5 times resulting in 5 HMM models per state number. The mean of state time course similarity between pairwise iterations was calculated to measure the model consistency47. The Calinski–Harabasz score was used to examine the performance of HMM states to distinguish the latent clusters of PCA time courses48. It was defined as the ratio between the across-cluster dispersion and within-cluster dispersion. Additionally, the model stability was measured by free energy49. Hence, 18 HMM states were generated composing state ROI-wise mean activation maps, state transition matrix, and state time courses on rsfMRI data without and with optogenetic stimulation. The state transition probability matrices for each experimental condition (OG- and OG+) were separately estimated by applying the model trained on the concatenated dataset, which combined rsfMRI data from both conditions. The difference matrix of state transition probabilities was then obtained by directly subtracting the OG- matrix from the OG+ matrix. Due to the use of a concatenated dataset to estimate state transition probabilities, no statistical correction was performed when comparisons were made between experimental conditions (Fig. 5A). In addition, we also examined the transition probability of the state present at the time of single-pulse stimulation. To do so, we normalized the number of transitions from each state by the total number of stimulation events that occurred during said state, yielding a transition probability matrix (Supplementary Fig. S11A). The fractional occupancy was defined as the proportion of the time spent in each state. In parallel, models with 14- and 22-states were also trained for comparison. To examine whether individual animal variability affects model output, we also trained the 18-state HMM model on different datasets. Here, each dataset excluded one animal from the original groupings, generating 17 and 11 different dataset combinations for OG- and OG+ conditions, respectively. Pairwise similarity of the state time courses was then computed between the newly trained and the original 18-state HMM models.

To dissect the dynamics of rsfMRI networks for higher spatial resolution compared to previous HMM commonly used state ROI-wise activation maps, we generated state voxel-wise activation maps based on the pre-processed rsfMRI data. Single-animal state voxel-wise mean activation maps were obtained by averaging the BOLD amplitudes of each voxel across the time points when the corresponding state was activated. Subsequently, the averaged voxel-wise state mean activation maps were thresholded via group-level masks generated from one-sample group-level t-tests using nonparametric inference with Gaussian random field (GRF, voxel-level p<0.05 and cluster-level p<0.001). Group-level masks were generated using the DPABI toolbox (https://rfmri.org/DPABI)135. Spatial similarity between the state voxel-wise mean activation map and the time-lapsed BOLD amplitude map was quantified using the two-dimensional Pearson correlation coefficient (MATLAB function: corr2), after applying median filtering (2×2 neighborhood) for spatial smoothing. To examine the dynamics at higher temporal resolution, state time courses were temporally decomposed into entering, plateau, and exiting periods based on probability thresholds and/or slopes of state time courses. Specifically, the state entering or exiting periods were defined as the time when the slope of state time courses was positive or negative, respectively. Meanwhile, the probabilities of the state needed to exceed the 95th percentile of the state’s own baseline fluctuations to avoid including the baseline periods associated with low probability values, while remaining below a fixed upper bound of probability=0.9. This approach also ensured that the lower threshold for each state time course can be adaptively defined based on its baseline fluctuations. The state plateau period was defined as the time when the probabilities of the state were above 0.9. Finally, decomposed state probability time courses were aligned to the BOLD time courses to obtain the substate voxel-wise mean activation maps. Temporal organization of rsfMRI network states. To examine the hierarchical organization of state transitions, Ward’s hierarchical clustering analyses136 were conducted based on (i) the cross-correlations between each pair of state time courses (maximum lag=5 TRs)51, with the maximum correlation coefficient used as an index of temporal similarity (Supplementary Fig. S3A), and (ii) the correlations of state fractional occupancies across animals45, quantifying the extent to which pairs of states co-varying in proportion of time spent across animals (Supplementary Fig. S3B).

Electrophysiological signal processing
Electrophysiological signal preprocessing procedures were performed before ECoG analysis (Fig. 1E). MRI gradient artifacts in the raw electrophysiological signals were removed using the fMRI Artifact Slice Template Removal algorithm (FASTR) through FMRIB (https://fsl.fmrib.ox.ac.uk/eeglab/fmribplugin/; the artifact timing was specified using the location of slice-timing events)137. ECoG data of both Cg and S1 channels were then band-pass filtered (0.1–150 Hz) and downsampled to 1024 Hz for subsequent data analyses using MATLAB R2022a. ECoG data quality was verified through visual inspection to exclude data showing periodic or gradient-related noise and to ensure the presence of discernible neural oscillatory events within canonical frequency bands (e.g., slow and gamma oscillations; Supplementary Fig. S35).

Stimulus-evoked neural activities were removed based on Averaged Artifact Subtraction (AAS)137 to examine spontaneous neural activity influenced by single-pulse stimulation (Supplementary Fig. S21). Specifically, the template of evoked activity was obtained by averaging the ECoG data across 100 ms windows after the onset of stimulation. The rescaled template was then computed by multiplying the template with the regression coefficients obtained from regressing the signal within the window against the template. The ECoG data for post-processing analysis was finally acquired by subtracting the rescaled template signal from the pre-processed ECoG data.

To measure the local Cg and S1 and their inter-regional neural activity synchronizations across both temporal and frequency domains, we calculated the power spectrograms of Cg and S1, and Cg-S1 coherence spectrograms138 using EEGLAB toolbox (Swartz Center for Computational Neuroscience, University of California San Diego, USA). Afterward, the power and coherence spectrograms were aligned to the onset of each state. To identify the significant state-aligned intra- and inter-regional neural activity synchronizations, a post hoc permutation test based on the power or coherence spectrograms against a baseline period (12.5 s to 11.5 s before and 4.5 s to 12.5 s after state onsets) with FDR correction (p-value<0.05; shuffled 10,000 times) was performed to generate the mask of significance, and small clusters (<100 pixels) were removed. To examine the influence of single-pulse stimulation on the state-aligned intra- and inter-regional neural activity synchronizations, we performed a post hoc permutation test based on the differential state-aligned power or coherence spectrograms (OG+ - OG-; p-value<0.05 or <0.01; shuffled 10,000 times) to generate the mask of significance and removed small clusters (<100 pixels). To further confirm these consistently altered time-frequency clusters upon stimulation, we then performed a comparison of cluster-wise averages using paired t-tests between rsfMRI data without and with single-pulse stimulation conditions.

To assess how the local neural activity synchronization contributes to the overall stability of rsfMRI network dynamics upon stimulation, ECoG power calculation and neural oscillatory event detection were performed. For ECoG power calculation, the ECoG data was firstly band-pass filtered to preserve the targeted frequency range, such as slow oscillatory (0.2–1 Hz), delta (1–4 Hz), spindle (7–15 Hz), and slow gamma (25–45 Hz) frequency bands. To extract the gross frequency band, the pre-processed ECoG data were filtered using a zero-phase Butterworth band-pass filter with cutoff frequencies of 0.2 and 45 Hz. Then, we computed the root mean square (RMS) values of filtered ECoG data with a window size of 125 ms and a step of 62.5 ms. To highlight the ECoG power co-fluctuated with the onsets of states, the RMS values were normalized to mean of time windows devoid of period when ECoG power consistently fluctuated, i.e., 9 s to 12 s before the state onsets. We further quantified the ECoG power fluctuations by calculating the area under the curve (AUC) for both gross and band-specific RMS values (Supplementary Fig. S28A). For neural event detection, slow, delta, spindle, and gamma neural oscillatory events were detected based on established methods43,139. Specifically, the ECoG data was firstly band-pass filtered at the corresponding frequency bands. Subsequently, to detect slow oscillatory and delta events, all zero crossings of bandpass-filtered ECoG data were calculated. Event duration was determined by the time between two successive positive-to-negative zero crossings. Slow oscillations were determined when the duration was between 0.8 s to 3 s and the amplitude exceeded the mean (μ) plus double standard deviation (σ), while delta events were determined when the duration was between 0.1 s to 1 s and amplitude exceeded μ+2×σ. To detect spindle and gamma events, the envelope was first calculated via Hilbert transform. Spindle events were detected when the envelope exceeded the duration threshold defined as μ+1.5×σ for at least 500 ms but not longer than 3 s and the amplitude threshold defined as μ+2.5×σ for at least one data point. Further, the number of cycles within single spindle events needed to be ranged from 4 to 25. Gamma events were identified when the envelope exceeded the duration threshold defined as μ+3.5×σ for at least 30 ms but not longer than 300 ms. Representative detected slow, delta, spindle, and gamma events were shown in Supplementary Fig. S35. Gross ECoG power fluctuations and neural event occurrence rates were then aligned to HMM state onsets, where the probabilities of the state were larger than 2/3. Comparisons between experimental conditions of each transient state were done by performing paired t-test with FDR correction, while comparisons across state transitions (transition into SN-FP or DMN-RL) and experimental conditions of bifurcating states (states 6 and 8) were performed through two-way ANOVA with Turkey correction.

Multivariate Granger causality analysis
Here, analysis140,141 was conducted on ECoG signals such as band-limited ECoG power (in local Cg and S1) and inter-trial coherence (inter-regional between Cg-S1), and the state time courses. To examine how neural activity synchronizations at different frequency bands influence rsfMRI network dynamics, all data were aligned to the onset of each state. 24-s windows centered on each state onset (12 s before and after) were extracted and averaged within each animal. To examine how these synchronizations contribute to neural activity-triggered reconfiguration, data were aligned to single-pulse stimulation onset. Here, 14-s windows centered on stimulation onset were extracted and averaged within each animal. Following the respective data preprocessing steps, the resulting ECoG signals were subsequently smoothed (moving average with window size=5 data points) and downsampled to match the temporal resolution of the fMRI signals (i.e., one sample per second). We then fitted multivariate vector autoregression (VAR) models to the band-limited ECoG power/coherence signals and state time courses. Granger causality analysis was then applied to assess whether local and long-range neural activity synchronizations across different frequency bands could predict the temporal dynamics of rsfMRI networks and network reconfiguration.

Granger causality analysis
The gross ECoG power fluctuations were similarly smoothed and downsampled to the same sampling rate as the fMRI signals. The VAR model was applied between ECoG power fluctuations and state time courses for each animal. To capture the most relevant temporal dependencies, we fitted VAR models with lags ranging from 1 to 5 s. The optimal lag was determined based on the lowest akaike information criterion (AIC) value142. We then performed the Granger causality test to assess whether neural power fluctuations significantly predicted changes in state time courses. The z-scores obtained from Granger causality analysis were averaged across animals to obtain group-level z-scores, which were subsequently converted to group-level p-values. Since we aimed to examine how thalamo-cortical oscillations contribute to rsfMRI network dynamics, we focused on whether ECoG power fluctuations predicted state time courses, rather than the reverse.

Histology, immunohistochemistry, and confocal imaging
To confirm the specific expression of ChR2-mCherry in the VPM excitatory neurons, histology was performed according to the previously published procedure42,43. Upon completion of all rsfMRI and rsfMRI-ECoG experiments, animals were anesthetized with pentobarbital and then transcardially perfused with ice-cold 4% paraformaldehyde (PFA) in phosphate buffered saline (PBS). The brains were then equilibrated in 20% sucrose in PBS at 4°C overnight. Axial sections (40 μm) were prepared on a freezing microtome (model 860, AO Scientific Instruments). Consecutive sections were mounted and examined with a laser confocal microscope (Carl Zeiss LSM780). For immunohistochemistry, free-floating sections were processed with 5% normal goat serum and 0.3% Triton X-100 in PBS with primary antibodies against rabbit polyclonal to CaMKIIα (1:400; Abcam) at 4°C for 24 h. After washing with PBS, sections were then incubated for 2 h at room temperature with secondary antibodies Alexa Fluor 647 conjugate goat anti-rabbit IgG and Alexa Fluor 488 conjugate goat anti-guinea pig IgG (both 1:500; Molecular Probe). Slices were then washed and mounted using FluoroShield mounting medium with DAPI (Abcam). Double or triple immunofluorescence was assessed with a laser confocal microscope (Carl Zeiss LSM780/LSM980).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript by Xie et al. presents an impressive set of experiments integrating optogenetic stimulation, resting-state\nfMRI, and electrophysiology to examine how brain-wide network dynamics respond to a short stimulation of thalamic\nneurons. The authors use Hidden Markov Modeling (HMM) to dissect transient network states and their transitions, and\nprovide evidence that the somatosensory thalamus can drive reconfiguration of large-scale networks. Overall, the data are of\nhigh quality and the methodological rigor is commendable. The authors’ approach is powerful and the results are potentially\nimpactful. Nevertheless, there are a few conceptual and methodological points that would benefit from\nfurther clarification.\n1. To help readers appreciate the contribution of this work, it may be helpful for the authors to more clearly articulate the key\nneuroscientific and technical questions addressed, and to clarify how the study builds upon or extends existing findings -\nparticularly in light of prior literature showing DMN suppression and sensorimotor activation following stimulation.\n2. The manuscript introduces an important question of whether rsfMRI networks can facilitate rapid neural information flow\n(line 74), which represents a significant knowledge gap in my opinion. The analyses seen subsequently in the manuscript,\nhowever, focused on showing the neuronal events preceding fMRI-derived brain states. Since the authors have access to\nthis beautiful dataset with concurrent electrophysiological and fMRI data, I wonder whether it is worthwhile to tap into this\ntopic to better echo one of the goals set by the authors originally.\n3. The term “rapid reconfiguration” was used throughout the manuscript, although the timing reported appears to initiate from\n2 seconds and normalize at 10 seconds. While this timing is in agreement with a related study (PMID: 36791185, where\nthose authors failed to provide a definitive explanation for the delay), it still falls within the hemodynamic delay window.\nThus, I wonder if the use of the term “rapid” could be further clarified, also in light of the 1 second fMRI temporal resolution.\n4. Studies dissecting brain states often place less emphasis on the spatial features of specific brain states. The authors\nshould be commended for showing detailed patterns representing brain states, as this is a step closer to linking findings with\nestablished brain networks. I noticed that the authors chose a relatively large number - 18 brain states. The authors may\nwant to provide further justification for selecting this number and evaluate whether the main features and state transition\npatterns persist under a different number of states. Understanding the relevance of these findings in the context of potential\noverfitting or underfitting would strengthen the study’s rigor.\n5. On the technical aspect, the choice of 18 brain states using HMM may also raise concern regarding redundancy and\ninterpretability. With the use of initial 18 states, it seems worthwhile to confirm that these are functionally distinct by\ncomputing similarity (e.g., Pearson correlation) between state activation maps and cluster similar states into meta-states.\nThis would help identify whether there is functional overlap among states and whether timepoints could be reassigned more\nparsimoniously. Alternatively, describing all 18 states in detail may dilute the clarity of the manuscript. One potential strategy\nis to focus on the dominant or most functionally relevant states during stimulation and non-stimulation conditions, and\naddress which states exhibit significantly higher lifetime, fractional occupancy, or occurrence rates under each condition.\n6. While it is helpful to display spatial patterns graphically in the main figures, I believe it would also be useful to show state\ntransition probabilities in a matrix format in the main figures, as done in some supplemental figures (e.g., S3).\n7. The definition of “DMN-related” states could benefit from more clarification. It appears that several background states also\nstrongly resemble DMN features.\n8. In Line 215, states 6 and 8 are identified as occurring between two major networks and their transition probabilities can be\nextracted. However, I am not sure this necessarily means that these brain regions “facilitate” switching. I would like to hear\nmore about the authors’ view on this point.\n9. Since the authors concatenated all data for analysis, do we know whether certain individual subjects contributed more\ndominantly to the results? That is, how does individual variability impact the reported findings?\n10. Brain state transitions are ever-changing (though following probabilistic patterns as shown). Does the brain state at the\ntime of thalamic stimulation matter for the outcome? Could the timing of stimulation with respect to current state dynamics\ninfluence the resulting reconfiguration?\n11. The term “hierarchical” is used to describe transitions, yet no formal analysis of hierarchy is presented. It would be better\nto avoid this unless further analysis is provided.\n12. The manuscript describes “efficient” neural processing. It is recommended to provide metric or model to support the use\nof this term.\n13. Regarding statistical analyses, Figure 5A (Difference Matrix) and 5C appear to involve multiple comparisons. Were\nthese corrected (e.g., FDR)? Please specify the correction method used.\n14. In Figure 3F–G, the label “OG− (n = 17)” suggests a no-stimulation condition, but the figure content appears inconsistent\nwith that label. Please clarify whether this panel reflects stimulation or baseline.\n15. The manuscript does not provide detailed methodology for how transition matrices for OG− and OG+ were computed.\nWere they derived from separate HMM models or the same model? Were transitions computed per animal and averaged, or\nacross the full concatenated data? Adding this information would improve reproducibility.\n16. The authors appear to imply causal inference between electrophysiological and fMRI findings. While the time-locked\ndisplay of these results is compelling, I feel additional analysis may be warranted to support the claim that low-frequency\nthalamic oscillations and inter-regional coherence “contribute” to fMRI connectivity and reconfiguration - rather than these\nbeing merely parallel observations.\n17. When comparing EEG power and state probabilities, the small delay between the peak electrophysiological power\nchanges and state timing appears to vary by state. Is there any potential explanation for this? For example, could it be due to\ndifferences in the hemodynamic response function across brain regions?", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_3": "The study combines optogenetics, fMRI, and electrophysiology to investigate resting-state network dynamics, which is\ntechnically impressive. The discovery that single thalamic inputs can rapidly reconfigure resting-state networks (~2s) adds\nnew insight on thalamic inputs modulating the brain intrinsic networks. The use of Hidden Markov Modeling (HMM) to\ncharacterize transient network states and transitions is methodologically robust. Potentially the findings on how thalamic\ninputs modulate DMN and sensorimotor networks could inform research on disorders like schizophrenia or Alzheimer’s.\nMajor concerns and suggestions:\nThe study uses light isoflurane anesthesia, which may alter brain intrinsic networks and neuronal activity compared to\nawake states. I would suggest discussing how anesthesia might affect the observed network reconfiguration (e.g.,\nsuppressed arousal mechanisms) and compare with prior awake rodent/human studies.\nThe optogenetic stimulation is elegant but raises questions on stimulation range, please provide histological verification of\nstimulation sites post-experiment. If any Off-Target virus injection that could be valuable control, or include a control group\nwith non-ChR2 animals or discuss limitations. A good control design will be very important for unexpected effects, such as\nsurgery/fiber damage effects/laser thermal stimulation and so on since fMRI signals may be so sensitive to be changeable to\nall introduced facts.\nThe study reports wide-band (local) and narrow-band (inter-regional) synchronizations but could better explain: How do slow\noscillations (<1Hz) specifically drive BOLD changes? Does optogenetic stimulation directly alter slow oscillations, or are\nchanges secondary to other frequencies? Add cross-frequency coupling analysis or Granger causality to clarify directed\ninfluences.\nSome key analyses lack detailed statistics (e.g., effect sizes for state probability changes, cluster correction methods for\nfMRI). Report exact p-values (not just thresholds) for major findings. Clarify multiple comparison corrections (e.g., FDR vs.\ncluster-based). Include confidence intervals for temporal metrics (e.g., 2s reconfiguration latency).\nThe DMN suppression and SN enhancement are intriguing but could be contextualized better: How might this rapid\nreconfiguration translate to perception/attention in awake animals? Contrast with task-evoked fMRI studies showing similar\nDMN/SN trade-offs. Link findings to theories of \"network switching\" (e.g., Menon’s triple-network model).\nHMM Parameters: Why were 18 states chosen? Provide model selection criteria (e.g., free energy plots). ECoG\nLimitations:Surface recordings may miss deep structures (e.g., thalamus). Acknowledge this. Also, provide the ECoG\nmethod details for a quality verification and reproducibility. Share code/HMM models for reproducibility. The concept of\ndecomposing states into entering, plateau, and exiting substates is interesting for analyzing within-state dynamics. The\nMethods section describes the HMM process well but could slightly elaborate on precisely how the boundaries for these\nsubstates were defined based on the state probability time courses (e.g., using thresholds on the probability or its\nderivative?).\nThe Discussion compellingly links the observed network reconfiguration to altered neural synchrony, particularly in low\nfrequencies. Could the Discussion speculate slightly further on the potential causal cascade? Does the thalamic input\nprimarily affect sensory cortex synchrony, which then propagates, or does it potentially influence key modulatory nodes (like\nBF/HTh implicated in states 6/8 ) more directly to orchestrate the shift? The text currently implies synchrony changes lead\nthe reconfiguration (lines 426-427 ), but acknowledging the potential role of hub regions like BF/HTh in mediating this could\nadd nuance.\nMinor Suggestions:\nAdd schematic of hypothesized circuit mechanism (e.g., thalamus → BF → DMN suppression). Improve color contrast in\nspectrograms (Fig. 6) for readability. Define \"internally/externally oriented processing\" early in the Introduction. Avoid\noverstating \"first to show\" claims—cite prior work on fast fMRI dynamics (e.g., Yu et al., 2016).\nThis is a technically outstanding study with novel insights into rapid network dynamics. With revisions addressing the above\npoints, it would make a strong contribution to Nature Neuroscience. I recommend major revision.", "reviewer_4": "In this study, authors investigated how brain-wide rsfMRI networks enable rapid and efficient neural information flow and\nprocessing. By using a single-pulse optogenetic stimulation of the somatosensory thalamus, authors showed activity-\ntriggered rapid reconfiguration of rsfMRI networks on the order of a few seconds (~ 2 s to begin reconfiguring; ~ 10 s to\nnormalize). Combining simultaneous electrophysiology and fMRI measurements, they further revealed that this\nreconfiguration was driven by distinct synchronizations of neural activity within and between networks occurring across\nspecific time and frequency windows. This is a very important study with a tremendous amount of data. The analysis is\ncomplex but rigorous. I only have a few minor comments.\n1. The results of the study are largely descriptive, and it might be difficult for readers to clearly grasp the picture of brain state\ntransitions reported herein. I wonder, in addition to laying out all observations, whether authors can consider a way to\nsummarize their findings in a way that might be easier for readers from outside of the field to understand?\n2. How stable are the results at a different probability threshold?\n3. The finding that a single somatosensory thalamic input at VPM largely maintains the overall stability of brain dynamics is\nconsistent with the findings in a recent study (PMID: 35179203), despite the differences in the model of stimulations and\nanimal’s physiological condition (awake vs anesthetized). Authors may consider discussing these findings in a larger\nframework."}, "author_response": "RESPONSES AND REVISION OF NATURE COMMUNICATIONS\nMANUSCRIPT # NCOMMS-25-18143\nBrain-wide Resting-state fMRI Network Dynamics Elicited by Activation of Single Thalamic Input\nDear Reviewers and Editor,\nWe appreciate that all reviewers recognized the novelty and significance of our work, i.e., the\ncombination of single-pulse optogenetics, simultaneous resting-state functional MRI (rsfMRI) with\nelectrophysiology recordings, and dynamic network analysis, to reveal previously unidentified (i)\nbrain-wide rsfMRI network dynamics capable of rapidly and sensitively support the processing of\nsingle somatosensory thalamic input, and (ii) their underlying neural mechanisms.\nIn this revision, we conducted additional analyses and revised the manuscript substantially. In\nparticular, we addressed concerns regarding the contributions of neural activity synchronizations\non rsfMRI network dynamics and their reconfiguration, as well as the robustness of our model\noutputs across different state numbers and training datasets. Additionally, we also included the\ndata and codes used in this study in a public repository.\nWe would like to take this opportunity to thank the four reviewers and the editor for their valuable\ncritiques and suggestions, which we believe have substantially strengthened this manuscript.\nBest Regards,\nAll Authors\nReviewer #1\nR1-G This manuscript by Xie et al. presents an impressive set of experiments integrating optogenetic\nstimulation, resting-state fMRI, and electrophysiology to examine how brain-wide network dynamics\nrespond to a short stimulation of thalamic neurons. The authors use Hidden Markov Modeling (HMM) to\ndissect transient network states and their transitions, and provide evidence that the somatosensory\nthalamus can drive reconfiguration of large-scale networks. Overall, the data are of high quality and the\nmethodological rigor is commendable. The authors’ approach is powerful and the results are potentially\nimpactful. Nevertheless, there are a few conceptual and methodological points that would benefit from\nfurther clarification.\nWe thank the reviewer for his/her encouragement and views on the importance of this study.\nclearly articulate the key neuroscientific and technical questions addressed, and to clarify how the study\nbuilds upon or extends existing findings - particularly in light of prior literature showing DMN\nsuppression and sensorimotor activation following stimulation.\na. Thank you for the suggestions to better signify the study.\nRsfMRI network dynamics are constantly modulated by a wide array of spontaneous neural\nactivities arising from diverse sources1-3, making it challenging to investigate the mechanisms\nby which these dynamics can rapidly process specific neural inputs. We address this by\nemploying a single-pulse optogenetic stimulation paradigm in ventral posteromedial thalamic\nnucleus (VPM), which isolates a well-defined hierarchical somatosensory thalamocortical neural\nactivity propagation pathway and avoids the complex, long-lasting neural cascades induced by\nnatural sensory stimuli.\nFor instance, a simplistic sensory stimulus used in previous studies of rsfMRI network dynamics,\ne.g., an oriented bar in the visual field4, whisker deflection5, an odor6, or an auditory tone7, can\ntrigger the sequential occurrence of synaptic excitation and inhibition in their corresponding\nsensory thalamic nuclei and cortices. Further, sensory stimuli also induce prolonged and complex\ninformation flows from multiple parallel pathways. For example, the study referred to by the\nreviewer in R1-3 below reported that a 50-ms acoustic oddball stimulus induced sustained\nGCaMP signal suppression in DMN-related regions and enhancement in salience networks,\nlasting from 0.5 to 8 s8. Similarly, in our prior work, a 100-ms flashing light evoked neural\nresponses in the visual cortex that persisted for 500 ms9 after cessation of stimulus.\nIn contrast, single 10-ms optogenetic pulses in VPM, as demonstrated in both our current (Fig.\nS20) and previous work9, only evoked brief neural activity lasting < 50 ms in S1, avoiding the\ncomplex cascades caused by prolonged or repetitive stimulation. Such a short evoked neural\nactivity is akin to a delta-function input (i.e., an idealized signal useful for analyzing dynamic\nbehaviors of any complex system) that offers a spatiotemporally controlled framework for probing\nthe rapid reconfiguration of large-scale rsfMRI networks.\nb. We would also like to clarify how our findings extend prior works reporting default mode network\n(DMN) suppression and sensorimotor network (SN) enhancement in response to external tasks or\nsustained stimulation10-12. While these studies demonstrate similar DMN suppression and SN\nactivation, they typically involve complex, prolonged neural activity across multiple pathways,\nmaking it difficult to isolate whether and how rsfMRI networks can reconfigure rapidly in\nresponse to single and brief neural events to support fast information processing. Note that rsfMRI\n2\nnetwork dynamics were characterized by two major transition structures, DMN-dominant\nrecurrent loop (DMN-RL) and SN-driven feedforward pathways (SN-FP; Fig. 2C). Our study\nrevealed that rsfMRI networks can be reconfigured rapidly within 2 s by a single thalamic input.\nSpecifically, we found that when any of the DMN-RL states were active at the time of single\noptogenetically evoked somatosensory thalamic input (see also detailed response to R1-10 and\ncorresponding Fig. R1-10 below), transitions to maintain the loop structure were suppressed to\nenhance those to state 1 or 2 of SN-FP. Further, state 8, which is located at the intersection\nbetween DMN-RL and SN-FP, mediated the rsfMRI network reconfiguration (Fig. 5D and Fig.\nR1-10B). Such network reconfiguration was driven by low-frequency neural activity\nsynchronizations in primary somatosensory (S1) and cingulate (Cg) cortices (see detailed analyses\nto examine the relationship between neural activity and rsfMRI dynamics in R1-16 below).\nNotably, state 8 exhibited activities in the insula (Ins) and orbitofrontal cortex (OFC), key regions\nimplicated in the rat salience network13,14, and in the basal forebrain (BF; Figs. 2C and 4D, and\nFig. S1). Salience network has been closely linked to bottom-up attention15, which is rapidly\ndeployed in response to sensory stimuli16-18, while BF activity can bias behavioral states toward\nexternally-oriented (e.g., bottom-up attention) or internally-oriented (e.g., grooming) processes19-\n21. Based on the theoretical framework proposed by the Menon triple-network model15, we\npropose that the neural activities within such salience network-related regions with prominent\ninputs from the BF mediated the switching between task-positive (i.e., SN) and task-negative (i.e.,\nDMN) networks, thereby prioritizing the rapid processing of a single bottom-up somatosensory\nthalamic input.\nTogether, our study uncovers the rsfMRI network dynamics of DMN suppression and diversion to\nSN, and the neural activities that drove such a rapid reconfiguration.\nc. In this revision, as per the discussions in R1-1a above, we added statements in the Introduction\nsection to (i) better articulate the key neuroscientific questions addressed and (ii) accentuate the\nadvantages of deploying single-pulse optogenetics to dissect brain network dynamics.\nIn this revision, we also incorporated new insights provided in our study in the Discussion section\nin light of prior experimental and theoretical works on DMN suppression and sensorimotor\nactivation following stimulation as discussed in R1-1b above.\nIn this revision, to help readers appreciate the contribution of the work (as pointed out also by\nReviewer 4 in R4-1 below), we added a new Fig. 8 (see also Fig. R4-1 below) to illustrate the key\nfindings of the study. The figure was cited in the Discussion section paragraph summarizing the\nwork.\nneural information flow (line 74), which represents a significant knowledge gap in my opinion. The\nanalyses seen subsequently in the manuscript, however, focused on showing the neuronal events\npreceding fMRI-derived brain states. Since the authors have access to this beautiful dataset with\nconcurrent electrophysiological and fMRI data, I wonder whether it is worthwhile to tap into this topic to\nbetter echo one of the goals set by the authors originally.\na. As also addressed in our response to your queries in R1-16 below, we performed Granger\ncausality analyses to show that local and long-range low-frequency (< 15 Hz) cortical neural\nactivity synchronizations preceded and significantly predicted rsfMRI network reconfiguration\nupon single-pulse stimulation. These findings indicate that rsfMRI network dynamics are closely\ndriven by underlying neural activity synchronizations.\nHowever, we would like to clarify that our use of the term “facilitate” was not intended to indicate\nthat rsfMRI networks were the origin or causal driver of neural information flow. Rather, we\npropose that the dynamic structure and rapid reconfigurability of rsfMRI networks enable the\npropagation of neural information flow through specific anatomical and functional pathways,\nthereby supporting fast processing of a single thalamic input.\nPrior human and animal studies showed that propagation of spontaneous brain activity can reflect\ncharacteristics of task-evoked activity22-25. For example, spontaneous memory reactivations have\nbeen observed to propagate across distributed networks that resemble those recruited during task\nperformance23,24. In our data, we observed that, in the absence of stimulation, rsfMRI state\ndynamics (e.g., state 1  2) mirrored the hierarchical information flow in somatosensory\npathways26-28 (Fig. 3F and 3G), suggesting that rsfMRI networks inherently reflect structured\nroutes for task- or stimulus-relevant information flow. Upon stimulation, we observed an increase\nof occurrence rates of SN-FP structure and a suppression of DMN-RL (Fig. 5C), indicating that\nactivity-triggered reconfiguration of rsfMRI networks rapidly and selectively engaged neural\npathways for externally-oriented information processing while suppressing those for internally-\noriented information processing.\nb. Given that our results did not demonstrate that rsfMRI dynamics causally influenced neural\ninformation flow, the use of the term “facilitate” could be misinterpreted as implying that rsfMRI\nnetwork dynamics were the source or driver of information flow.\nIn this revision, we revised the Introduction, Results and Discussion sections to more precisely\ndescribe how rsfMRI networks support/accommodate fast information flow and processing.\nappears to initiate from 2 seconds and normalize at 10 seconds. While this timing is in agreement with a\nrelated study (PMID: 36791185, where those authors failed to provide a definitive explanation for the\ndelay), it still falls within the hemodynamic delay window. Thus, I wonder if the use of the term “rapid”\ncould be further clarified, also in light of the 1 second fMRI temporal resolution.\na. Thank you for the comment. We agree that the use of the term \"rapid reconfiguration\" of rsfMRI\nnetworks requires further clarification, considering the observed timing (~ 2 seconds to initiate\nand ~ 10 seconds to normalize) and the limitations associated with hemodynamic delays and the\n1-second temporal resolution of fMRI data acquisition.\nWhile our observed timing of rsfMRI network reconfiguration is within the hemodynamic delay\nwindow, we posit that it could be faster as previously noted in the Discussion section. This is\nbecause hemodynamic responses in the presence of an appreciable neural population activity can\nbe captured within 400 ms when measured using a line-scan fMRI sequence with millisecond and\nmicrometer resolutions (i.e., TR of 100 ms and in-plane 150 μm)29 and optical imaging30. Hence,\nthe delay in hemodynamic response to neural activity is the predominant limiting factor in\ndetermining the temporal resolution of fMRI measurements. As such, we contend that the response\nof rsfMRI network states to a single thalamic input were likely shorter than what we observed with\na GE-EPI sequence of TR = 1s and a spatial resolution of 0.5 x 0.5 x 1.0 mm3.\n\"Rapid\" reconfiguration refers to the prompt neural activity-triggered reconfiguration (~ 2 s to\nbegin reconfiguring; ~ 10 s to normalize; Fig. 4A) of multiple large-scale rsfMRI networks that\nwe were able to discern with a single thalamic input. This contrasts with previous studies of\nrsfMRI network dynamics (i.e., also discussed in R1-1a above), which did not quantitatively\ncharacterize the timescales of such dynamics31-33. These studies were mostly qualitative whereby\nthe timescales for rsfMRI reconfiguration were observed at > 10 s and/or persisted for hours,\nlikely due to relatively longer neural perturbations (e.g., > 10 s or ~ 2 mins). Note that the ~6 s\ndelay observed between peak BOLD fMRI activity in DMN and peak fiber photometry GCaMP\nsignals in DMN regions in the absence of any stimulus (i.e., study referenced by reviewer8), while\nin agreement with the delays between rsfMRI states and neural activity synchronizations (Fig.\n6A), do not pertain to the reconfiguration of rsfMRI networks. The study however showed an ~ 2\ns onset delay for peak/trough GCaMP signals/neural activities in certain cortical DMN regions\nfollowing a 50-ms acoustic oddball stimulus. In fact, we documented instantaneous changes to\nneural activity synchronizations at S1 and Cg (i.e., peaked within 200ms; Fig. S19) that preceded\nthe rsfMRI network reconfiguration onset delay of 2 s (see detailed discussions in R1-16c and\ncorresponding Table R1-16-3 below), demonstrating how rapid our characterizations were.\n5\nb. In this revision, as per discussions in R1-3a above, we clarified in the Introduction section that it\nis presently challenging to demonstrate that rsfMRI network can be rapidly reconfigured to\nprocess specific neural input. This is to better distinguish our work from previous attempts and to\ncontextualize the rapid timescales of rsfMRI reconfiguration characterized following activation of\na single optogenetically evoked somatosensory thalamic input.\nIn this revision, we exclusively used “rapid” in all sections of the manuscript to refer to the\ntimescale documented for rsfMRI network reconfiguration.\nIn this revision, we now used the term \"fast/prompt\" throughout the manuscript to describe neural\ninformation processing occurring on the millisecond or second time scale instead of “rapid”.\nR1-4 Studies dissecting brain states often place less emphasis on the spatial features of specific brain\nstates. The authors should be commended for showing detailed patterns representing brain states, as this\nis a step closer to linking findings with established brain networks. I noticed that the authors chose a\nrelatively large number - 18 brain states. The authors may want to provide further justification for\nselecting this number and evaluate whether the main features and state transition patterns persist under a\ndifferent number of states. Understanding the relevance of these findings in the context of potential\noverfitting or underfitting would strengthen the study’s rigor.\na. We thank the reviewer for highlighting the importance of spatial features in brain state\ncharacterization and for recognizing our efforts in detailing the spatial activation patterns of\nrsfMRI network states.\nb. Regarding the choice of using 18 states in our Gaussian Hidden Markov model (HMM), we agree\nthat it is essential to justify this decision and ensure the robustness of our findings. We selected 18\nbased on a combination of quantitative criteria evaluating HMM performance, including model\nconsistency (i.e., similarity between state time courses), clustering performance (i.e., Calinski-\nHarabasz score), and stability (i.e., free energy), as shown in Fig. S2 and previously described in\nthe Methods section. Notably, prior human34 and rodent35 studies using Gaussian HMM in rsfMRI\nanalyses have also employed larger numbers of states (e.g., 14 or 18).\nc. To further address concerns about potential model underfitting or overfitting, we trained two\nHMMs with lower and higher number of states (i.e., 14 and 22) on the identical rsfMRI dataset.\nRegarding the HMM performance, the 14-state model outperformed our 18-state model in\nclustering performance but showed reduced model consistency and stability (Fig. S2). Conversely,\nthe 22-state model exhibited greater model stability but poorer clustering performance and\nreduced model consistency. We then compared these models with our original 18-state model and\n6\nfound that spatial patterns (Fig. R1-4-1) and characteristics of the states, such as their\nrelatively short lifetimes and high occurrence rates (Fig. R1-4-2), as well as the key transition\nstructures, including the SN-FP and DMN-RL pathways (Fig. R1-4-3), were preserved across\nmodels with different states numbers. We also observed that some spatial patterns of individual\nstates and state transition probabilities showed slight variations, whereby the model with fewer\nstate numbers (e.g., 14) tended to underfit the data by merging distinct spatial patterns and\ntransitions, while the model with more states (e.g., 22) introduced redundant or infrequently\noccurring states, indicating potential overfitting.\nFor instance, state 1 in the 14-state model, which exhibited widespread activity in the\nsensorimotor cortices and striatal regions, resembled a combination of the spatial patterns of states\n1, 2, and 4 from the 18-state model (Fig. R1-4-1C vs. R1-4-1B). Despite this, the SN-FP\ntransition structure was preserved, showing a sequence of transitions from a state exhibiting\npositive sensorimotor cortical activity (state 1), to a state involving activity in the basal forebrain\nand hypothalamus (state 3), and ultimately to states characterized by negative sensorimotor\nactivity (state 2; Fig. R1-4-3B). Similarly, the DMN-RL in the 14-state model retained its loop\nstructure involving robust and widespread activity in DMN-related regions and sensorimotor\ncortices. In the 22-state model, SN-FP and DMN-RL state transition structures remained\ndiscernable in state transition space. However, more states exhibited low occurrence rates (< 17\ntimes/hr) or long lifetime (> 6 s) and had to be excluded from further analysis (states 16-22; Fig.\nR1-4-2C). Additionally, one extra background state (state 13) emerged, displaying similar spatial\npatterns of another background state (state 11; Fig. R1-4-1D). Therefore, our 18-state model\nstrikes a balance between data under- and overfitting by capturing distinct and meaningful rsfMRI\nnetwork dynamics (e.g., state 1  2 reflecting the hierarchical somatosensory activity propagation,\nwhich was not captured in the 14-state model) without introducing excessive redundant or long-\nlasting/infrequently occurring states, as observed in the 22-state model.\nd. In this revision, as discussed in R1-4a above, we added brief descriptions of the quantitative\nevaluation of HMM model performance to the Results section to better highlight the rationale of\nselecting 18 states. Note that justifications were previously only made in the Methods section.\nIn this revision, we added the additional analyses discussed in R1-4c above to the Results section.\nWe also integrated Figs. R1-4-1, R1-4-2, and R1-4-3 as new Figs. S5 and S6 in the revised\nSupplementary Information. These results collectively demonstrated that models with lower or\nhigher state numbers might lead to potential underfitting or overfitting. Moreover, our main\nfindings of rsfMRI network dynamics, including the state transition structures, spatial patterns,\nand characteristics, remained consistent across models with different state numbers, confirming\n7\nthe robustness of our HMM model and the reliability of our results.\nFig. R1-4-1. New 14-state and 22-state HMM models: Comparison of state spatial patterns.\n(A) Definition of atlas-based ROI. Voxel-wise mean activation maps of transient rsfMRI network\nstates showed similar patterns across model with different numbers of states, i.e., (B; also shown\nin Supplementary Figure S1C) original and models with (C) 14 and (D) 22 states.\nThe bar charts showed that states denoted by orange, green, and purple bars exhibited relatively\ntransient characteristics with high occurrence rates. In contrast, states represented by grey bars,\nwhich showed relatively long lifetimes, low occurrence rates, and/or low fractional occupancy,\nwere excluded from subsequent analyses.\nstructures. Transition space (Right) captured 90th percentile of state transitions (Left). Our results\nshowed that the properties of states and transition structures showing in 18-state model (A; also\nshown in Figure 5A and 5B) were largely preserved in (B) 14-state and (C) 22-state model,\ndemonstrating the robustness of our Gaussian HMM model.\nredundancy and interpretability. With the use of initial 18 states, it seems worthwhile to confirm that\nthese are functionally distinct by computing similarity (e.g., Pearson correlation) between state activation\nmaps and cluster similar states into meta-states. This would help identify whether there is functional\noverlap among states and whether timepoints could be reassigned more parsimoniously. Alternatively,\ndescribing all 18 states in detail may dilute the clarity of the manuscript. One potential strategy\nis to focus on the dominant or most functionally relevant states during stimulation and non-stimulation conditions, and\naddress which states exhibit significantly higher lifetime, fractional occupancy, or occurrence rates under each condition.\na. Thank you for raising this concern. Although 18 rsfMRI network states were identified, four (i.e.,\nstates 15–18) were discarded due to relatively long lifetimes, high fractional occupancies, and/or\nlow occurrence rates.\nAs suggested, we conducted additional Pearson correlation analysis across the 14 rsfMRI network\nstate activation maps to evaluate their spatial similarity. We found high spatial similarities (i.e.,\nabsolute value of correlation coefficient, |CC| > 0.5) within our three defined state transition\nstructures of SN-FP, DMN-RL and background states (Fig. R1-5). For instance, states 2 and 5\n(SN-FP; |CC| = 0.8361), states 7–10 (DMN-RL; |CC| > 0.6893), states 11 and 12 (background\nstates; |CC| = 0.5867), and states 13 and 14 (background states; |CC| = 0.5194) shared dominant\nBOLD activation patterns. However, spatial similarities were generally high between SN-FP and\nDMN-RL as well, demonstrating that such metrics alone is likely insufficient to determine\nwhether rsfMRI network states were functionally distinct from one another.\nb. We then conducted hierarchical clustering analysis (see also response to the reviewer’s query in\nR-11 below) to determine if we can better distinguish the rsfMRI networks functionally. Here,\nstates that closely follow one another in a temporal sequence will be clustered together, which can\nreflect neural information flow in brain networks during resting-state. In brief, rsfMRI network\nstates were clustered into three state transition structures, namely SN-FP (states 1–5), DMN-RL\n(states 6–10), and background states (states 11–14; see Fig. R1-11 below), confirming our\nexisting assignment of rsfMRI network state transition structures (Figs. 2, 4 and 5).\nc. In this revision, we decided not to integrate discussions in R1-5a above and Fig. R1-5 below in\nthe manuscript as the classification of transition structures (SN-FP, DMN-RL, and background\nstates) based on their temporal characteristics of and activation maps, discussed in R1-5b above,\nwas more effective in clustering rsfMRI network states. Note that the Results and Methods\nsections were revised to reflect the addition of hierarchical clustering analysis.\nsection to improve clarity and reduce potential redundancy.\nIn this revision, we also consolidated two sentences describing the spatial patterns of background\nstates (states 11-14) in the Results section into a single brief statement to summarize their key\nactivity features. This is done to better focus the study’s description of states 1–10, as they\nconstitute the SN-FP and DMN-RL transition structures, whose probabilities and occurrence rates\nwere significantly altered by single-pulse optogenetic stimulation.\nFig. R1-5. New analysis of rsfMRI network states: Some rsfMRI network states exhibited\nhigh spatial similarity between them. (A) The matrix of correlation coefficient between state\nmean activation maps after spatial smoothing (median filtering with 2 × 2 neighborhood) showed\nthe spatial similarity between rsfMRI network states. (B) Voxel-wise mean activation maps of\nstates with high pairwise spatial similarity (absolute values of correlation coefficient, |CC| > 0.5),\nshowed overlapping or functionally related activities patterns across certain states.\nR1-6 While it is helpful to display spatial patterns graphically in the main figures, I believe it would also\nbe useful to show state transition probabilities in a matrix format in the main figures, as done in some\nsupplemental figures (e.g., S3).\na. Thank you for your suggestion. We have already presented the state transition probability matrices\nin Fig. 5A of the original manuscript, where we compared the transition probability changes\n12\nbetween no-stimulation (OG−) and stimulation (OG+) conditions. Displaying the state transition\nspace and spatiotemporal evolution in Figs. 2C (OG−) and 4D (OG+) was to visually integrate\nboth state transitions and the spatial patterns of substates. This approach provides an intuitive\nunderstanding of how network states evolve spatially and temporally.\nb. In this revision, to better guide readers, we added the following clarification to the captions of\nFigs. 2 and 4: \"The corresponding state transition matrix is shown in Fig. 5A.\".\nR1-7 The definition of “DMN-related” states could benefit from more clarification. It appears that several\nbackground states also strongly resemble DMN features.\na. Thank you for this comment. We agree that the definition of “DMN-related” states warrants\nfurther clarification.\nb. In our study, “DMN-related” refers to rsfMRI network states whose activation patterns were\nconfined within a subset of key regions commonly associated with the DMN, such as auditory\n(Aud), insular (Ins), temporal association (TeA), cingulate (Cg), prelimbic (PrL), retrosplenial\n(RSC), and orbitofrontal cortices (OFC)36.\nWe acknowledge that the spatial patterns of the background states (11-14) exhibited DMN and/or\nhippocampus (HP) activity. However, their transition characteristics were different from those of\nthe DMN-RL states (6-10). Specifically, the background states have no dominant inputs from any\nrsfMRI network and they primarily transit to specific SN-FP states. In contrast, DMN-RL states\nexhibit a loop structure with high sequential transition probabilities, indicating a more stable and\npersistent transition structure. This distinction in transition structure may reflect functional\ndifferences in the type of internally-oriented information processing supported, as discussed\npreviously in the Discussion section. For instance, DMN-RL states might be associated with\ninternally-oriented information processing and global activity involving sustained internal\ncoordination37, while background states may play a role in utilizing existing memory and\ninternally-oriented information to guide and support externally-oriented sensory information\nprocessing37-39. While both DMN-RL states (6–10) and background states (11–14) exhibited\nDMN-related activity, they were classified separately to distinguish their transition characteristics\n(see also our detailed response in R1-5b above and R1-11 below).\nc. In the revision, we clarified the definition of “DMN-related” in the Results and Discussion\nsections and replaced all “DMN regions” with “DMN-related regions.”\nIn this revision, we also explicitly stated in the Results section that our classification of transition\nstructures (SN-FP, DMN-RL, and background states) was based on both state activation maps and\n13\ntransition characteristics.\nR1-8 In Line 215, states 6 and 8 are identified as occurring between two major networks and their\ntransition probabilities can be extracted. However, I am not sure this necessarily means that these brain\nregions “facilitate” switching. I would like to hear more about the authors’ view on this point.\na. We agree that the term \"facilitate switching\" may not best describe the transitions between DMN-\nRL and SN-FP. In our state space (Fig. 5B), we found that states 6 and 8 were located at the\nintersection of these two major transition structures. The observation that these states appeared to\nmediate/modulate transitions (i.e., act as an intermediary or bridge) does not necessarily imply\nthat they actively facilitate switching.\nRather, states 6 and 8 likely represent key transitional points where rsfMRI networks switch\nbetween internally-oriented and externally-oriented processing. Previous studies have shown that\nbasal forebrain (BF) activity, whose robust presence was noted in states 6 and 8, can modulate\nboth the DMN20,40 and SN41,42, indicating that these two rsfMRI network states can easily mediate\ntransitions towards DMN-RL or SN-FP. However, our study does not provide direct evidence that\nthe neuromodulatory regional activities prominent in states 6 and 8 are capable of actively\nfacilitating network structure switching.\nb. In this revision, we avoided using the term “facilitate” to describe the role of states 6 and 8 in\nrsfMRI network state transitions between DMN-RL and SN-FP in the Results and Discussion\nsections. Instead, we used “mediate” to describe the role of states 6 and 8 in switching between\nDMN (task negative) and SN (task positive) states.\nR1-9 Since the authors concatenated all data for analysis, do we know whether certain individual subjects\ncontributed more dominantly to the results? That is, how does individual variability impact the reported\nfindings?\na. Thank you for bringing our attention to the potential impact of individual variability on our\nreported findings. We acknowledge that our existing quantifications of animal variability (see\ndetailed description in R1-9c below) are likely insufficient.\nb. To further address reviewer’s concern, we trained new HMM models with the same parameters\n(e.g., state number = 18 and PCA retaining 70% variance) on different datasets. Here, each dataset\nexcluded one animal from the original, generating 17 and 11 different dataset combinations for OG-\nand OG+ conditions respectively. We then computed the pairwise similarity of the state time\ncourses between the newly trained HMM and the original from the manuscript.\nThe results showed that removing any single animal from the OG− or OG+ group did not\n14\nsubstantially alter the model's output (Fig. R1-9). The similarities between the retrained models and\nthe original model remained high and consistent, with most (i.e., ~65%) falling within mean ± 1\nstandard deviation, indicating Gaussian distribution and that no single subject skewed the model's\noutcome.\nFig. R1-9. New HMM models trained on modified datasets: The output of the HMM model\nwas not skewed by individual animals. Each dot represents the similarity (%) between the state\ntime courses of the original model trained on the full dataset (OG−: n = 17; OG+: n = 11) and a\nmodel trained on a modified dataset in which one animal was excluded from the full dataset. Error\nbar indicates ± 1 standard deviation (SD).\nc. As suggested also by Reviewer 3 to clarify the variability in our samples (see detailed responses\nto R3-4 and the corresponding Fig. R3-4 below), we further calculated the effect sizes (Cohen’s d)\nand 95% confidence intervals (CIs) for state probability changes following stimulation, compared\nto baseline. In brief, we showed that all SN-FP and DMN-RL states (states 1-10) exhibited large\neffect sizes (Cohen’s d > 0.8), indicating consistent alterations to their state probability profiles\nupon stimulation (see Fig. R3-4B below). Furthermore, the 95% CIs for these changes were\nrelatively narrow, indicating robust and consistent findings across animals (see Fig. R3-4A below).\nd. Note that we did examine the consistency of the original model's results across animals via several\ngroup-level analyses and presentation of individual subject data. For example:\ni) For voxel-wise state activation maps, we performed a group-level one-sample t-test with\nGaussian random field correction (GRF, voxel-level p < 0.05, followed by cluster-level p <\n0.001), demonstrating that the observed activation patterns were present across the majority of\nsubjects (Figs. 2C and 4D).\nii) We displayed the state properties (e.g., lifetime, occurrence rate, fractional occupancy; Figs.\n2A and 4C) and transition characteristics (e.g., occurrence rate and duration of transitions for\nSN-FP and DMN-RL; Fig. 5C) of each animal. These values were evenly distributed around\ntheir means, without any apparent skewing in the quantitation of these parameters.\n15\niii) We statistically compared state probabilities following single-pulse stimulation against (group-\nlevel paired t-test with FDR correction) to demonstrate significant and consistent alterations\nacross animals (Fig. 4B).\niv) For the BOLD activation maps (Fig. 3B) elicited by single-pulse optogenetic stimulation,\nthey represent overall activations across all animals and do not reflect inter-subject variability.\ne. In this revision, we incorporated the discussions in R1-9b above in the Results section to\ndemonstrate minimal impact of individual animal variability on the output of the HMM model.\nWe added Fig. R1-9 above as a new Fig. S9 in the Supplementary Information. Methodological\ndescription of the new analysis here was also incorporated in the Methods section.\nIn this revision, discussions in R1-9c above were incorporated into the Results section with\ncorresponding updates to the presentation of state probability profiles and alterations following\nstimulation in Fig. 4A and 4B.\nR1-10 Brain state transitions are ever-changing (though following probabilistic patterns as shown). Does\nthe brain state at the time of thalamic stimulation matter for the outcome? Could the timing of stimulation\nwith respect to current state dynamics influence the resulting reconfiguration?\na. To address this, we examined whether stimulation applied during different states affected their\nsubsequent transitions.\nFor each of the 14 states active at the time of single-pulse VPM stimulation, we quantified the\ncumulative occurrences of transitions to its immediate subsequent states (Fig. R1-10A). To then\ncompute the probabilities of these transitions, we normalized the number of transitions from each\nstate by the total number of stimulation events that occurred during said state, yielding a transition\nprobability matrix. We then constructed a difference matrix by subtracting the baseline (OG−) state\ntransition matrix from the matrix upon single-pulse stimulation (OG+), allowing us to isolate\ninstantaneous, stimulation-induced changes in state transition probabilities (Fig. R1-10B).\nb. Transition probabilities were altered immediately upon stimulation, with most showing increased\ntransition probabilities into states 1 and 2, except for states 1 and 11, whose transition probabilities\ninto state 2 decreased (see the difference matrix in Fig. R1-10B below). These alterations likely led\nto the documented instantaneous increase in probabilities of states 1 and 2 (Fig. 4A), and the\noccurrence rate of the SN-FP structure (Fig. 5C). Conversely, transitions into and within DMN-RL\nstates markedly decreased (transition probability change > 20%; i.e., state 4  6, state 8  9,\nstate 9  10, state 10  6). These decreased transitions likely resulted in the disruption of DMN-RL\nstructure (i.e., reduced transition probabilities to DMN-RL states; Fig. 5A), followed by diversion\n16\ntoward SN-FP states (i.e., increased transition probabilities to states 1 and 2; Fig. 5A).\nTaken together, rsfMRI network state transitions were more likely to rapidly shift toward states 1\nand 2, initiating the SN-FP structure regardless of the specific ongoing state at the time of\nstimulation. Notably, when stimulation occurred during an active DMN-RL state (i.e., 6/7/8/9/10),\nsubsequent transitions to maintain the loop structure were suppressed, indicating a switch from\ninternally-oriented activity to an externally-oriented sensory information processing.\nFig. R1-10. New analysis to examine the influence of current state dynamics at the time of\nstimulation on rsfMRI network reconfiguration: The transitions into states 1 and 2 were\nfacilitated, while transitions into and within DMN-RL were suppressed, regardless of the\nongoing state at the time of thalamic single-pulse stimulation. (A) Cumulative occurrences of\neach state at the time of thalamic stimulation (Left), color-coded by the identity of the immediate\nsubsequent states (1–14). States 1 and 2 were the most frequently transited state compared to\nothers. Transition probabilities following single-pulse stimulation (Right) were computed by\nnormalizing the number of transitions from each state at the time of stimulation by the total\nnumber of stimulation events occurring during said state. (B) State transition matrices of rsfMRI\ndata in the absence of (Left; OG-) and following (Middle; OG+) single-pulse optogenetic\n17\nstimulation, including their differences (Right) to delineate altered transition probabilities at the\ntime of single-pulse stimulation with respect to the current state dynamics.\nc. In this revision, we incorporated the new analyses and discussions on the timing of optogenetic\nstimulation relative to ongoing rsfMRI network states in the Results and Methods sections, as\ndiscussed in R1-10a and R1-10b above. We added Fig. R1-10 above as a new Fig. S10 to the\nSupplementary Information.\nR1-11 The term “hierarchical” is used to describe transitions, yet no formal analysis of hierarchy is\npresented. It would be better to avoid this unless further analysis is provided.\na. Thank you for pointing this out. We agree that the use of the term \"hierarchical\" to describe state\ntransitions should be supported by formal analysis. In response, we have conducted additional\nhierarchical clustering analyses to quantitatively characterize the organization of state transitions.\nb. We performed Ward’s hierarchical clustering43 on a state time course correlation matrix, which\nwas derived from cross-correlation analysis between each pair of state time courses (maximum lag\n= 5 TRs; Fig. R1-11A). For each state pair, we selected the maximum correlation coefficient as an\nindex of temporal similarity. Here, higher values indicate that two states closely followed each\nother, reflecting strong temporal coupling. We revealed that states within each transition\nstructure (i.e., SN-FP, DMN-RL, and background states) clustered together, demonstrating\nthat they were more likely to transition among each other (Fig. R1-11A). An exception was\nstate 1 in the non-stimulation condition, which clustered with background states. As the initial\nstate of the SN-FP structure, state 1 received no dominant input in resting-state (i.e., transition\nprobability exceeding the 90th percentile; Fig. 5B) akin to background states. As discussed\npreviously in the Discussion section of the manuscript, background states associated with memory\nand internally-oriented information processing can modulate externally-oriented processing37-39.\nc. In addition, we computed pairwise correlations of fractional occupancy (FO) values across\nanimals to assess the temporal variability of states at the group level44 (Fig. R1-11B). If two states\nshowed high FO correlation across subjects, they tended to co-vary in the proportion of time spent\nacross animals. The results again showed that states within SN-FP, DMN-RL, and background\nstates were clustered together across animals, as reflected by higher fractional occupancy\ncorrelations, except state 1 under non-stimulation condition (Fig. R1-11B).\nclustering analyses revealed structured organization of rsfMRI network state transitions. (A)\nWard’s hierarchical clustering (Top) based on maximum cross-correlations among state time\ncourses (Bottom; lag ≤ 5 TRs) showed that states within the SN-FP, DMN-RL, and background\nstates clustered together under both OG− (n = 17) and OG+ (n = 11) conditions. (B) Ward’s\nhierarchical clustering (Top) of state fractional occupancy (FO) correlations (Bottom) across\nanimals showed similar clustering, with states in the same functional structure exhibiting higher\nFO correlations. Our results showed that states within the same transition structures exhibited\nsimilar temporal dynamics.\n19\nd. In this revision, we revised the Results section based on discussions in R1-11b and R1-11c above\nto clarify the use of the term “hierarchical” when describing state transitions. Further, Fig. R1-11\nwas added as a new Fig. S3 in the Supplementary Information and methodological description of\nthe hierarchical analyses were added to the Methods section. The hierarchical clustering analyses\nbased on both temporal dynamics and across-subject occupancy patterns provide formal support\nfor a hierarchical organization of rsfMRI state transitions.\nIn this revision, to avoid confusion with the term “hierarchical,” as it was also used in the Results\nsection to describe neural activity propagation along somatosensory thalamo-cortical pathways,\nwe amended statements throughout the manuscript. Specifically, we changed to “hierarchical\nstructure/pathways” to describe such neural activity propagation initiated upon single-pulse\noptogenetic stimulation of the VPM.\nR1-12 The manuscript describes “efficient” neural processing. It is recommended to provide metric or\nmodel to support the use of this term.\na. Thank you for pointing this out to improve the readability of the manuscript. We did not perform\nany quantitative analysis to define \"efficient\" in the context of neural processing. The primary\nfocus of our manuscript is to reveal and dissect rapid brain-wide rsfMRI network dynamics.\nb. In this revision, we removed any mention of efficient neural processing in the Introduction and\nDiscussion sections, as it does not pertain to the central hypothesis of our study.\nR1-13 Regarding statistical analyses, Figure 5A (Difference Matrix) and 5C appear to involve multiple\ncomparisons. Were these corrected (e.g., FDR)? Please specify the correction method used.\na. Thank you for bringing to our attention the lack of specific information detailed for the statistical\nanalysis of state transition probability matrices. To clarify, the difference matrix shown in Fig. 5A\nwas not subjected to statistical correction. A single transition matrix for each dataset (OG− or\nOG+) was estimated from the model trained on the concatenated data from both conditions (OG-\n& OG+), representing the overall state transitions across all trials and animals (see also our\nresponses in R1-15 below for detailed descriptions of our methodology when training the HMM\nmodel). The difference matrix was then computed by simply subtracting the OG− transition\nmatrix from the OG+ transition matrix. Hence, statistical comparisons made between these two\nconditions would not require correction (e.g., FDR). In the difference matrix, the green boxes\nhighlight the top 10% of state transition probability changes following single-pulse stimulation.\nb. For the analysis of DMN-RL transition structures, we performed two-way ANOVA to compare\ntheir occurrence rates and durations across different starting states and experimental conditions\n20\n(i.e., OG- and OG+; Fig. 5C). Following two-way ANOVA, we conducted post-hoc comparisons\nusing Sidak's test (originally annotated in the caption) which corrects for multiple comparisons\nand controls the family-wise error rate (FWE; the probability of making one or more false\ndiscoveries when performing multiple statistical tests).\nc. In this revision, as discussed in R1-13a above, we updated the caption of Fig. 5A to better define\nthe representation of green boxes in the difference matrix. We also included a more detailed\ndescription regarding how the difference matrix was derived and why statistical correction was not\napplied in the Methods section.\nIn this revision, we amended the caption of Fig. 5C to clarify the correction methods used when\ncomparing state occurrence rates and durations as discussed in R1-13b above.\nR1-14 In Figure 3F–G, the label “OG− (n = 17)” suggests a no-stimulation condition, but the figure content appears inconsistent\nwith that label. Please clarify whether this panel reflects stimulation or baseline.\na. Thank you for drawing our attention to the contents of Fig. 3F & G. To clarify, the data primarily\nused for this analysis was from the OG− condition. Data from OG+ condition was partially used\nfor comparison purposes. For instance:\ni) Fig. 3F (left) and Fig. 3G both represented data only from the OG− condition, specifically\nshowing state activation maps and the occurrence rates of state transitions for OG−.\nii) Fig. 3F (right) showed a spatial correlation between the state activation maps from OG− data\nand the time-lapsed BOLD amplitude maps obtained from OG+ data.\nThus, while the label \"OG− (n = 17)\" in the heading indicates the no-stimulation condition animal\ngroup, the content of the figure involves analysis primarily based on OG− data, with a spatial\ncorrelation performed using data from both OG− and OG+.\nb. In this revision, we updated the figure caption of Fig. 3 to clarify this distinction.\nR1-15 The manuscript does not provide detailed methodology for how transition matrices for OG− and\nOG+ were computed. Were they derived from separate HMM models or the same model? Were\ntransitions computed per animal and averaged, or across the full concatenated data? Adding this\ninformation would improve reproducibility.\na. Thank you for your pointed suggestions to improve manuscript’s readability.\nWe trained the HMM model using concatenated data that combined animals from both the OG−\nand OG+ conditions. This approach allowed us to capture the overall dynamics of the 14 relevant\ntransient state transitions across animals and both conditions in a unified model. Next, the model\nwas applied to two sets of concatenated data (one set that consists of animals from OG- condition\nand the other from OG+) to estimate the respective transition probability matrices.\nDue to the application of a unified model to two datasets with distinct conditions, we were able to\nexamine how a single thalamic neural activity input affected the configurations of the 14 transient\nrsfMRI network states and their transitions. Concatenated rsfMRI data temporally has been\nshown to significantly improve the test-retest reliability of functional connectivity measures45,46.\nWe posit that when concatenated dataset is implemented in HMM models, it enables more robust\nand stable estimation of transition probabilities, which is particularly important when\ncharacterizing complex network dynamics and low-probability transitions.\nb. In this revision, we included additional methodological details discussed in R1-15a above in the\nMethods section. This will improve the reproducibility of our analysis and clarify the process by\nwhich the transition probability matrices were computed.\nR1-16 The authors appear to imply causal inference between electrophysiological and fMRI findings.\nWhile the time-locked display of these results is compelling, I feel additional analysis may be warranted\nto support the claim that low-frequency thalamic oscillations and inter-regional coherence “contribute” to\nfMRI connectivity and reconfiguration - rather than these being merely parallel observations.\na. Thank you for your critique here and suggestions to further strengthen our study’s claims\nregarding the relationship between electrophysiological and rsfMRI findings. We performed\nGranger causality analysis47,48 between gross (i.e., spectrally non-specific, 0.2-45 Hz)\nelectrocorticography (ECoG) power fluctuations and rsfMRI network state time courses to\nexamine the causal influence of local gross neural activity on rsfMRI network dynamics. We then\nconducted a multivariate Granger causality analysis to explore the influence of different band-\nlimited ECoG power and coherence fluctuations (i.e., slow, delta, spindle, and gamma frequency\nbands) on rsfMRI network dynamics and their reconfiguration upon stimulation.\nb. Granger causality analysis: In preprocessing, ECoG data were first processed to remove single-\npulse-evoked neural activities (Fig. S20), and gross ECoG power fluctuations were subsequently\ncomputed as the root mean square (RMS) of the bandpass-filtered signals (0.2–45 Hz), as\ndescribed previously in the Methods section. The resulting ECoG power fluctuations were first\nsmoothed (moving average with window size = 5 data points) and downsampled to the same\nsampling rate as the rsfMRI signals (i.e., one sample per second). The vector autoregression\n(VAR) model was applied between ECoG power fluctuations and state time courses for each\n22\nanimal. To capture the most relevant temporal dependencies, we fitted VAR models with lags\nranging from 1 s to 5 s. The optimal lag was determined based on the lowest akaike information\ncriterion (AIC) value49. We then performed the Granger causality test to assess whether neural\npower fluctuations significantly predicted changes in state time courses. The z-scores obtained\nfrom Granger causality analysis were averaged across animals to obtain group-level z-scores,\nwhich were subsequently converted to group-level p-values. Since we aimed to examine how local\nneural activity synchronizations in Cg and S1 “contribute” to large-scale rsfMRI network\ndynamics, we focused on whether ECoG power fluctuations predicted state time courses, rather\nthan the reverse.\nTable R1-16-1. New analyses to examine whether neural power fluctuations significantly\npredicted changes in state time courses: Local gross neural activity synchronizations showed\nsignificant Granger causality effect on rsfMRI network state time courses. Tables showed the\ngroup-level uncorrected p-values derived from Granger causality tests between ECoG power\nfluctuations in Cg (n = 6) and S1 (n = 4) and state time courses in the absence (A; OG-) and\npresence (B; OG+) of single-pulse optogenetic stimulation. The results indicated that preceding\nlocal neural activity synchronizations mostly predicted rsfMRI network dynamics.\nWe found that local gross ECoG power fluctuations in Cg and S1 were significant Granger-\ncauses of state time courses (Cg: states 4, 6-8; S1: states 2, 4, 6-10), irrespective of whether\noptogenetic stimulation was applied (see Table R1-16-1 above). Note that the Granger causality\neffect was more significant in most DMN-RL states (states 6-8), showing lower p-values (p-value\n< 0.0001) than SN-FP states (states 1-5). This finding indicated that local neural activity\nsynchronizations in Cg of the DMN and S1 of the SN can drive rsfMRI network dynamics,\n23\nespecially states associated with internally-oriented information processing.\nc. Multivariate Granger causality analysis: Here, analysis was conducted on ECoG signals such\nas band-limited ECoG power (in local Cg and S1) and inter-trial coherence (inter-regional\nbetween Cg-S1), and the state time courses. Note that single-pulse-evoked neural activities were\nremoved from the ECoG data (Fig. S20). First, to examine whether neural activity\nsynchronizations at different frequency bands influence rsfMRI network dynamics, all data were\naligned to the onset of each state. 24-second windows centered on each state onset (12 s before\nand after) were extracted and averaged within each animal. Second, to examine how these\nsynchronizations contribute to neural activity-triggered reconfiguration, data were aligned to the\noptogenetic single-pulse stimulation onset. Here, 14-second windows centered on each\nstimulation onset were extracted and averaged within each animal. Following the respective data\npre-processing steps, the resulting ECoG signals were subsequently smoothed and downsampled\nto match the temporal resolution of the rsfMRI signals. We then fitted multivariate VAR models\nto the band-limited ECoG power/coherence signals and state time courses. Finally, Granger\ncausality analysis was performed to assess whether ECoG power and inter-regional coherence in\ndifferent frequency bands could significantly predict the temporal dynamics of rsfMRI networks\nand network reconfiguration.\nIn the absence of single thalamic input, local neural synchronizations in both low- (< 15 Hz, e.g.,\nslow, delta, and spindle oscillations) and high-frequency bands (> 15 Hz, e.g., gamma oscillations)\nwere Granger-causes of the majority of state time courses (see Table R1-16-2 below). Low- and\nhigh-frequency power in Cg predicted the occurrence of states 7, 9, and 10, while in S1, these\nfrequency components predicted state 6. Some states were only predicted by local low-frequency\npower in Cg and S1 (Cg: states 2-4; S1: states 2, 4, 7, 8, and 10), particularly delta and/or spindle\noscillations, or by high-frequency gamma oscillations in Cg (states 1 and 6). Meanwhile, long-\nrange synchronization between Cg and S1 only showed significant Granger causality effect on the\ntime courses of states 3 and 5 at low frequency bands. Following single somatosensory thalamic\ninput, we observed that stimulus-induced alterations in both local and long-range low-frequency\nneural synchronizations were Granger-causes of the reconfiguration of rsfMRI network state time\ncourses (see Table R1-16-3 below). Specifically, local low-frequency power in Cg predicted\nalterations to states 1, 7, and 8, while in S1 it predicted states 2 and 5-7. Inter-regional coherence\nbetween Cg and S1 predicted states 1 and 2, especially at slow oscillations (0.2-1 Hz).\ndifferent frequency bands influence rsfMRI network dynamics: Local wide-band and long-\nrange narrow-band (i.e., low-frequency) neural activity synchronizations were significant\nGranger-causes of rsfMRI network dynamics in the absence of single thalamic input. The\ntable showed p-values derived from Granger causality tests between state-aligned ECoG power in\nCg and S1 and Cg-S1 coherence, and state time courses in the absence of single-pulse optogenetic\nstimulation (OG-; p-values were FDR-corrected across frequency bands).\nfrequency bands contribute to activity-triggered reconfiguration of rsfMRI networks: The\nalterations of local and long-range low-frequency neural activity synchronizations were\nsignificant Granger-causes of the reconfiguration of rsfMRI networks upon single thalamic\ninput. The table showed p-values derived from Granger causality tests between stimulus-aligned\nECoG power in Cg and S1 and Cg-S1 coherence, and state time courses (OG+; p-values were\nFDR-corrected across frequency bands).\nTogether, these findings provided statistical support for the contribution of (i) local wide-band\nand long-range narrow-band (i.e., low-frequency) neural activity synchronizations on\nrsfMRI network dynamics, and (ii) local and long-range low-frequency synchronizations on\nnetwork reconfiguration, especially slow oscillations in long-range synchronization.\nd. In this revision, we incorporated the main findings of the Granger causality analysis discussed in\nR1-16b and R1-16c above into the Results and Discussion sections. In parallel, Tables R1-16-1,\nR1-16-2, and R1-16-3 were added as new Fig. S32, Fig. S17, and Fig. S24, respectively, to\nSupplementary Information. These additions further support our claim that local and long-range\n26\nneural activity synchronizations not only co-occur with but also Granger-causally contribute to\nrsfMRI network dynamics and reconfiguration.\nIn this revision, the corresponding methodological details of the Granger causality analyses\ndescribed in R1-16b and R1-16c above have been added to the Methods section.\nR1-17 When comparing EEG power and state probabilities, the small delay between the peak\nelectrophysiological power changes and state timing appears to vary by state. Is there any potential\nexplanation for this? For example, could it be due to differences in the hemodynamic response function\nacross brain regions?\na. Thank you for drawing our attention to this matter. Indeed, while ECoG power fluctuations were\nconsistently time-locked to rsfMRI state onsets, we observed that the timing of peak power of\nelectrophysiological activity relative to state onset varied across states (Fig. 7A and Fig. S26).\nThis variation was particularly evident in DMN-RL states (states 6–10), where we quantified this\nand found significant correlations between the peak latencies (both positive and negative) of gross\nECoG power fluctuations and the state transition sequence within the DMN-RL structure (6  7\n 8  9  10; Fig. 7B and Fig. S20D). In the Discussion section, we have suggested that the\ntiming differences of ECoG signals observed in DMN-RL states might relate to phase differences\nin underlying global signals as measured with rsfMRI50.\nb. Moreover, ECoG power changes in the same local region (e.g., Cg or S1) could be shaped by the\ninformation flow originating from functionally connected brain regions across multiple rsfMRI\nnetworks51,52. Therefore, when a region is active across multiple states, variations in its afferent\ninputs may lead to distinct temporal profiles. For example, in state 1 (i.e., SN-FP), S1 activity\nmay be driven by thalamic input and modulated by RSC and Cg activities, whereas in DMN-RL\nstates, S1 likely receives modulatory inputs from widespread brain regions such as the basal\nforebrain and striatal regions (Fig. 2C). Additionally, the local neural activity synchronization in\nCg across the various rsfMRI network states can be shaped by different inputs (e.g., thalamus or\nbasal forebrain), and thus its electrophysiological response may also vary in timing.\nTogether, the differences in rsfMRI network state timing relative to ECoG power fluctuations are\nlikely influenced by state-dependent differences in neural information flow at the recorded cortical\nregions. Although hemodynamic response function has been shown to differ across regions53,54,\nhere it is likely caused by such variation in local neural activities during a specific state.\nc. In this revision, we incorporated our discussions in R1-17b above into the Discussion section to\nspeculate on the observed varied delays between ECoG power fluctuations and state onsets.\nR2-G I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nWe are grateful for your inputs and the time you dedicated to co-reviewing our manuscript.\nReviewer #3\nR3-G The study combines optogenetics, fMRI, and electrophysiology to investigate resting-state network\ndynamics, which is technically impressive. The discovery that single thalamic inputs can rapidly\nreconfigure resting-state networks (~2s) adds new insight on thalamic inputs modulating the brain\nintrinsic networks. The use of Hidden Markov Modeling (HMM) to characterize transient network states\nand transitions is methodologically robust. Potentially the findings on how thalamic inputs modulate\nDMN and sensorimotor networks could inform research on disorders like schizophrenia or Alzheimer’s.\nThis is a technically outstanding study with novel insights into rapid network dynamics. With revisions\naddressing the above points, it would make a strong contribution to Nature Neuroscience. I recommend\nmajor revision.\nWe thank the reviewer for the constructive comments and for sharing your views on the strength\nand importance of our study. Please see our detailed responses to the points you raised below.\nMajor concerns and suggestions:\nR3-1 The study uses light isoflurane anesthesia, which may alter brain intrinsic networks and neuronal\nactivity compared to awake states. I would suggest discussing how anesthesia might affect the observed\nnetwork reconfiguration (e.g., suppressed arousal mechanisms) and compare with prior awake\nrodent/human studies.\na. Thank you for raising the important point regarding the use of light isoflurane anesthesia in our\nstudy. In fact, we did discuss the potential influence of anesthesia on rsfMRI network dynamics,\nincluding the comparison with prior awake rodent and human studies, in the second paragraph of\nthe Discussion section in the original manuscript. However, our discussions of light isoflurane\neffects were not sufficiently comprehensive and did not extend to suppressed arousal mechanisms.\nb. Key features of rsfMRI network dynamics such as hierarchically organized transitions between\nrsfMRI network states encompassing sensory and motor-related networks, as well as higher-order\n28\ncognition-related networks (Fig. 2B) were similarly observed in previous awake rat55 and\nhuman44,55 studies. Moreover, rsfMRI networks generated through conventional seed-based\nanalysis (Fig. S1B) exhibited similar interhemispheric connectivity in sensory cortical, DMN-\nrelated (e.g., Cg), limbic, and striatal regions to those observed in awake56 or sleeping57 rodents.\nHowever, we fully acknowledge that light anesthesia will inevitably influence such network\ndynamics. States with robust activity in the basal forebrain (BF) and hypothalamus (HTh), two\nregions closely associated with arousal58-60, such as states 6 and 8 (Fig. 2C and 4D) stand to be\nthe most profoundly affected as BF neural activity is attenuated under anesthesia conditions61.\nThese could include alterations in the occurrence rates and certain transition probabilities of\nnetwork states due to differences in arousal levels in rodent62 and human63 brains. Note that states\n6 and 8 were positioned at the intersection of two major transition structures, i.e., sensorimotor\nnetwork-driven feedforward pathways (SN-FP) and default mode network-dominant recurrent\nloop (DMN-RL), and appeared to mediate transitions between externally- and internally-oriented\ninformation processing. We posit that rsfMRI network state transitions and reconfiguration in\nawake animals will have more prominent BF engagement during network reconfiguration. This\ncan manifest as stronger BF activity in the rsfMRI network state activation maps and higher\nprobabilities of state switching between SN-FP and DMN-RL structures in response to\nstimulation. These differences highlight the potential influence of arousal levels on large-scale\nnetwork dynamics and underscore the importance of future awake animal studies to further\nelucidate these mechanisms.\nc. In this revision, we acknowledged and further elaborated on how anesthesia could influence\nnetwork states potentially related to arousal by referring to previous studies, particularly the BF\nactivity vs. light anesthesia as discussed in R3-1b above, in the Discussion section.\nWe also explicitly recognized the need for future awake animal study in order to fully interrogate\nrsfMRI dynamics.\nR3-2 (technical concerns) The optogenetic stimulation is elegant but raises questions on stimulation\nrange, please provide histological verification of stimulation sites post-experiment. If any Off-Target\nvirus injection that could be valuable control, or include a control group with non-ChR2 animals or\ndiscuss limitations. A good control design will be very important for unexpected effects, such as\nsurgery/fiber damage effects/laser thermal stimulation and so on since fMRI signals may be so sensitive\nto be changeable to all introduced factors.\na. We acknowledge the reviewer’s concerns and are fully aware of the numerous factors in\noptogenetic stimulation experiments that can affect the BOLD fMRI signal measured.\noptogenetic fMRI studies9,26,33,64 to ensure minimal damage to brain tissues due to fiber\ninsertion and that the fiber tip was located in the targeted stimulation location (i.e., ventral\nposteromedial thalamic nucleus, VPM).\nTo verify the stimulation sites, high-resolution T2-weighted anatomical MRI images of the fiber-\nimplanted animal brains clearly illustrated the perpendicular fiber tract, and the fiber tip location\nat the middle of VPM (see Fig. R3-2-1A below). Based on our previous works9,26,33,64 and that of\nothers65, the light illumination range in the our study was rather small (200 μm and 350 μm at\n50% and 10% of initial light intensity, respectively). Therefore, the illumination of blue light was\nconstrained only in VPM without spreading to neighboring thalamic nuclei. Histological images\nacquired from a representative rat post-optogenetic fMRI experiment (see Fig. R3-2-2 below and\nFig. 1A) also clearly showed the optical fiber tract, which was perpendicular to the surface of the\ncortex and the tip positioned directly above the center of the VPM. Note that the brightness and\ncontrast of the histology image here was enhanced x2.5 to better visualize the optical fiber tract.\nThe immunohistochemical staining of the brain slices also confirmed the colocalization of the\noptogenetic viral vector (i.e., ChR2) in CaMKIIa excitatory thalamocortical VPM neurons (Fig.\n1A) and their corresponding monosynaptic axonal projections to layer IV of the primary\nsomatosensory cortex (S1; Fig. R3-2-2). Altogether, the accurate fiber implantation location,\nsmall light illumination range, and on-target optogenetic viral vector expression\ndemonstrated that we targeted VPM with high precision and specificity.\nFig. R3-2-1. New illustration of fiber implantation precision at ventral posteromedial\nthalamic nucleus (VPM) & new simulation of light-induced thermal deposition in brain\ntissue at fiber tip: Precise targeting of VPM in the thalamus with optical fiber and minimal\nincrease in temperature at fiber tip caused a 10 ms single-pulse blue light. (A) T2-weighted\nanatomical image of 0.125 mm² isotropic (Left) displayed the overall optical fiber implantation\ntract, and the enlarged view (Right) showed the fiber tip's position at the middle of the VPM,\nwhich was delineated by the green line. The fiber tip is indicated by a red arrow, where the fiber\n30\ntract demarcated by continuous darkened pixels are discontinued by lighter pixels. (B)\nTemperature changes predicted by the Pennes bio-heat model66 for 473 nm light from an optical\nfiber (diameter = 450 μm, NA = 0.5, pulse width = 10 ms, light intensity = 40 mW/mm2).\nFig. R3-2-2. New histological verification of optical fiber tract in a representative rat post-\noptogenetic experiment: Post-experiment confocal image showed ChR2-mCherry expression\nin ipsilateral/right VPM excitatory neurons and optical fiber placement. The brightness and\ncontrast of the confocal image were enhanced x2.5 to visualize the optical fiber tract (Left;\nindicated by white arrows), which appeared as a vertically oriented line of with denser cellular\ndistribution compared to surrounding tissue and a subtle indentation at the cortical surface. The\ntract was positioned directly above the center of the VPM.\nb. The rsfMRI functional connectivity maps shown in Fig. S1B, including our previous rsfMRI\nworks with optical fiber implantation33,67, demonstrated comparable homotopic bilateral\nconnectivity as in healthy control animals shown in studies by us68 and others36,69,70, indicating\nthat the surgical procedures for fiber insertion did not disrupt interhemispheric rsfMRI\nconnectivity. In our study, ipsilateral VPM, LGN, dHP, and S1 were regions where tissue could\nbe damaged by the optical fiber during implantation and disrupt local hemodynamic responses (if\nany). However, the single-pulse optogenetic stimulation evoked BOLD activations in these\nregions (Fig. 3B and Fig. S8A), demonstrating minimal damage to brain tissues traversed by the\noptical fiber during insertion. Note that the activations in these regions surrounding the optical\nfiber were not due to light leakage as expression of ChR2 was absent except at VPM, the targeted\nlocation (see Fig. R3-2-2 above). Together, these results confirmed minimal influence of\nsurgery/fiber implantation on interhemispheric rsfMRI connectivity and local BOLD\nresponses to optogenetic stimulation.\nc. To eliminate the possibility of undesired direct visual stimulation, we wrapped the exposed optical\nfiber and fiber cannula with opaque heat shrinkable sleeves as per our previous optogenetic\n31\nstudies9,26,33,64, thereby preventing light leakage. This precautionary measure has been described\npreviously in the Methods section. Furthermore, in our previous study26, local field potential\n(LFP) recordings from control animals (i.e., no optogenetic viral vector expression) subjected to\nthe identical VPM stimulation setup (as shown in Fig. 1B) showed no light-evoked cortical and\nVPM activity, confirming that the widespread cortical responses observed in this study can\nonly be evoked in animals with optogenetic expression in the VPM.\nd. To further examine the thermal effect of light illumination, we modeled the heat spread around the\ntip of the optical fiber using the well-known Pennes bio-heat equation66, which approximated\nthermal spread as a diffusion process. Due to the sparse and short nature of our 10-ms single-pulse\noptogenetic stimulation paradigm at light intensity of 40 mW/mm2, we found that even the highest\nillumination-induced heating was below 0.06°C (see Fig. R3-2-1B above). This temperature\nchange is much lower than the 0.2-2°C heating described in past experimental measurements and\nmodeling calculations following commonly used sustained illumination of brain tissue, which\nhave been shown to influence both neuronal circuit processes and hemodynamic activity66,71,72.\nMoreover, our previous works have shown that light intensity of 40mW/mm2 do not induce\npseudo-BOLD activations and deactivations at the fiber tip location9,26,64. Hence, our single-pulse\nstimulation experiment does not suffer from undesired light-induced thermal effects.\ne. In this revision, we updated Fig. 1B following Fig. R3-2-1A above to better display the implanted\noptical fiber tract and the fiber tip location at VPM. We also added Fig. R3-2-2 above in Fig. 1A\nto provide a histological verification of the optical fiber tract in the post-experiment confocal\nimage. Discussions in R3-2a and R3-2b above were added in the Results and Methods sections to\nclarify the minimal effects of surgery for fiber implantation to the measured BOLD fMRI signals.\nIn this revision, we added the modeling of the heat spread around the tip of optical fiber upon\nsingle-pulse optogenetic stimulation (see Fig. R3-2-1B above) as a new Fig. S33 in\nSupplementary Information. Discussions as in R3-2d above were added in Methods section to\nclarify that the measured BOLD fMRI signals were not affected by thermal heat deposited at tip of\nthe optical fiber.\nR3-3 The study reports wide-band (local) and narrow-band (inter-regional) synchronizations but could\nbetter explain: How do slow oscillations (<1Hz) specifically drive BOLD changes? Does optogenetic\nstimulation directly alter slow oscillations, or are changes secondary to other frequencies? Add cross-\nfrequency coupling analysis or Granger causality to clarify directed influences.\na. Thank you for drawing our attention to issues with clarity regarding the relationship between\nelectrophysiological activity and BOLD rsfMRI changes. As suggested by the reviewer, we\n32\nperformed multivariate Granger causality analyses to assess the directed influence of local (Cg\nand S1 power) and inter-regional (Cg–S1 coherence) electrocorticography (ECoG)\nsynchronizations across distinct frequency bands (i.e., slow, delta, spindle, and gamma) on rsfMRI\nnetwork state time courses, both at rest and upon stimulation (see also our response to R1-16\nabove). In brief, our results revealed that in the absence of stimulation, local wide-band and\nlong-range narrow-band (i.e., low-frequency, < 15 Hz) neural activity synchronizations were\nsignificant Granger-causes of rsfMRI network dynamics (see Table R1-16-2 above).\nFollowing stimulation, alterations in local and inter-regional low-frequency synchronizations,\nparticularly at slow oscillations (0.2–1 Hz), exhibited significant Granger-causal influence on\nnetwork reconfiguration (see Table R1-16-3 above), indicating a role of these slow rhythms in\nshaping BOLD dynamics.\nb. In this revision, as briefly described in R3-3a above and detailed discussions in R1-16c, we\nincorporated the multivariate Granger causality analyses into the Results and Discussion sections\nto clarify how the observed ECoG power fluctuations across different frequency bands contribute\nto rsfMRI network dynamics both in the absence and presence of single-pulse stimulation. In\nparallel, Tables R1-16-2 and R1-16-3 above were added as new Fig. S17 and Fig. S24, respectively, to\nSupplementary Information. These additions strengthen our claim that local and long-range\nneural activity synchronizations particularly at slow oscillations contribute to rsfMRI network dynamics and reconfiguration (Fig. 6).\nIn this revision, the corresponding methodological details of the multivariate Granger causality\nanalyses have been added to the Methods section.\nR3-4 Some key analyses lack detailed statistics (e.g., effect sizes for state probability changes, cluster\ncorrection methods for fMRI). Report exact p-values (not just thresholds) for major findings. Clarify\nmultiple comparison corrections (e.g., FDR vs. cluster-based). Include confidence intervals for temporal\nmetrics (e.g., 2s reconfiguration latency).\na. We appreciate the reviewer’s critique here and the suggestions made to include more detailed\nstatistical analysis and to clarify existing ones.\nb. We calculated and plotted 95% confidence intervals for the rsfMRI network state probability\nprofiles to provide a clearer representation of the variability (see Fig. R3-4A below). The\nconfidence intervals were small relative to the observed alterations in state probabilities upon\nstimulation, particularly for states 1-10, 13, and 14, indicating stable estimates with low variability\nacross animals. Cohen’s d values were also computed to quantify the effect sizes for state\nprobability changes. Specifically, we calculated Cohen’s d by comparing the state probability\n33\nevery second following single-pulse stimulation with the baseline (i.e., probability averaged over a\n3 s window preceding stimulation). Cohen’s d values > 0.8, 0.5 – 0.8, and 0.2 – 0.5 are indicative\nof large, medium, and small effects, respectively73.\nOur quantitation showed that states 1–10, 13, and 14 exhibited large effect sizes, indicating that\nthese states underwent substantial alterations in state probability profiles upon stimulation\ncompared to baseline (see Fig. R3-4B below). Notably, states 1 and 2, and 6, which were the\ninitial states of SN-FP (1 and 2) and DMN-RL (6) transition structures, showed large effect sizes\nas early as 1 second upon stimulation. Importantly, most states that exhibited large effects (i.e.,\nstates 1, 2, 4–7, and 9) did so within 2 seconds, indicating that single thalamic input rapidly\nreconfigure rsfMRI networks. In contrast, states 11 and 12 exhibited medium and small effect\nsize, respectively, indicating weaker changes upon stimulation.\nFig. R3-4. New statistical analysis to demonstrate effect sizes of state probability alterations:\nState probability profiles were rapidly altered by single somatosensory thalamic stimulation,\nwith most states exhibiting large effect sizes compared to baseline within 2 seconds. (A)\nAveraged state probability profiles of rsfMRI data in the presence (Top; OG+, n = 11) and\nabsence (Bottom; OG-, n = 17) of single-pulse optogenetic stimulation. Error bar denotes 95%\nconfidence interval at each time point across animals. (B) Effect sizes between state probability\nprofiles after single-pulse stimulation and baseline period (3 seconds before stimulation). The\nresult revealed that probability profiles of most states exhibited large effects (1–10, 13, and 14;\nlarge effect: Cohen’s d value > 0.8) compared to baseline, with the majority reaching large effects\nwithin 2 seconds (states 1, 2, 4–7, and 9).\nc. For the voxel-wise substate activation maps (Figs. 2C and 4D), we employed Gaussian Random\nField (GRF) correction, applying an initial voxel-wise threshold of p < 0.05, followed by a\ncluster-level threshold of p < 0.001, to identify statistically robust clusters. However, cluster\ncorrection methods were not applied to the BOLD activation maps upon single-pulse stimulation\n(Fig. 3B) or the time-lapsed BOLD amplitude maps (Fig. 3E), for the following reasons.\nFor the BOLD activation maps (Fig. 3B) elicited by single-pulse optogenetic stimulation, group-\nlevel statistical correction was not performed, as these maps were derived from concatenated data\nacross all animals. The coherence map represents overall activations across all animals and does\nnot reflect inter-subject variability, thereby we did not apply group-level statistical correction.\nFor the time-lapsed BOLD amplitude maps (Fig. 3E), significant BOLD signal amplitudes\nfollowing single-pulse stimulation were computed using paired t-test compared to baseline across\nstimulation events. Due to the large number of sample size (3192 optogenetic stimulation events)\nand the large variability of BOLD responses upon sparse single-pulse stimulation (i.e., followed\nlognormal distribution as documented for spontaneous neuronal firing rates74; Fig. 3D), statistical\ntesting yielded no significant voxels after FDR correction. Consequently, we report uncorrected\nresults at a voxel-wise threshold of p < 0.01, without applying cluster-level correction.\nd. For state-aligned ECoG power and coherence spectrograms without stimulation (OG-; Fig. 6A),\nwe applied FDR correction to the significance masks to identify neural synchronizations that were\nsignificantly different from baseline. This correction has been noted in the corresponding figure\ncaptions. However, for the differential power and coherence spectrograms (OG+ - OG-; Fig. 6B\nand 6C), no statistical correction was applied to the masks, as most clusters disappeared following\nFDR correction. Nonetheless, the quantification in Fig. 6C (i.e., the right-most panels) showed\nconsistent trends of altered neural activity synchronizations, such as clusters b, d, f, g, h, i, and j\nduring state 6  7, demonstrating the robustness of these observed changes.\ne. In this revision, we replaced the plot in Fig. 4A with the updated version from Fig. R3-4A above,\nwhich differs only in that the error bars represent 95% confidence intervals instead of standard\nerrors of the mean to convey statistical reliability. Effect sizes shown in Fig. R3-4B was added to\nFig. 4B. We also incorporated the corresponding statement addressing R3-4b above into the\nResults section.\nIn this revision, we reported the exact p-values in Fig. 3G, 5C, 5D, and 6C, rather than using\nthresholded notations (e.g., * for p < 0.05), and included the exact p-values in the caption of Fig.\n7B, to allow for more precise interpretation of statistical significance.\nIn this revision, as per discussions in R3-4c above, we revised the captions of several figures for\nclarity. For instance, we specified that no correction methods were applied for the BOLD\nactivation maps in Fig. 3B and the time-lapsed BOLD amplitude maps in Fig. 3E. The statistical\ntests used in Fig. 3E was added.\nIn this revision, we added statements in the captions of Fig. 6B and 6C regarding the absence of\n35\nstatistical correction as discussed in R3-4d above.\nR3-5 The DMN suppression and SN enhancement are intriguing but could be contextualized better: How\nmight this rapid reconfiguration translate to perception/attention in awake animals? Contrast with task-\nevoked fMRI studies showing similar DMN/SN trade-offs. Link findings to theories of \"network\nswitching\" (e.g., Menon’s triple-network model).\na. Thank you for the comments and pointed suggestions. We agree that the observed DMN\nsuppression and SN enhancement warrant further contextualization (as also pointed by Reviewer 1\nin R1-1 above) in relation to established models of network switching, and network/brain regions\nrelated to attention.\nb. SN and DMN are typically regarded as task-positive and task-negative networks, respectively,\nwith their activity being enhanced and suppressed in response to external sensory stimuli or task\ndemands10-12. Prior studies often can only focus on longer timescales ranging from several seconds\nto minutes when examining such DMN/SN trade-offs due to the specific cognitive or sensory\ntasks involving complex and sustained neural inputs across multiple pathways8,75-79. The limitation\nhere makes it challenging to isolate whether and how rsfMRI networks can reconfigure rapidly in\nresponse to single and brief neural events to support fast information processing. Our study\nrevealed that such network reconfiguration can be rapidly and sensitively triggered by a single\nthalamic input, with alterations observed within 2 s. Specifically, we found that when any of the\nDMN-RL states were active at the time of single optogenetically evoked somatosensory thalamic\ninput (see also detailed response to R1-10 and corresponding Fig. R1-10 above), transitions to\nmaintain the loop structure were rapidly suppressed to enhance those to state 1 or 2 of SN-FP (Fig.\n5A, 5C and 5D). Further, our study found that the switching between SN-FP and DMN-RL\ntransition structures was mediated by states 6 and 8, which were located at the intersection\nbetween these two transition structures (Fig. 5C and 5D). Such network reconfiguration was\ndriven by low-frequency neural activity synchronizations in primary somatosensory (S1) and\ncingulate (Cg) cortices (Fig. 6 and see also the analyses and discussions to examine the\nrelationship between neural activity and rsfMRI network dynamics in R3-3 above).\nc. Notably states 6 and 8 exhibited robust activities in insular cortex (InS) and orbitofrontal cortex\n(OFC; Fig. 2C and Fig. S1C), which are regions implicated in the salience network of\nanesthetized rat13,14, and in BF. Salience network has been linked to bottom-up attention\nmechanism15,80, which is engaged rapidly in the presence of external sensory stimuli16-18.\nMeanwhile, activity of BF excitatory or inhibitory neurons in rodents can engage sensorimotor\nand DMN-like cortical regions to bias behavioral states toward externally-oriented or internally-\n36\noriented processes19-21. These observations suggest that the rapid reconfiguration of rsfMRI\nnetworks could also be associated with arousal or attentional mechanisms. Moreover, Menon’s\ntriple-network model, a widely recognized theoretical framework describing the dynamic\ninterplay among three core brain networks15, proposes that the salience network plays a central\nrole in switching between task-positive networks (e.g., the central executive network in humans or\nsensorimotor network in rodents) and task-negative networks (e.g., DMN), thereby facilitating fast\nresponses to both externally- and internally-oriented information processing demands13,15,81,82.\nBased on Menon’s triple-network model, we postulate that activities in salience network-related\nregions with prominent inputs from the BF, as observed in the states 6 and 8, mediated switching\nbetween task-positive (i.e., SN) and task-negative networks (i.e., DMN), prioritizing prompt\nprocessing of a bottom-up single somatosensory thalamic input.\nd. We would also like to draw the reviewer’s attention to our existing discussion in the second\nparagraph of the Discussion section that consciousness levels can influence rsfMRI network\ndynamics, including state activation maps, occurrence rates, and transition probabilities, although\nthe overall rsfMRI network dynamics were preserved. Arousal and attention are known to enhance\nsensory processing, while sleep and quiet wakefulness support learning and memory83. Therefore,\ncompared to sleep, awake animals may exhibit a more pronounced increase in the probability of\nstates with positive SN activity (e.g., state 1) and higher SN-FP occurrence rates following\nthalamic input, whereas DMN-RL and background states may become more dominant under sleep\nconditions.\ne. In this revision, we have incorporated the link between our findings and Menon’s triple-network\ntheoretical framework as discussed in R3-5b and R3-5c above, respectively, into the Discussion\nsection. Specifically, we highlighted the potential role of states 6 and 8, which exhibited activity\nin BF and salience network-related regions, in mediating the dynamic switching between task-\npositive and task-negative networks.\nR3-6 HMM Parameters: Why were 18 states chosen? Provide model selection criteria (e.g., free energy\nplots).\na. We thank the reviewer for raising this important point. As described previously in the Methods\nsection, the number of states for the Hidden Markov Model (HMM) was determined to be 18\nbased on relatively high model consistency (i.e., quantified as the mean pairwise similarity of state\ntime courses across iterations), clustering performance (i.e., measured by the Calinski-\nHarabasz score), and model stability (i.e., indicated by low free energy). The corresponding plots of these\ncriteria have been shown in Fig. S2. Notably, our choice of 18 states aligns with prior human34\n37\nand rodent35 studies using Gaussian HMM in rsfMRI analyses, which have commonly employed a\nsimilar number of states (e.g., 14 or 18).\nFig. R3-6. New 14-state and 22-state HMM models: Representative state spatial patterns\nderived from the 14- and 22-state models showed merging and splitting of states,\nrespectively, compared to the 18-state model.\nb. To further validate our choice of state numbers, we analyzed HMMs trained with different\nnumbers of states such as 14 and 22 (see also the responses to R1-4 above). We compared the\nresulting state spatial patterns (Fig. R1-4-1), characteristics (Fig. R1-4-2), and state transition\nmatrices (Fig. R1-4-3) across models with different state numbers. When using fewer states (e.g.,\n14), we observed that certain states and their associated dominant transitions were merged into a\nsingle state (e.g., original states 1, 2, and 4 were combined into state 1 in 14-state model; see Fig.\nR3-6 above). Conversely, increasing the number of states (e.g., 22) resulted in the splitting of\nindividual states in 18-state model into multiple states exhibiting activities in similar brain\nregions. For instance, states 12 in the 18-state model were subdivided into states 11 and 13 in the\n22-state model, showing positive activities in hippocampus (see Fig. R3-6 above). Importantly,\nthe key state spatial features, characteristics (i.e., state lifetime, occurrence rate, and\nfractional occupancy), and transition structures were preserved across models albeit with\ndifferent denomination. These findings support the robustness of the 18-state model and the\nreliability of our results.\nc. In this revision, we added descriptions of model selection criteria to the Results section as\ndiscussed in R3-6a above to better highlight the rationale for our state number selection.\nPreviously, such criteria were only mentioned in the Methods section,\nIn this revision, we incorporated this new analysis and clarified the findings, as discussed briefly\nin R3-6b above, in the Results section. Please also refer to our detailed response to R1-4 above.\nR3-7 ECoG Limitations: Surface recordings may miss deep structures (e.g., thalamus). Acknowledge.\na. Thank you for pointing this out. We acknowledge the limitation that surface ECoG recordings\ncannot measure neural activity in deep brain structures, such as the thalamus. In our study, ECoG\nrecordings were only conducted in key cortical regions of the DMN (Cg) and SN (S1).\n38\nb. In this revision, to clearly acknowledge the technical constraint of our ECoG recordings, we added\na statement in the Results section to explicitly note that electrophysiological recordings were\nlargely sensitive to cortical neural activity, specifically in the Cg and S1 regions, and therefore do\nnot provide direct measurements to subcortical structures.\nR3-8 Also, provide the ECoG method details for a quality verification and reproducibility.\na. Thank you for drawing our attention to issues with clarity in our methodological description of\nECoG recordings. Regarding the gradient artifact removal procedure, we applied the published\nfMRI Artifact Slice Template Removal (FASTR) algorithm via FMRIB84, with artifact timing\ndetermined based on slice-timing event locations. To ensure the quality of the denoised ECoG\nsignals, we performed visual inspection to confirm the absence of periodic or consistent gradient-\nrelated noise. In addition, the presence of detectable neural oscillatory events across multiple\nfrequency bands (e.g., slow, delta, spindle, and gamma oscillations) was also used as an essential\nvalidation criterion (Fig. S34).\nb. In this revision, we expanded the Methods section on ECoG data preprocessing to improve\ntransparency and reproducibility. Specifically, we added more detailed information on the MRI\ngradient artifact removal procedure and quality control procedures applied to the ECoG signals\nfollowing gradient artifact removal, as discussed in R3-8a above.\nR3-9 Share code/HMM models for reproducibility.\na. Thank you for the suggestion. We agree that sharing code and models is essential for\ntransparency and scientific rigor.\nAll procedures for training the HMM and extracting state properties, such as state ROI-wise\nactivation maps, state time courses, and transition matrices, were conducted using the publicly\navailable HMM-MAR toolbox (https://github.com/OHBA-analysis/HMM-MAR), which we have\ncited previously in the Methods section.\nb. In this revision, we uploaded custom software codes used for HMM model training and\nestimation, state transition space computation and visualization, state decomposition, and\ngenerating substate voxel-wise activation maps (without statistical thresholding), as well as the\ntrained and estimated HMM models used in our study to a public repository\n(https://github.com/bispmri/Optogenetic_rsfMRI).\nFor statistical inference, the statistical masks for rsfMRI network activation maps were generated\nusing the DPABI toolbox (https://rfmri.org/DPABI). In this revision, we updated the Methods\nsection to include these links and tool references to ensure clarity and reproducibility.\nanalyzing within-state dynamics. The Methods section describes the HMM process well but could\nslightly elaborate on precisely how the boundaries for these substates were defined based on the state\nprobability time courses (e.g., using thresholds on the probability or its derivative?).\na. We appreciate the reviewer’s suggestion to clarify how the boundaries for the entering, plateau,\nand exiting substates were defined in our analysis of within-state dynamics.\nThe boundaries for the entering, plateau, and exiting substates were determined based on\nthresholded probability values and/or slopes of state time courses. Here is a more detailed\nexplanation of how these boundaries were defined:\ni) Entering Substate: This is the period when the state time course increased, indicating that the\nbrain was transitioning into the state. Typically, this is when the slope of the state time course\nbecame positive. In addition to the slope, we imposed a threshold on state probability, where\nthe state probability exceeds the 95th percentile of the state’s own baseline fluctuations but\nremains below a fixed upper bound of probability = 0.9 (see Fig. R3-10 below). This approach\nensured the exclusion of low-probability baseline periods while including sufficient data points\nduring entering the states, allowing the threshold for each state time course to be adaptively\ndefined based on its own baseline fluctuations.\nFig. R3-10. New illustration: State decomposition into substates was based on probability\nand/or slopes of the state time courses. The grey and orange dashed lines in the state\ndecomposition panel (Left) denote the lower and upper thresholds, respectively, used for\nidentifying entering and exiting substates. The plateau period was identified when the state\nprobability exceeded 0.9.\nii) Plateau Substate: This is the period when state probability remained high and stable.\nindicating that the system was fully engaged in the state. This threshold was chosen to\nrepresent a strong and sustained engagement with the state.\niii) Exiting Substate: This is the period when the state probability decreased, indicating that the\nsystem was transitioning out of the state. Here, the slope of the state time course became\nnegative. Similar to the entering substate, probability must be above the 95th percentile of the\nstate’s baseline fluctuations while remaining below the upper limit of 0.9.\nb. In this revision, we clarified the definitions and criteria of the boundaries for identifying entering, plateau,\nwhich were based on the probability thresholds and/or slopes of the state time\ncourses, as discussed in R3-10a above, in the Methods section.\nIn this revision, we updated the state decomposition panel in Fig. 1C as shown in Fig. R3-10\nabove and the corresponding captions to improve the clarity of substate boundary identification.\nR3-11 The Discussion compellingly links the observed network reconfiguration to altered neural\nsynchrony, particularly in low frequencies. Could the Discussion speculate slightly further on the\npotential causal cascade? Does the thalamic input primarily affect sensory cortex synchrony, which then\npropagates, or does it potentially influence key modulatory nodes (like BF/HTh implicated in states 6/8)\nmore directly to orchestrate the shift? The text currently implies synchrony changes lead the\nreconfiguration (lines 426-427), but acknowledging the potential role of hub regions like BF/HTh in\nmediating this could add nuance.\na. We appreciate your suggestion to further consider the causal cascade underlying network\nreconfiguration. Previous anatomical tracing studies have consistently shown that VPM input\nprimarily targets the primary sensorimotor cortices, with subsequent propagation to subcortical\nstructures85-87. This feedforward propagation is also supported by our findings, as demonstrated by\nthe early activation of sensory regions following stimulation (Fig. 3C and Movie S1).\nOur simultaneous rsfMRI and ECoG measurements further demonstrated that single-pulse\nthalamic stimulation induced rapid changes in local and inter-regional coherence across cortical\nregions, particularly at low-frequency bands (Fig. 6). Importantly, the Granger causality analysis\nadded in this revision indicated the contribution of low-frequency cortical synchronizations on\nsubsequent network reconfiguration (see detailed discussions in R1-16 and R3-3 above).\nb. Furthermore, cortical neurons and BF neurons maintain extensive reciprocal connections88,\nenabling bidirectional desynchronization processes that regulate brain state transitions83. BF has\nbeen shown to modulate large-scale resting-state fMRI fluctuations89 (as also discussed in R3-1b\n41\nand R3-5c above). It can coordinate complex network dynamics via cell-type-specific projections\nto regions associated with DMN and sensorimotor networks, thereby influencing internally or\nexternally oriented behavioral preferences20,21. Thus, cortical regions may also provide\ndescending inputs to the BF, which in turn facilitates rapid rsfMRI network reconfiguration.\nHowever, we are unable to discern the exact/primary cascades due to the absence of neural\nrecordings in the subcortical hub regions of BF/HTh.\nc. In this revision, we have added the discussions to speculate on the potential cascade as in R3-11a\nand R3-11b above in the Discussion section to acknowledge the role of BF during activity-\ntriggered rsfMRI network reconfiguration. We propose that input from the VPM is first\ntransmitted to cortical regions, where changes in low-frequency neural activity synchronizations\ncan directly influence reconfiguration of network states. These cortical synchronizations can then\npotentially modulate BF activity, which in turn mediates the switching from DMN to\nsensorimotor networks.\nMinor Suggestions:\nR3-12 Add schematic of hypothesized circuit mechanism (e.g., thalamus → BF → DMN suppression).\na. Thank you for the suggestion here.\nb. This hypothesized mechanism that we first propose involves a single thalamic input inducing\nalterations in S1/Cg neural synchronizations at low-frequency bands, which in turn\ncontribute to DMN suppression and SN activation (i.e., 1: thalamus  cortical activity\nsynchronizations  rsfMRI network reconfiguration; shown in the left panel of Fig. R3-12\nbelow; see also discussions in R3-11 above).\nThe second proposed mechanism is more speculative due to the absence of electrophysiological\nrecordings in BF and subcortical regions. Here, alterations in S1/Cg low-frequency cortical neural\nsynchronizations modulate BF activity, which then contribute to network reconfiguration upon\nsingle somatosensory thalamic input (i.e., 2: thalamus  cortex  BF/HTh  rsfMRI network\nreconfiguration; shown in the right panel of Fig. R3-12 below). Future studies will be required to\ndirectly test this circuit mechanism.\nrsfMRI network reconfiguration upon single somatosensory thalamic input. Single thalamic\ninput altered local Cg and S1, and their long-range neural activity synchronizations at low-\nfrequency bands (Left; < 15 Hz; OG ) underlying rapid rsfMRI network reconfiguration. This\n+\nreconfiguration involved suppression and diversion of internally-oriented information processing\n(i.e., DMN-RL) through a state (state 8) also exhibiting BF and hypothalamus (HTh) activities to\nenhance sensory information processing (i.e., SN-FP). Additionally, such thalamic input can be\nrelayed through the cortex to neuromodulatory regions, such as BF and HTh, to further modulate\nswitching between DMN-RL and SN-FP (Right).\nc. In this revision, we added a schematic illustrating the hypothesized circuit mechanism of rsfMRI\nnetwork reconfiguration (Fig. R3-12 above) as a new figure panel in Fig. 8C.\nR3-13 Improve color contrast in spectrograms (Fig. 6) for readability.\na. Thank you for the suggestion.\nb. In this revision, we improved the readability of the spectrograms in Fig. 6 by modifying the mask\nfor statistically significant regions from green to black, thereby increasing the contrast with the\ncolor of the spectrogram.\nIn this revision, we also updated the spectrograms in Figs. S14, S18, S19, S21-S23, and S25 in\nthe Supplementary Information to reflect identical changes made to the main manuscript figure.\na. Thank you for your comment. RsfMRI networks are thought to be associated with these two\nmodes of neural information processing. Specifically, externally-oriented information processing\nrefers to information flow driven by environmental stimuli and task demands, typically involving\nsensory processing that follows a feedforward flow of information along hierarchical sensory\npathways, where brain regions sequentially encode increasingly higher-order features of the\nexternal environment90. This process predominantly engages the sensorimotor networks. In\ncontrast, internally-oriented processing relies on internal cognitive resources, such as memory\nretrieval, planning, and self-referential thought, and is primarily supported by the default mode\nnetwork (DMN)91,92.\nb. In this revision, we clarified the definitions of “internally- and externally-oriented information\nprocessing”, as discussed in R3-14a above, in first paragraph of the Introduction section.\nR3-15 Avoid overstating \"first to show\" claims—cite prior work on fast fMRI dynamics (e.g., Yu et al.,\n2016).\na. Thanks for drawing our attention to this issue. We have cited the relevant prior work by Yu et\nal29., which provides important context on fast optogenetically-evoked hemodynamic responses\nwithin 400 milliseconds. As clarified in our response to Reviewer 1 in R1-3 above, our use of the\nterm “rapid” specifically refers to the reconfiguration of brain-wide rsfMRI networks in\nresponse to a single and brief neural input, rather than to the localized hemodynamic\nresponse alone. This distinction is important, as typical rsfMRI network dynamics are often\nqualitatively observed on timescales ranging from several tens of seconds to minutes31-33.\nb. In this revision, we revised the Introduction section to clearly define the concept of rapid large-\nscale rsfMRI network dynamics, distinguishing it from localized hemodynamic responses.\nR4-G In this study, authors investigated how brain-wide rsfMRI networks enable rapid and efficient\nneural information flow and processing. By using a single-pulse optogenetic stimulation of the\nsomatosensory thalamus, authors showed activity-triggered rapid reconfiguration of rsfMRI networks on\nthe order of a few seconds (~ 2 s to begin reconfiguring; ~ 10 s to normalize). Combining simultaneous\nelectrophysiology and fMRI measurements, they further revealed that this reconfiguration was driven by\ndistinct synchronizations of neural activity within and between networks occurring across specific time\nand frequency windows. This is a very important study with a tremendous amount of data. The analysis is\ncomplex but rigorous. I only have a few minor comments.\nWe sincerely thank the reviewer for the supportive and encouraging comments. We appreciate\nyour recognition of the significance of our study and the rigor in our analysis. Your positive\nevaluation of our integrative approach, which combines single-pulse optogenetics, resting-state\nfMRI, and electrophysiology to examine rapid brain-wide network reconfiguration, has been\nhighly motivating. We are also grateful for your suggestions, which helped us further clarify and\nstrengthen the manuscript.\nR4-1 The results of the study are largely descriptive, and it might be difficult for readers to clearly grasp\nthe picture of brain state transitions reported herein. I wonder, in addition to laying out all observations,\nwhether authors can consider a way to summarize their findings in a way that might be easier for readers\nfrom outside of the field to understand?\na. Thank you for this helpful suggestion to improve readability of the manuscript to non-experts. We\nacknowledge that the study contains numerous descriptions due to the varied analyses and\ncorresponding quantitation to examine rsfMRI network dynamics (Figs. 2, 3 and 4),\nthe\ntimescales of rapid rsfMRI reconfiguration (Figs 4 and 5) and neural activities that drove such\ndynamics and reconfiguration (Figs. 6 and 7).\nb. In this revision, we decided to include a new schematic (see Fig. R4-1 below) as a new Fig. 8 that\nsummarizes the major findings of our study demonstrating rapid rsfMRI network reconfiguration\nupon a single somatosensory thalamic input. This illustration is intended to provide a clearer and\nmore accessible overview of our main findings. The figure was cited in the appropriate paragraph\nsummarizing our work in the Discussion section.\nnetwork dynamics and the rapid activity-triggered network reconfiguration upon single\nsomatosensory thalamic input. (A) In the absence of stimulation (OG), two dominant state\n- transition structures were identified, i.e., DMN-dominant recurrent loop (DMN-RL) and\nsensorimotor network-driven feedforward pathways (SN-FP; Left). The switching between these\nstructures were mediated by a rsfMRI network state (state 6) characterized by prominent activity\nin the basal forebrain (BF) and hypothalamus (HTh). This process is intrinsically dynamic because\nthe brain must continuously adapt the behaviors to environmental changes, switching between\nbrain states for internally- (DMN-RL) and externally-oriented (SN-FP) information processing90-\n92, as demonstrated by the rapid (~ 2 s) suppression of DMN-RL states to strengthen SN-FP states\nfollowing activation of a single optogenetically evoked somatosensory thalamic input (Right). (B)\nLocal wide-band and long-range narrow-band cortical neural activity synchronizations were found\nto drive the observed rsfMRI network dynamics, where such synchronizations were significant\npredictors of the presence of various rsfMRI network states. (C) Single somatosensory thalamic\ninput subsequently altered local Cg and S1, and their long-range neural activity synchronizations,\nparticularly at low-frequency bands (< 15 Hz; OG ), which triggered network reconfiguration\n+\ninvolving suppression and diversion of internally-oriented information processing (DMN-RL) to\n46\nenhance sensory information processing (SN-FP; Left). Additionally, thalamic input can be\nrelayed through the cortex to neuromodulatory regions, such as BF and HTh88, to further modulate\nswitching between DMN-RL and SN-FP (Right), though this pathway remains speculative due to\nthe absence of electrophysiological recordings at BF/HTh.\nR4-2 How stable are the results at a different probability threshold?\na. Thank you for raising this important point.\nTo evaluate the stability of our results under different transition probability thresholds, we\nexamined the state transition structure using both 95th (see Fig. R4-2A below) and 85th percentiles\nas thresholds (see Fig. R4-2B below). Note that the 90th percentile was used as a threshold in our\noriginal analysis (Fig. 5B). As expected, a higher probability threshold reduced the number of\nidentified state transitions, while a lower threshold increased them. However, across the 95th and\n85th percentile thresholds, the major transition structures remained consistent regardless of\nexperimental conditions.\nSpecifically, the DMN-RL states still exhibited a robust loop-like transition pattern (state 6  7\n 8  9  10  6), and the key SN-FP transition from state 1 to state 2 was preserved in the\nabsence of single thalamic input (OG-; left panels in Fig. R4-2A and R4-2B). As demonstrated in\nthe original manuscript, the DMN-RL was suppressed and diverted to the SN-FP through state 8\nupon single thalamic input (Fig. 5). The state 8  9 transition was no longer observed in the 95th-\npercentile state transition space (right panel in Fig. R4-2A), indicating that its transition\nprobability fell below the top 5% of all transitions upon stimulation. This is consistent with our\nobservation with the original probability threshold, where state 8 showed an appreciable decrease\n(-49%) in transitioning to state 9 (Fig. 5D). However, the state 8  1 transition (right panel in\nFig. R4-2A) was masked, despite a substantial increase (+61%) following stimulation, indicating\ndecreased sensitivity due to an excessively high probability threshold.\nThe background states still exhibited only dominant outputs and received inputs almost evenly\nfrom other states, with the exception of state 11, which received a dominant input from state 14 in\nthe 85th-percentile transition space due to the lowered probability threshold (Fig. R4-2B). Note\nthat in our study, we only focused on SN-FP and DMN-RL rsfMRI network states as they were\nreconfigured upon a single thalamic input (Fig. 4A).\nOverall, our key findings based on the 90th percentile threshold applied to the state transition space\nwere not affected. In fact, it showed a good balance to sensitize the detection of significant\nalterations to state transitions.\nthresholds: State transition spaces of rsfMRI data at different probability thresholds\ndemonstrate preserved major transition structures. (A) and (B) showed the transition space\ncapturing the 95th and 85th percentile of transition probabilities in the OG- (Left) and OG+ (Right)\nconditions, respectively. Both thresholds largely preserved the DMN-RL loop structure and the\nkey SN-FP transition from state 1 to state 2. Notably, the transition probability and occurrence\nrate of state 8  9 were dramatically reduced upon single thalamic input (see Figure 5A and 5D),\nand consequently, this transition was not observed in the 95th-percentile state transition space.\nBackground states exhibited only dominant outputs, with the exception of state 11, which received\ninput from state 14 in the 85th-percentile space. Transitions involving states without dominant\ninput are indicated by black dashed arrows, corresponding to probabilities at or below the 95th or\n85th percentiles.\nb. In this revision, as discussed in R4-2a above, the existing probability threshold using the 90th\npercentile strikes a good balance to detect significant alterations of state transitions. Further, the\nstate transition matrices in Fig. 5A already provided the full set of transition probabilities. We felt\n48\nthat visualizing state transition spaces at different thresholds, though illustrative of pathway\nrobustness, would be redundant. Hence, we decided not to include Fig. R4-2 and the\ncorresponding descriptions in the revised manuscript for simplicity.\nR4-3 The finding that a single somatosensory thalamic input at VPM largely maintains the overall\nstability of brain dynamics is consistent with the findings in a recent study (PMID: 35179203), despite the\ndifferences in the model of stimulations and animal’s physiological condition (awake vs anesthetized).\nAuthors may consider discussing these findings in a larger framework.\na. Thank you for sharing this study and highlighting its relevance. We agree that our findings can be\nfurther contextualized within a broader framework of spontaneous neural activity and its\ninteraction with extraneous inputs (i.e., either sensory stimulus or optogenetic stimulation).\nb. The referenced study (PMID: 35179203)93 demonstrated that brain-wide ongoing rsfMRI BOLD\nactivity, reflected in large-scale co-fluctuations across regions, was largely independent of\nexternal stimulation and contributes significantly to cross-trial variability in evoked responses. As\nnoted by the reviewer, this is consistent with our findings that variability is an intrinsic feature of\nspontaneous network dynamics, and that the single-pulse thalamic input largely preserves, rather\nthan overrides, this variability (Fig. 3D), in contrast to a study with prolonged stimulation that\nshowed a reduction in neural variability25. Moreover, both in awake and anesthetized conditions,\nthese large-scale patterns remained relatively stable in the presence of stimulation or task\ndemands24,93. These observations suggest that the rsfMRI networks are highly robust and\nconserved across different brain states and task conditions.\nc. Importantly, our findings showed that rsfMRI network dynamics mirrored the neural activity\npropagation of stimulus-evoked neural activity (Fig. 3F and 3G) along somatosensory thalamo-\ncortical pathways, suggesting that ongoing rsfMRI BOLD activity fluctuations reflect hierarchical\nstructural pathways that were engaged during tasks. Consistent with prior human and animal\nstudies, spontaneous neural activity has been shown to reflect task-evoked activity22-25. For\nexample, in humans, spontaneous memory reactivations have been observed to propagate across\ndistributed brain-wide networks that closely resemble those engaged during sensory, motor, and\ncognitive tasks23,24. In rodents, resting-state cortical activity patterns have been shown to encode\nand partially predict sensory and movement-related signals, reflecting their integration for\nbehavior24,25. Our results suggest that the stable architecture of resting-state functional\nnetworks allowed existing networks and their transitions to be rapidly and selectively\nrecruited/suppressed in response to task demands or stimulation, thereby facilitating fast\nprocessing of any neural inputs without necessitating the formation of new network\nstructures.\nd. In this revision, we have incorporated the insights in R4-3a and R4-3b above into the Discussion\nsection.\nReferences\n1. Hipp, J.F., Hawellek, D.J., Corbetta, M., Siegel, M. & Engel, A.K. Large-scale cortical correlation\nstructure of spontaneous oscillatory activity. Nature neuroscience 15, 884-890 (2012).\n2. Krishnan, G.P., Gonzalez, O.C. & Bazhenov, M. Origin of slow spontaneous resting-state neuronal\nfluctuations in brain networks. Proc Natl Acad Sci U S A 115, 6858-6863 (2018).\n3. Vafaii, H., et al. Multimodal measures of spontaneous brain activity reveal both common and\ndivergent patterns of cortical functional organization. Nat Commun 15, 229 (2024).\n4. Anderson, J.S., Carandini, M. & Ferster, D. Orientation tuning of input conductance, excitation, and\ninhibition in cat primary visual cortex. Journal of neurophysiology 84, 909-926 (2000).\n5. Okun, M. & Lampl, I. Instantaneous correlation of excitation and inhibition during ongoing and\nsensory-evoked activities. Nat Neurosci 11, 535-537 (2008).\n6. Poo, C. & Isaacson, J.S. Odor representations in olfactory cortex: \"sparse\" coding, global inhibition,\nand oscillations. Neuron 62, 850-861 (2009).\n7. Kline, A.M., Aponte, D.A., Tsukano, H., Giovannucci, A. & Kato, H.K. Inhibitory gating of\ncoincidence-dependent sensory binding in secondary auditory cortex. Nat Commun 12, 4610 (2021).\n8. Chao, T.-H.H., et al. Neuronal dynamics of the default mode network and anterior insular cortex:\nIntrinsic properties and modulation by salient stimuli. Science Advances 9, eade5732 (2023).\n9. Leong, A.T.L., Wang, X., Wong, E.C., Dong, C.M. & Wu, E.X. Neural activity temporal pattern\ndictates long-range propagation targets. Neuroimage 235, 118032 (2021).\n10. Menon, V. 20 years of the default mode network: A review and synthesis. Neuron 111, 2469-2487\n(2023).\n11. Fox, M.D., et al. The human brain is intrinsically organized into dynamic, anticorrelated functional\nnetworks. Proceedings of the National Academy of Sciences 102, 9673-9678 (2005).\n12. Hinz, R., et al. Bottom-up sensory processing can induce negative BOLD responses and reduce\nfunctional connectivity in nodes of the default mode-like network in rats. Neuroimage 197, 167-176\n(2019).\n13. Xu, N., et al. Functional Connectivity of the Brain Across Rodents and Humans. Front Neurosci 16,\n816331 (2022).\n14. Tsai, P.J., et al. Converging Structural and Functional Evidence for a Rat Salience Network. Biol\nPsychiatry 88, 867-878 (2020).\n15. Menon, V. Large-scale brain networks and psychopathology: a unifying triple network model. Trends\nCogn Sci 15, 483-506 (2011).\n16. Baluch, F. & Itti, L. Mechanisms of top-down attention. Trends Neurosci 34, 210-224 (2011).\n17. Bradley, C., Nydam, A.S., Dux, P.E. & Mattingley, J.B. State-dependent effects of neural stimulation\non brain function and cognition. Nat Rev Neurosci (2022).\n18. Condylis, C., et al. Context-Dependent Sensory Processing across Primary and Secondary\nSomatosensory Cortex. Neuron 106, 515-525 e515 (2020).\n19. Sarter, M., Givens, B. & Bruno, J.P. The cognitive neuroscience of sustained attention: where top-\ndown meets bottom-up. Brain research reviews 35, 146-160 (2001).\n20. Zou, Y., et al. Cell-type-specific optogenetic fMRI on basal forebrain reveals functional network\nbasis of behavioral preference. Neuron 112, 1342-1357. e1346 (2024).\n21. Lozano-Montes, L., et al. Optogenetic Stimulation of Basal Forebrain Parvalbumin Neurons Activates\nthe Default Mode Network and Associated Behaviors. Cell Rep 33, 108359 (2020).\n22. Cole, M.W., Bassett, D.S., Power, J.D., Braver, T.S. & Petersen, S.E. Intrinsic and task-evoked\nnetwork architectures of the human brain. Neuron 83, 238-251 (2014).\n23. Huang, Q., et al. Replay-triggered brain-wide activation in humans. Nat Commun 15, 7185 (2024).\n24. Dimakou, A., Pezzulo, G., Zangrossi, A. & Corbetta, M. The predictive nature of spontaneous brain\nactivity across scales and species. Neuron 113, 1310-1332 (2025).\n25. Liu, M., Liang, Y., Song, C., Knopfel, T. & Zhou, C. Cortex-wide spontaneous activity non-linearly\nsteers propagating sensory-evoked activity in awake mice. Cell Rep 41, 111740 (2022).\n26. Leong, A.T., et al. Long-range projections coordinate distributed brain-wide neural activity with a\nspecific spatiotemporal profile. Proc Natl Acad Sci U S A 113, E8306-E8315 (2016).\n27. Chuang, K.H., Li, Z., Huang, H.H., Khorasani Gerdekoohi, S. & Athwal, D. Hemodynamic transient\nand functional connectivity follow structural connectivity and cell type over the brain hierarchy. Proc\nNatl Acad Sci U S A 120, e2202435120 (2023).\n28. Jung, W.B., Im, G.H., Jiang, H. & Kim, S.G. Early fMRI responses to somatosensory and optogenetic\nstimulation reflect neural information flow. Proc Natl Acad Sci U S A 118(2021).\n29. Yu, X., et al. Sensory and optogenetically driven single-vessel fMRI. Nat Methods 13, 337-340\n(2016).\n30. Chen, B.R., Bouchard, M.B., McCaslin, A.F., Burgess, S.A. & Hillman, E.M. High-speed vascular\ndynamics of the hemodynamic response. Neuroimage 54, 1021-1030 (2011).\n31. Bricault, S., et al. Peripheral contributions to resting state brain dynamics. Nat Commun 15, 10820\n(2024).\n32. Menon, V., et al. Optogenetic stimulation of anterior insular cortex neurons in male rats reveals\ncausal mechanisms underlying suppression of the default mode network by the salience network. Nat\nCommun 14, 866 (2023).\n33. Wang, X., Leong, A.T.L., Chan, R.W., Liu, Y. & Wu, E.X. Thalamic low frequency activity\nfacilitates resting-state cortical interhemispheric MRI functional connectivity. Neuroimage 201,\n115985 (2019).\n34. Stevner, A.B.A., et al. Discovery of key whole-brain transitions and dynamics during human\nwakefulness and non-REM sleep. Nat Commun 10, 1035 (2019).\n35. Salvan, P., et al. Frequency modulation of entorhinal cortex neuronal activity drives distinct\nfrequency-dependent states of brain-wide dynamics. Cell Rep 37, 109954 (2021).\n36. Lu, H., et al. Rat brains also have a default mode network. Proc Natl Acad Sci U S A 109, 3979-3984\n(2012).\n37. Kaefer, K., Stella, F., McNaughton, B.L. & Battaglia, F.P. Replay, the default mode network and the\ncascaded memory systems model. Nat Rev Neurosci 23, 628-640 (2022).\n38. Preston, A.R. & Eichenbaum, H. Interplay of hippocampus and prefrontal cortex in memory. Current\nbiology 23, R764-R773 (2013).\n39. Yeshurun, Y., Nguyen, M. & Hasson, U. The default mode network: where the idiosyncratic self\nmeets the shared social world. Nat Rev Neurosci 22, 181-192 (2021).\n40. Nair, J., et al. Basal forebrain contributes to default mode network regulation. Proc Natl Acad Sci U S\nA 115, 1352-1357 (2018).\n41. Zhu, F., Elnozahy, S., Lawlor, J. & Kuchibhotla, K.V. The cholinergic basal forebrain provides a\nparallel channel for state-dependent sensory signaling to auditory cortex. Nat Neurosci 26, 810-819\n(2023).\n42. Ananth, M.R., Rajebhosale, P., Kim, R., Talmage, D.A. & Role, L.W. Basal forebrain cholinergic\nsignalling: development, connectivity and roles in cognition. Nat Rev Neurosci 24, 233-251 (2023).\n43. Ward Jr, J.H. Hierarchical grouping to optimize an objective function. Journal of the American\nstatistical association 58, 236-244 (1963).\n44. Vidaurre, D., Smith, S.M. & Woolrich, M.W. Brain network dynamics are hierarchically organized\nin time. Proc Natl Acad Sci U S A 114, 12827-12832 (2017).\n45. Cho, J.W., Korchmaros, A., Vogelstein, J.T., Milham, M.P. & Xu, T. Impact of concatenating fMRI\ndata on reliability for functional connectomics. Neuroimage 226, 117549 (2021).\n46. Elliott, M.L., et al. General functional connectivity: Shared features of resting-state and task fMRI\ndrive reliable and heritable individual differences in functional brain networks. Neuroimage 189, 516-\n532 (2019).\n47. Granger, C.W. Investigating causal relations by econometric models and cross-spectral methods.\nEconometrica: journal of the Econometric Society, 424-438 (1969).\n48. Hamilton, J.P., Chen, G., Thomason, M.E., Schwartz, M.E. & Gotlib, I.H. Investigating neural\nprimacy in Major Depressive Disorder: multivariate Granger causality analysis of resting-state fMRI\ntime-series data. Mol Psychiatry 16, 763-772 (2011).\n49. Deshpande, G., LaConte, S., James, G.A., Peltier, S. & Hu, X. Multivariate Granger causality analysis\nof fMRI data. Human brain mapping 30, 1361-1373 (2009).\n50. Raut, R.V., et al. Global waves synchronize the brain’s functional systems with fluctuating arousal.\nScience advances 7, eabf2709 (2021).\n51. Cole, M.W., Ito, T., Bassett, D.S. & Schultz, D.H. Activity flow over resting-state networks shapes\ncognitive task activations. Nat Neurosci 19, 1718-1726 (2016).\n52. Ito, T., Hearne, L.J. & Cole, M.W. A cortical hierarchy of localized and distributed processes revealed\nvia dissociation of task activations, connectivity changes, and intrinsic timescales. Neuroimage 221,\n117141 (2020).\n53. Aguirre, G.K., Zarahn, E. & D'Esposito, M. The variability of human, BOLD hemodynamic\nresponses. Neuroimage 8, 360-369 (1998).\n54. Handwerker, D.A., Ollinger, J.M. & D'Esposito, M. Variation of BOLD hemodynamic responses\nacross subjects and brain regions and their effects on statistical analyses. Neuroimage 21, 1639-1651\n(2004).\n55. Ma, Z. & Zhang, N. Temporal transitions of spontaneous brain activity. Elife 7(2018).\n56. Liu, Y., et al. An open database of resting-state fMRI in awake rats. Neuroimage 220, 117094 (2020).\n57. Yu, Y., et al. Sleep fMRI with simultaneous electrophysiology at 9.4 T in male mice. Nature\nCommunications 14(2023).\n58. Liu, X., et al. Subcortical evidence for a contribution of arousal to fMRI studies of brain activity.\nNat Commun 9, 395 (2018).\n59. Liu, X., Leopold, D.A. & Yang, Y. Single-neuron firing cascades underlie global spontaneous brain\nevents. Proc Natl Acad Sci U S A 118(2021).\n60. Yamagata, T., et al. The hypothalamic link between arousal and sleep homeostasis in mice. Proc Natl\nAcad Sci U S A 118(2021).\n61. Thomas, K., et al. GABAergic neurons in basal forebrain exert frequency-specific modulation on\nauditory cortex and enhance attentional selection of auditory stimuli. Commun Biol 8, 149 (2025).\n62. Gutierrez-Barragan, D., et al. Unique spatiotemporal fMRI dynamics in the awake mouse brain. Curr\nBiol 32, 631-644 e636 (2022).\n63. Huang, Z., Zhang, J., Wu, J., Mashour, G.A. & Hudetz, A.G. Temporal circuit of macroscale dynamic\nbrain activity supports human consciousness. Science advances 6, eaaz0087 (2020).\n64. Wang, X., et al. Functional MRI reveals brain-wide actions of thalamically-initiated oscillatory\nactivities on associative memory consolidation. Nat Commun 14, 2195 (2023).\n65. Yizhar, O., Fenno, L.E., Davidson, T.J., Mogri, M. & Deisseroth, K. Optogenetics in neural systems.\nNeuron 71, 9-34 (2011).\n66. Stujenske, J.M., Spellman, T. & Gordon, J.A. Modeling the Spatiotemporal Dynamics of Light and\nHeat Propagation for In Vivo Optogenetics. Cell Rep 12, 525-534 (2015).\n67. Chan, R.W., et al. Low-frequency hippocampal-cortical activity drives brain-wide resting-state\nfunctional MRI connectivity. Proc Natl Acad Sci U S A 114, E6972-E6981 (2017).\n68. Zhou, I.Y., et al. Brain resting-state functional MRI connectivity: morphological foundation and\nplasticity. Neuroimage 84, 1-10 (2014).\n69. Chuang, K.H. & Nasrallah, F.A. Functional networks and network perturbations in rodents.\nNeuroimage 163, 419-436 (2017).\n70. Liang, Z., Liu, X. & Zhang, N. Dynamic resting state functional connectivity in awake and\nanesthetized rodents. Neuroimage 104, 89-99 (2015).\n71. Arias-Gil, G., Ohl, F.W., Takagaki, K. & Lippert, M.T. Measurement, modeling, and prediction of\ntemperature rise due to optogenetic brain stimulation. Neurophotonics 3, 045007 (2016).\n72. Owen, S.F., Liu, M.H. & Kreitzer, A.C. Thermal constraints on in vivo optogenetic manipulations.\nNat Neurosci 22, 1061-1065 (2019).\n73. Sullivan, G.M. & Feinn, R. Using Effect Size-or Why the P Value Is Not Enough. J Grad Med Educ\n4, 279-282 (2012).\n74. Buzsáki, G. & Mizuseki, K. The log-dynamic brain: how skewed distributions affect network\noperations. Nature Reviews Neuroscience 15, 264-278 (2014).\n75. Vidaurre, D., Llera, A., Smith, S.M. & Woolrich, M.W. Behavioural relevance of spontaneous,\ntransient brain network interactions in fMRI. Neuroimage 229, 117713 (2021).\n76. Song, H., Park, B.Y., Park, H. & Shim, W.M. Cognitive and Neural State Dynamics of Narrative\nComprehension. J Neurosci 41, 8972-8990 (2021).\n77. Song, H., Shim, W.M. & Rosenberg, M.D. Large-scale neural dynamics in a shared low-dimensional\nstate space reflect cognitive and attentional dynamics. Elife 12(2023).\n78. Vidaurre, D., et al. Discovering dynamic brain networks from big data in rest and task. Neuroimage\n180, 646-656 (2018).\n79. Taghia, J., et al. Uncovering hidden brain state dynamics that regulate performance and decision-\nmaking during cognition. Nat Commun 9, 2505 (2018).\n80. Nghiem, T.E., et al. Space wandering in the rodent default mode network. Proc Natl Acad Sci U S A\n121, e2315167121 (2024).\n81. Swick, D., Ashley, V. & Turken, U. Are the neural correlates of stopping and not going identical?\nQuantitative meta-analysis of two response inhibition tasks. Neuroimage 56, 1655-1665 (2011).\n82. Chen, T., Cai, W., Ryali, S., Supekar, K. & Menon, V. Distinct Global Brain Dynamics and\nSpatiotemporal Organization of the Salience Network. PLoS Biol 14, e1002469 (2016).\n83. Lee, S.H. & Dan, Y. Neuromodulation of brain states. Neuron 76, 209-222 (2012).\n84. Allen, P.J., Josephs, O. & Turner, R. A method for removing imaging artifact from continuous EEG\nrecorded during functional MRI. Neuroimage 12, 230-239 (2000).\n85. Jung, W.B., Jiang, H., Lee, S. & Kim, S.G. Dissection of brain-wide resting-state and functional\nsomatosensory circuits by fMRI with optogenetic silencing. Proc Natl Acad Sci U S A 119(2022).\n86. Shepherd, G.M.G. & Yamawaki, N. Untangling the cortico-thalamo-cortical loop: cellular pieces of a\nknotty circuit puzzle. Nat Rev Neurosci 22, 389-406 (2021).\n87. Theyel, B.B., Llano, D.A. & Sherman, S.M. The corticothalamocortical circuit drives higher-order\ncortex in the mouse. Nat Neurosci 13, 84-88 (2010).\n88. Tian, J., et al. Dissection of the long-range projections of specific neurons at the synaptic level in the\nwhole mouse brain. Proc Natl Acad Sci U S A 119, e2202536119 (2022).\n89. Turchi, J., et al. The Basal Forebrain Regulates Global Resting-State fMRI Fluctuations. Neuron 97,\n940-952 e944 (2018).\n90. Parker, P.R.L., Brown, M.A., Smear, M.C. & Niell, C.M. Movement-Related Signals in Sensory\nAreas: Roles in Natural Behavior. Trends Neurosci 43, 581-595 (2020).\n91. Cohen, D., Nakai, T. & Nishimoto, S. Brain networks are decoupled from external stimuli during\ninternal cognition. Neuroimage 256, 119230 (2022).\n92. Anticevic, A., et al. The role of default network deactivation in cognition and disease. Trends Cogn\nSci 16, 584-592 (2012).\n93. Zhang, Q., et al. Brain-wide ongoing activity is responsible for significant cross-trial BOLD\nvariability. Cereb Cortex 32, 5311-5329 (2022).\nRESPONSES AND REVISION OF NATURE COMMUNICATIONS\nMANUSCRIPT # NCOMMS-25-18143A\nBrain-wide Resting-state fMRI Network Dynamics Elicited by Activation of Single Thalamic Input\nReviewer #1\nR1-G The authors have done a fantastic job revising the manuscript. I have no further questions.\nWe thank the reviewer for his/her compliment and prior suggestions to improve the manuscript.\nReviewer #2\nR2-G I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nWe are grateful for your inputs and the time you dedicated to co-reviewing our manuscript.\nReviewer #3\nR3-G The revised manuscript \"Brain-wide Resting-state fMRI Network Dynamics Elicited by Activation\nof Single Thalamic Input\" by Xie et al. addresses critical questions about the temporal dynamics and neural\nmechanisms of resting-state fMRI (rsfMRI) networks. The authors have made substantial revisions to the\noriginal submission, significantly strengthening the study's methodology, clarity, and impact. Below, I\ndetail my assessment of the revisions and provide recommendations for final improvements.\nMajor Strengths of the Revision:\nMethodological Rigor and Transparency:\nThe authors have clarified the Hidden Markov Modeling (HMM) approach, including state decomposition\ninto substates (entering/plateau/exiting) and the rationale for selecting 18 states (Fig. S2). This transparency\nenhances reproducibility.\nThe simultaneous rsfMRI-electrocorticography (ECoG) experiments (Figs. 6–7) now robustly link neural\nactivity synchronizations (intra-regional wide-band; inter-regional narrow-band) to rsfMRI dynamics. The\n1\ninclusion of permutation tests and FDR correction strengthens statistical claims.\nAddressing Temporal Resolution Concerns:\nThe study convincingly demonstrates rapid rsfMRI reconfiguration (~2 s onset, ~10 s normalization) using\nsparse single-pulse optogenetics. The coherence analysis of 3192 stimulation events (Fig. 3) and temporal\ndownsampling (0.1 s resolution) mitigate limitations of EPI sequence timing (TR = 1 s).\nNeural Mechanism Elucidation:\nThe revised manuscript clarifies how single thalamic inputs suppress internally oriented DMN-RL states\nwhile boosting externally oriented SN-FP states via neuromodulatory hubs (BF/HTh) (Figs. 4–5). The\nECoG data (Fig. 6) elegantly show that low-frequency synchronizations (<15 Hz) drive reconfiguration.\nSuggestion: While the added discussion on BF modulation under anesthesia is excellent, explicitly contrast\nthe ~2s reconfiguration latency in this study with awake rodent/human fMRI dynamics (e.g., Gutierrez-\nBarragan et al., Curr Biol 2022; Song et al., eLife 2023). Highlight whether anesthesia slows network\nswitching compared to awake states.\nRationale: This contextualizes the \"rapid\" reconfiguration claim and clarifies anesthesia-specific\nlimitations.\nAnesthesia Considerations:\nThe authors now explicitly address anesthesia effects (Fig. S8), showing neural oscillation densities (~1%\nisoflurane) comparable to natural sleep. This contextualizes findings within awake/sleep literature (Lu et\nal., 2012; Yu et al., 2023).\nOptogenetic Controls:\nSuggestion: Add histological verification of non-target regions (e.g., adjacent thalamic nuclei) in the\nsupplement to confirm absence of off-target ChR2 expression.\nRationale: The thermal modeling (Fig. S33) and fiber placement validation are robust, but direct histology\nwould preempt concerns about viral spread.\nOverall Recommendation:\nAccept after minor revisions. The manuscript now provides compelling evidence that rsfMRI networks can\nrapidly reconfigure to support neural information processing, with mechanistic insights into local/global\nneural synchronizations. The revisions have adequately addressed prior concerns, and the study advances\nthe field by bridging optogenetics, fMRI, and electrophysiology.\nWe sincerely thank the reviewer for the constructive feedback and are encouraged by the reviewer’s\ncomments on our revisions to address his/her concerns. We also appreciate the valuable suggestions\n2\non contextualizing the network reconfiguration latency under anesthesia and adding histological\nverification, to which we have provided detailed responses below.\nRemaining Concerns and Recommendations:\nR3-1 While the added discussion on BF modulation under anesthesia is excellent, explicitly contrast\nthe ~2s reconfiguration latency in this study with awake rodent/human fMRI dynamics (e.g., Gutierrez-\nBarragan et al., Curr Biol 2022; Song et al., eLife 2023). Highlight whether anesthesia slows network\nswitching compared to awake states.\nRationale: This contextualizes the \"rapid\" reconfiguration claim and clarifies anesthesia-specific\nlimitations.\na. Thank you for the suggestions to contrast our rsfMRI network reconfiguration latency with findings\nfrom other fMRI studies. Isolating and quantifying the latency of network responses to a single\nneural input has remained challenging in previous studies, largely due to the use of rsfMRI\nmeasurements alone1 or the use of sensory/cognitive tasks that evoked prolonged and complex\ninformation flows across multiple parallel pathways2-6. For example, the study by Gutierrez-\nBarragan et al., 2022, which was conducted without external stimulation or task engagement, did\nnot explicitly report the latency of network reconfiguration1. The study by Song et al., 2023, in\nawake humans identified a slower reconfiguration timescale (~3–6 s following narrative event\nboundaries) than ours (~2 s following a single neural input), likely due to the complexity of the\ntasks, which induced neural activities arising from diverse sources, including cognitive and\nattentional fluctuations2. In this study, we were able to directly quantify the latency of rsfMRI\nnetwork reconfiguration by employing a single-pulse optogenetic stimulation paradigm targeting\nthe excitatory neurons of ventral posteromedial thalamic nucleus (VPM), thereby isolating a well-\ndefined somatosensory thalamo-cortical pathway and avoiding the complex, long-lasting neural\ncascades induced by natural sensory stimuli/tasks6,7.\nb. Notably, the influence of anesthesia on the speed of rsfMRI network reconfiguration has not been\ndirectly quantified in previous rodent/human studies. Nevertheless, some evidence postulated that\nanesthesia might slow the large-scale network reconfiguration by reducing the number of state\ntransitions per session8,9. This reduction suggests decreased flexibility of network switching and a\ndiminished capability of networks to rapidly respond and process neural inputs. However, direct\nevidence regarding the influence of anesthesia on the latency of activity-triggered rsfMRI network\nreconfiguration is lacking, highlighting the need for future studies to address this issue.\nc. In this revision, we explicitly clarified in the Discussion the challenges faced by prior rsfMRI\n3\nstudies in isolating and quantifying the latency of rapid network reconfiguration in response to\nspecific neural inputs, and we cited the work of Song et al., 2023 in the corresponding sentence, as\ndiscussed in R3-1a above. The study by Gutierrez-Barragan et al., 2022 has already been cited in\nthe second paragraph of original Discussion section to discuss the influence of anesthesia on rsfMRI\nnetwork dynamics.\nIn this revision, as discussed in R3-1b above, we contextualized in the Discussion section that there\nis currently no direct evidence demonstrating the effect of anesthesia on the speed of rsfMRI\nnetwork reconfiguration, and emphasized the need for future studies to examine this effect.\nR3-2 (concern of optogenetic controls) Suggestion: Add histological verification of non-target regions\n(e.g., adjacent thalamic nuclei) in the supplement to confirm absence of off-target ChR2 expression.\nRationale: The thermal modeling (Fig. S33) and fiber placement validation are robust, but direct histology\nwould preempt concerns about viral spread.\na. Thank you for the suggestion. We enlarged the histological image shown in the left panel of Fig.\n1A (see Fig. R3-2 below). ChR2 expressions were identified in the adjacent thalamic nuclei,\nincluding the posterior nucleus (PO) and ventral posterolateral nucleus (VPL), due to the expected\ndiffusion of the injected viral vector from the centre of VPM7,10-12. Although viral spread extended\nto the neighboring nuclei, accurate fiber placement ensured that our stimulation was centered on\nVPM (Fig. 1A and 1B). Importantly, this accurate targeting was further supported by the\noptogenetically-evoked BOLD responses at VPM, which were consistent with our previous\nstudies7,10, exhibiting the most robust BOLD responses in monosynaptic targets of VPM13, i.e., the\nipsilateral primary somatosensory barrel field (S1BF; Fig. 3B). In contrast to the downstream\nactivation targets of VPM, stimulation of PO predominantly activated ipsilateral S1BF and bilateral\nvisual subcortical regions (lateral geniculate nucleus, LGN and superior colliculus, SC) as shown\nin our previous work10. Meanwhile, VPL primarily conveyed bottom-up somatosensory input to the\nipsilateral primary somatosensory cortex of limb (S1Limb)14. Taken together, the accurate fiber\nplacement and overall downstream activation targets identified by BOLD fMRI support that our\nkey findings were specific to VPM stimulation.\nb. In this revision, as discussed in R3-2a above, we clarified that the injection target was made in the\nmiddle of VPM in the Methods section.\nIn this revision, as also discussed in R3-2a above, we added a statement in the Results section to\nhighlight that, in addition to accurate fiber placement, the BOLD activation maps provided further\nsupport that our stimulation setup was specific to VPM.\nFig. R3-2. New illustration of viral expression in VPM and its adjacent thalamic nuclei: ChR2\nexpression in ventral posteromedial nucleus (VPM) and neighboring thalamic nuclei.\nHistological image showed that the viral injection targeted the VPM, as it was located in the middle\nof the viral spread, with some additional spread into the posterior nucleus (PO) and ventral\nposterolateral nucleus (VPL) due to diffusion of the injected solution.\nR3-3 (concern of temporal resolution and neural dynamics) While the 0.1 s downsampling improves\ntemporal granularity, BOLD signals inherently reflect hemodynamic delays. The claim that rsfMRI\nresponses \"could be actually shorter\" than 2 s (p. 26) requires caution.\nRecommendation: Explicitly discuss the disconnect between fast neural events (e.g., gamma oscillations)\nand slow BOLD responses. Cite Logothetis (2012) to frame hemodynamic limitations.\na. Thank you for raising this important point. Note that we downsampled the BOLD responses into\n0.1 s intervals to enhance the temporal resolution for examining the delay of single-pulse-evoked\nneural activity propagation (Fig. 3C and 3E; Movie S1). However, BOLD fMRI data used for\nrsfMRI network analyses, such as HMM training and prediction, were kept at the acquired temporal\nresolution (i.e., repetition time, TR = 1 s). We fully agree that BOLD signals provide an indirect\nmeasure of neural activity through blood oxygenation changes, inherently constrained by\nhemodynamic delays, and thus cannot directly capture the fast dynamics of neural events15,16.\nIn the original Discussion section, our statement that rsfMRI responses “could be actually shorter”\nthan 2 s was intended to highlight the limitations of the temporal resolution of our imaging protocol\n(TR = 1 s) rather than to suggest equivalence with neural timescales. Most conventional fMRI\nacquisitions are on the order of seconds (e.g., with TR typically ranging from 1 to 5 s), resulting in\nthe BOLD response lagging behind neural oscillations by several seconds (≥1 TR)17-20. For\n5\nexample, in Logothetis et al., 201217, the lag between gamma events and localized peak event-\ntriggered BOLD responses is approximately 4 s with an EPI sequence parameter of TR = 2 s. In\ncontrast, a study employing high temporal resolution fMRI acquisition protocols (e.g., 100 ms) has\nreported onset delays of localized sensory-evoked or optogenetically-evoked hemodynamic\nresponses as short as 400 ms21. These studies suggest that while the rsfMRI network\nreconfiguration observed in our study (~ 2 s) is inevitably constrained by hemodynamic delays,\nusing higher temporal resolution fMRI acquisition may enable the detection of more rapid\nrsfMRI network responses.\nb. In this revision, we clarified the hemodynamic limitation on the rsfMRI network response time of\nour study in the Discussion section, as discussed in R3-3a above, and cited the referenced study to\nsupport this limitation.\nR3-4 (concern of anesthesia and \"resting state\" generalizability) Despite physiological stability (Fig.\nS8), isoflurane may alter DMN dynamics (Liang et al., 2015). The term \"resting-state\" should be qualified.\nRecommendation: Add a limitation paragraph acknowledging anesthesia-induced state differences versus\nawake resting states (e.g., Gutierrez-Barragan et al., 2022).\na. In our study, rsfMRI data were acquired under light isoflurane anesthesia (~1%) without task\nengagement or external stimulation, a condition typically used in rodents to prevent motion and\nminimize stress22. As noted in the original Discussion section, consciousness levels in rodents and\nhumans affected the occurrence rates and transition probabilities of rsfMRI network states, and\nanesthesia was shown to alter spatial patterns of states associated with cognition.\nBeyond these acknowledged effects, anesthesia was also reported to alter the lifetime of certain\nrsfMRI network states9. Moreover, animal studies, including that by Liang et al., 2015, showed that\nanesthesia weakened the mean activation of some rsfMRI network states associated with cognition,\nyet overall network activity can still be reliably detected23,24. For example, isoflurane was shown to\nattenuate the infralimbic cortex-related activities in rodents1,24. However, the effect of isoflurane on\nstates related to sensorimotor systems was not pronounced24.\nImportantly, the spatial patterns of rsfMRI network states and their overall transition structures has\nbeen shown to be largely preserved in anesthetized conditions. Specifically, our rsfMRI networks\n(Fig. S1B) exhibited interhemispheric connectivity in sensory cortical, DMN-related, limbic, and\nstriatal regions that were similar to those observed in awake23-25 or sleeping26 rodents. The overall\nrsfMRI network state transition structure remained largely conserved across awake and anesthetized\nconditions, although anesthesia may alter certain transition probabilities, as demonstrated by\nGutierrez-Barragan et al., 20221. Notably, our findings aligned with key features of rsfMRI network\n6\ndynamics reported in awake rat27 and human27,28 studies, including hierarchically organized\ntransitions between network states involving sensory and higher-order cognitive systems (Fig. 2B).\nFor detailed discussion, please refer to the second paragraph of original Discussion section.\nTogether, despite potential influences of anesthesia, the preservation of spatial patterns23,24,28 and\nhierarchical transition structures28 in rsfMRI networks indicates the generalizability of our findings\nbeyond anesthetized conditions. In particular, as rsfMRI networks associated with sensorimotor\nfunctions are less susceptible to anesthesia24, our network reconfiguration supporting fast\nprocessing of somatosensory thalamic input is likely generalizable beyond anesthetized conditions.\nb. In this revision, we separated the statement regarding the effects of anesthesia on the properties of\nrsfMRI network states in the second paragraph of the Discussion section into an independent\nsentence to emphasize the limitations imposed by anesthesia.\nIn this revision, we incorporated more specific descriptions in the Discussion section, such as the\nimpact of anesthesia on state lifetime and weakened mean activation of DMN-related state, as\ndiscussed in R3-4a above.\nIn this revision, we added a statement in the Discussion section to acknowledge the potential\ninfluence of anesthesia while noting the generalizability of our findings beyond anesthetized\nconditions, as discussed in R3-4a above.\nR3-5 (concern of overinterpretation of \"bifurcation role\") States 6/8 are proposed as \"bifurcation\npoints\" between DMN-RL and SN-FP (Fig. 5D). While their transition shifts are clear, causal evidence for\nneuromodulatory \"bifurcation\" is indirect.\nRecommendation: Tone down mechanistic language (e.g., \"suggesting a role\" vs. \"demonstrating\"). Cite\nNair et al. (2018) on BF-DMN modulation to support inference.\na. Thank you for the suggestion.\nb. In this revision, we modified the original phrase “highlighting their bifurcation role” in the Results\nsection to “suggesting their bifurcation role” to avoid overly mechanistic language.\nThe study by Nair et al., 2018, along with other relevant literature, was previously cited in the\noriginal Discussion section to support the modulatory effects of BF on both DMN29,30 and SN31,32.\nR3-6 (concern of ECoG limitations) Explicitly state that subcortical hubs (BF/HTh) were not recorded\nand that their role in the cascade remains inferred. Since thalamic input may directly affect BF (state 6/8\nhub), suggesting future experiments with simultaneous subcortical recordings would strengthen the causal\ncascade hypothesis.\na. Thank you for this suggestion.\nb. In this revision, we added a statement in the Discussion section to explicitly clarify that\nneuromodulatory regional activities, such as BF and HTh, were not directly recorded in our study.\nFuture investigations incorporating subcortical recordings may help elucidate the causal role of\nBF/HTh in the rapid rsfMRI network reconfiguration.\nR3-7 (concern of schematic clarity) In Fig. 8C, use distinct colors for feedforward (thalamocortical) vs.\nfeedback (cortico-BF) pathways.\nRationale: Clarifies speculative vs. supported mechanisms in the circuit diagram.\na. Thanks for your suggestion to improve clarity.\nb. In this revision, we modified Fig. 8C by changing the solid arrow representing the cortico-BF\npathway to a dashed arrow and explicitly indicated it as a speculative pathway in the corresponding\nfigure caption, in contrast to the relatively well-supported thalamocortical pathway.\nR3-8 (concern of statistical reporting) Figure 5C uses non-parametric tests (Mann-Whitney) but omits\neffect sizes. Outlier removal (ROUT, Q=10%) may bias occurrence-rate comparisons.\nThe statistics improvements are solid - effect sizes and CIs address variability concerns. But the cluster\ncorrection justification for Fig 3E remains problematic. They should either apply more sensitive methods\nlike TFCE or clearly state it's uncorrected as a study limitation.\nRecommendation: Report effect sizes (e.g., Cohen’s *d*) for key comparisons. Justify outlier thresholds in\nmethods.\na. Thank you for drawing our attention to this matter.\nTo avoid the potential bias introduced by the outlier removal procedure (Fig. 5C and 5D), we\nperformed additional statistical analyses using rsfMRI data without outlier removal as well as with\na lower threshold (i.e., robust regression and outlier removal method, ROUT; Q = 5%). These\nanalyses showed that our key results remained statistically significant across all conditions, although\nthe level of significance was reduced in most comparisons without outlier removal (i.e., larger p-\nvalues; see Fig. R3-8 below). Specifically, state 6 remained the main initiating state of completed\ndefault mode network-dominant recurrent loops (DMN-RLs) at resting state, exhibiting\nsignificantly higher occurrence rates than DMN-RLs initiated from states 7-10 across all outlier\nremoval thresholds. Moreover, single thalamic input robustly and significantly reduced the\noccurrence rate of state 6-initiated DMN-RL transition, while increasing the occurrence rates of\n8\nsensorimotor network-driven feedforward pathways (SN-FP) transition and state 8  1 transition.\nb. To complement the statistical comparisons, we also calculated their effect sizes (see Fig. R3-8\nbelow). These analyses revealed large effect sizes (Cohen’s d ≥ 0.8) in all key comparisons, despite\nreduced Cohen’s d values when no outlier removal was applied.\nFig. R3-8. New statistical analysis to examine occurrence rates of state transitions under\ndifferent outlier removal conditions with corresponding effect sizes: Statistical significance\nand large effect sizes of state transition occurrence rates were largely preserved across\ndifferent outlier removal thresholds. Occurrence rates of state transitions without (OG-) and with\n(OG+) a single thalamic input were shown under three conditions: (A) no outlier removal,\n(B) Q =\n5% outlier removal, and (C) Q = 10% outlier removal (n = 17 in OG- and n = 11 in OG+). In all\ncases, state 6 remained the main initiating state of DMN-RLs at resting state, while single thalamic\ninput reduced state 6–initiated DMN-RL transitions, and increased SN-FP and state 8  1 transition.\nStatistical comparisons of DMN-RL transition structures were conducted using two-way ANOVA\nwith FDR correction (Left). Comparisons of SN-FP transitions and transitions starting from state 6\n9\nand 8 between OG- and OG+ conditions were performed using the nonparametric Mann–Whitney\nU test (Middle and Right). Error bars indicate ± SEM; p and d represent p-values and Cohen’s d\nvalues, respectively.\nc. We would also like to draw the reviewer’s attention to our existing description in the caption of Fig.\n3E that the time-lapsed BOLD amplitude maps were presented uncorrected. This decision was based\non the large number of stimulation events (3192) and the large variability of BOLD responses upon\nsparse single-pulse stimulation, which followed lognormal distribution as documented for\nspontaneous neuronal firing rates33 (Fig. 3D), resulting in no significant voxels surviving cluster-\nlevel correction.\nd. In this revision, we replaced Figs. 5C and 5D with Fig. R3-8A above, which presents results\nderived from rsfMRI data without outlier removal, in order to avoid bias in the comparisons of\noccurrence rates introduced by the outlier removal method.\nIn this revision, we added effect sizes for the significant results presented in Fig. R3-8A above to\nFig. 5C and 5D to provide a more comprehensive statistical description of the group differences.\nIn this revision, as discussed in R3-8c above, we added a sentence in the Methods section to explain\nthat the maps in Fig. 3E are uncorrected, as this is due to our experimental design.\nMinor Suggestions:\nR3-9 Figure 3E–F: Clarify how spatial similarity was computed (e.g., Pearson correlation vs. spatial\noverlap).\na. Thank you for the suggestion. We quantified spatial similarity using the two-dimensional Pearson\ncorrelation coefficient (MATLAB function: corr2) between the state activation map (Fig. 3F) and\nthe time-lapsed BOLD amplitude map (Fig. 3E and Movie S1) after spatial smoothing (median\nfiltering with a 2 × 2 neighborhood).\nb. In this revision, we clarified that we quantified the spatial similarity using 2-D spatial Pearson\ncorrelation in the caption of Fig. 3F and in the Methods section, as discussed in R3-9a above.\nR3-10 Methods (ECoG): Specify band-pass filter cutoffs for \"gross\" power (0.2–45 Hz).\na. Thank you for the suggestion. In the original Methods section, we had specified that ECoG data\nwere band-pass filtered to extract distinct frequency bands, including the “gross” frequency band\ndefined as 0.2-45 Hz.\nb. In this revision, we separated the description of gross frequency band filtering into an individual\nsentence in the Methods section and explicitly stated that the gross frequency band was defined by\napplying a zero-phase Butterworth band-pass filter with cutoff frequencies at 0.2 and 45 Hz.\nspecies relevance.\na. We thank the reviewer for the suggestion. The architecture of large-scale rsfMRI networks appears\nto be largely conserved across species, both in terms of state spatial patterns and transition\nstructures. Specifically, many spatial patterns of rodent rsfMRI networks are homologous to those\nof human networks, including the DMN, SN, and salience network22. Moreover, our results revealed\na temporal hierarchy similar to that reported in human studies, where state transitions were\norganized into sensory and motor-related networks (i.e., SN-FP) and higher-order cognition-related\nnetworks (i.e., DMN-RL)28.\nNevertheless, cross-species differences in rsfMRI network dynamics are evident, particularly in the\nspatial patterns and the presence of certain rsfMRI network states. Specifically, although the spatial\npatterns of rodent rsfMRI networks generally resemble those observed in humans, certain networks\nmay differ22,34, such as medial and lateral frontal cortical connectivity. Furthermore, some networks\nare well characterized in humans, such as the dorsal attention network and language network, but\nremain elusive in rodents, likely due to specialized higher-order cognitive functions as a result of\nthe evolution of mammalian brains22,35. Despite these slight differences, animal models provide\na valuable platform for investigating rapid rsfMRI network reconfiguration and its associated\nfunctions, owing to their largely preserved rsfMRI network architecture22,35 and suitability\nfor integration with neuromodulation techniques (e.g., optogenetics)36.\nb. In this revision, we added the discussions in R3-11a above to highlight cross-species relevance and\ndifferences in the Discussion section.\nReviewer #4\nR4-G Authors have adequately addressed my comments.\nWe sincerely thank the reviewer for the comments and suggestions, which have helped to further\nclarify and strengthen the manuscript.\nThe authors would like to take this opportunity to thank the four reviewers and editor for their\nvaluable critiques and suggestions during the review of our original submission and subsequent\nrevision, which we believe have substantially strengthened this manuscript.\n1. Hipp, J.F., Hawellek, D.J., Corbetta, M., Siegel, M. & Engel, A.K. Large-scale cortical correlation\nstructure of spontaneous oscillatory activity. Nature neuroscience 15, 884-890 (2012).\n2. Krishnan, G.P., Gonzalez, O.C. & Bazhenov, M. Origin of slow spontaneous resting-state neuronal\nfluctuations in brain networks. Proc Natl Acad Sci U S A 115, 6858-6863 (2018).\n3. Vafaii, H., et al. Multimodal measures of spontaneous brain activity reveal both common and\ndivergent patterns of cortical functional organization. Nat Commun 15, 229 (2024).\n4. Anderson, J.S., Carandini, M. & Ferster, D. Orientation tuning of input conductance, excitation, and\ninhibition in cat primary visual cortex. Journal of neurophysiology 84, 909-926 (2000).\n5. Okun, M. & Lampl, I. Instantaneous correlation of excitation and inhibition during ongoing and\nsensory-evoked activities. Nat Neurosci 11, 535-537 (2008).\n6. Poo, C. & Isaacson, J.S. Odor representations in olfactory cortex: \"sparse\" coding, global inhibition,\nand oscillations. Neuron 62, 850-861 (2009).\n7. Kline, A.M., Aponte, D.A., Tsukano, H., Giovannucci, A. & Kato, H.K. Inhibitory gating of\ncoincidence-dependent sensory binding in secondary auditory cortex. Nat Commun 12, 4610 (2021).\n8. Chao, T.-H.H., et al. Neuronal dynamics of the default mode network and anterior insular cortex:\nIntrinsic properties and modulation by salient stimuli. Science Advances 9, eade5732 (2023).\n9. Leong, A.T.L., Wang, X., Wong, E.C., Dong, C.M. & Wu, E.X. Neural activity temporal pattern\ndictates long-range propagation targets. Neuroimage 235, 118032 (2021).\n10. Menon, V. 20 years of the default mode network: A review and synthesis. Neuron 111, 2469-2487\n(2023).\n11. Fox, M.D., et al. The human brain is intrinsically organized into dynamic, anticorrelated functional\nnetworks. Proceedings of the National Academy of Sciences 102, 9673-9678 (2005).\n12. Hinz, R., et al. Bottom-up sensory processing can induce negative BOLD responses and reduce\nfunctional connectivity in nodes of the default mode-like network in rats. Neuroimage 197, 167-176\n(2019).\n13. Xu, N., et al. Functional Connectivity of the Brain Across Rodents and Humans. Front Neurosci 16,\n816331 (2022).\n14. Tsai, P.J., et al. Converging Structural and Functional Evidence for a Rat Salience Network. Biol\nPsychiatry 88, 867-878 (2020).\n15. Menon, V. Large-scale brain networks and psychopathology: a unifying triple network model. Trends\nCogn Sci 15, 483-506 (2011).\n16. Baluch, F. & Itti, L. Mechanisms of top-down attention. Trends Neurosci 34, 210-224 (2011).\n17. Bradley, C., Nydam, A.S., Dux, P.E. & Mattingley, J.B. State-dependent effects of neural stimulation\non brain function and cognition. Nat Rev Neurosci (2022).\n18. Condylis, C., et al. Context-Dependent Sensory Processing across Primary and Secondary\nSomatosensory Cortex. Neuron 106, 515-525 e515 (2020).\n19. Sarter, M., Givens, B. & Bruno, J.P. The cognitive neuroscience of sustained attention: where top-\ndown meets bottom-up. Brain research reviews 35, 146-160 (2001).\n20. Zou, Y., et al. Cell-type-specific optogenetic fMRI on basal forebrain reveals functional network\nbasis of behavioral preference. Neuron 112, 1342-1357. e1346 (2024).\n21. Lozano-Montes, L., et al. Optogenetic Stimulation of Basal Forebrain Parvalbumin Neurons Activates\nthe Default Mode Network and Associated Behaviors. Cell Rep 33, 108359 (2020).\n22. Cole, M.W., Bassett, D.S., Power, J.D., Braver, T.S. & Petersen, S.E. Intrinsic and task-evoked\nnetwork architectures of the human brain. Neuron 83, 238-251 (2014).\n23. Huang, Q., et al. Replay-triggered brain-wide activation in humans. Nat Commun 15, 7185 (2024).\n24. Dimakou, A., Pezzulo, G., Zangrossi, A. & Corbetta, M. The predictive nature of spontaneous brain\nactivity across scales and species. Neuron 113, 1310-1332 (2025).\n25. Liu, M., Liang, Y., Song, C., Knopfel, T. & Zhou, C. Cortex-wide spontaneous activity non-linearly\nsteers propagating sensory-evoked activity in awake mice. Cell Rep 41, 111740 (2022).\n26. Leong, A.T., et al. Long-range projections coordinate distributed brain-wide neural activity with a\nspecific spatiotemporal profile. Proc Natl Acad Sci U S A 113, E8306-E8315 (2016).\n27. Chuang, K.H., Li, Z., Huang, H.H., Khorasani Gerdekoohi, S. & Athwal, D. Hemodynamic transient\nand functional connectivity follow structural connectivity and cell type over the brain hierarchy. Proc\nNatl Acad Sci U S A 120, e2202435120 (2023).\n28. Jung, W.B., Im, G.H., Jiang, H. & Kim, S.G. Early fMRI responses to somatosensory and optogenetic\nstimulation reflect neural information flow. Proc Natl Acad Sci U S A 118(2021).\n29. Yu, X., et al. Sensory and optogenetically driven single-vessel fMRI. Nat Methods 13, 337-340\n(2016).\n30. Chen, B.R., Bouchard, M.B., McCaslin, A.F., Burgess, S.A. & Hillman, E.M. High-speed vascular\ndynamics of the hemodynamic response. Neuroimage 54, 1021-1030 (2011).\n31. Bricault, S., et al. Peripheral contributions to resting state brain dynamics. Nat Commun 15, 10820\n(2024).\n32. Menon, V., et al. Optogenetic stimulation of anterior insular cortex neurons in male rats reveals\ncausal mechanisms underlying suppression of the default mode network by the salience network. Nat\nCommun 14, 866 (2023).\n33. Wang, X., Leong, A.T.L., Chan, R.W., Liu, Y. & Wu, E.X. Thalamic low frequency activity\nfacilitates resting-state cortical interhemispheric MRI functional connectivity. Neuroimage 201,\n115985 (2019).\n34. Stevner, A.B.A., et al. Discovery of key whole-brain transitions and dynamics during human\nwakefulness and non-REM sleep. Nat Commun 10, 1035 (2019).\n35. Salvan, P., et al. Frequency modulation of entorhinal cortex neuronal activity drives distinct\nfrequency-dependent states of brain-wide dynamics. Cell Rep 37, 109954 (2021).\n36. Lu, H., et al. Rat brains also have a default mode network. Proc Natl Acad Sci U S A 109, 3979-3984\n(2012).\n37. Kaefer, K., Stella, F., McNaughton, B.L. & Battaglia, F.P. Replay, the default mode network and the\ncascaded memory systems model. Nat Rev Neurosci 23, 628-640 (2022).\n38. Preston, A.R. & Eichenbaum, H. Interplay of hippocampus and prefrontal cortex in memory. Current\nbiology 23, R764-R773 (2013).\n39. Yeshurun, Y., Nguyen, M. & Hasson, U. The default mode network: where the idiosyncratic self\nmeets the shared social world. Nat Rev Neurosci 22, 181-192 (2021).\n40. Nair, J., et al. Basal forebrain contributes to default mode network regulation. Proc Natl Acad Sci U S\nA 115, 1352-1357 (2018).\n41. Zhu, F., Elnozahy, S., Lawlor, J. & Kuchibhotla, K.V. The cholinergic basal forebrain provides a\nparallel channel for state-dependent sensory signaling to auditory cortex. Nat Neurosci 26, 810-819\n(2023).\n42. Ananth, M.R., Rajebhosale, P., Kim, R., Talmage, D.A. & Role, L.W. Basal forebrain cholinergic\nsignalling: development, connectivity and roles in cognition. Nat Rev Neurosci 24, 233-251 (2023).\n43. Ward Jr, J.H. Hierarchical grouping to optimize an objective function. Journal of the American\nstatistical association 58, 236-244 (1963).\n44. Vidaurre, D., Smith, S.M. & Woolrich, M.W. Brain network dynamics are hierarchically organized\nin time. Proc Natl Acad Sci U S A 114, 12827-12832 (2017).\n45. Cho, J.W., Korchmaros, A., Vogelstein, J.T., Milham, M.P. & Xu, T. Impact of concatenating fMRI\ndata on reliability for functional connectomics. Neuroimage 226, 117549 (2021).\n46. Elliott, M.L., et al. General functional connectivity: Shared features of resting-state and task fMRI\ndrive reliable and heritable individual differences in functional brain networks. Neuroimage 189, 516-\n532 (2019).\n47. Granger, C.W. Investigating causal relations by econometric models and cross-spectral methods.\nEconometrica: journal of the Econometric Society, 424-438 (1969).\n48. Hamilton, J.P., Chen, G., Thomason, M.E., Schwartz, M.E. & Gotlib, I.H. Investigating neural\nprimacy in Major Depressive Disorder: multivariate Granger causality analysis of resting-state fMRI\ntime-series data. Mol Psychiatry 16, 763-772 (2011).\n49. Deshpande, G., LaConte, S., James, G.A., Peltier, S. & Hu, X. Multivariate Granger causality analysis\nof fMRI data. Human brain mapping 30, 1361-1373 (2009).\n50. Raut, R.V., et al. Global waves synchronize the brain’s functional systems with fluctuating arousal.\nScience advances 7, eabf2709 (2021).\n51. Cole, M.W., Ito, T., Bassett, D.S. & Schultz, D.H. Activity flow over resting-state networks shapes\ncognitive task activations. Nat Neurosci 19, 1718-1726 (2016).\n52. Ito, T., Hearne, L.J. & Cole, M.W. A cortical hierarchy of localized and distributed processes revealed\nvia dissociation of task activations, connectivity changes, and intrinsic timescales. Neuroimage 221,\n117141 (2020).\n53. Aguirre, G.K., Zarahn, E. & D'Esposito, M. The variability of human, BOLD hemodynamic\nresponses. Neuroimage 8, 360-369 (1998).\n54. Handwerker, D.A., Ollinger, J.M. & D'Esposito, M. Variation of BOLD hemodynamic responses\nacross subjects and brain regions and their effects on statistical analyses. Neuroimage 21, 1639-1651\n(2004).\n55. Ma, Z. & Zhang, N. Temporal transitions of spontaneous brain activity. Elife 7(2018).\n56. Liu, Y., et al. An open database of resting-state fMRI in awake rats. Neuroimage 220, 117094 (2020).\n57. Yu, Y., et al. Sleep fMRI with simultaneous electrophysiology at 9.4 T in male mice. Nature\nCommunications 14(2023).\n58. Liu, X., et al. Subcortical evidence for a contribution of arousal to fMRI studies of brain activity.\nNat Commun 9, 395 (2018).\n59. Liu, X., Leopold, D.A. & Yang, Y. Single-neuron firing cascades underlie global spontaneous brain\nevents. Proc Natl Acad Sci U S A 118(2021).\n60. Yamagata, T., et al. The hypothalamic link between arousal and sleep homeostasis in mice. Proc Natl\nAcad Sci U S A 118(2021).\n61. Thomas, K., et al. GABAergic neurons in basal forebrain exert frequency-specific modulation on\nauditory cortex and enhance attentional selection of auditory stimuli. Commun Biol 8, 149 (2025).\n62. Gutierrez-Barragan, D., et al. Unique spatiotemporal fMRI dynamics in the awake mouse brain. Curr\nBiol 32, 631-644 e636 (2022).\n63. Huang, Z., Zhang, J., Wu, J., Mashour, G.A. & Hudetz, A.G. Temporal circuit of macroscale dynamic\nbrain activity supports human consciousness. Science advances 6, eaaz0087 (2020).\n64. Wang, X., et al. Functional MRI reveals brain-wide actions of thalamically-initiated oscillatory\nactivities on associative memory consolidation. Nat Commun 14, 2195 (2023).\n65. Yizhar, O., Fenno, L.E., Davidson, T.J., Mogri, M. & Deisseroth, K. Optogenetics in neural systems.\nNeuron 71, 9-34 (2011).\n66. Stujenske, J.M., Spellman, T. & Gordon, J.A. Modeling the Spatiotemporal Dynamics of Light and\nHeat Propagation for In Vivo Optogenetics. Cell Rep 12, 525-534 (2015).\n67. Chan, R.W., et al. Low-frequency hippocampal-cortical activity drives brain-wide resting-state\nfunctional MRI connectivity. Proc Natl Acad Sci U S A 114, E6972-E6981 (2017).\n68. Zhou, I.Y., et al. Brain resting-state functional MRI connectivity: morphological foundation and\nplasticity. Neuroimage 84, 1-10 (2014).\n69. Chuang, K.H. & Nasrallah, F.A. Functional networks and network perturbations in rodents.\nNeuroimage 163, 419-436 (2017).\n70. Liang, Z., Liu, X. & Zhang, N. Dynamic resting state functional connectivity in awake and\nanesthetized rodents. Neuroimage 104, 89-99 (2015).\n71. Arias-Gil, G., Ohl, F.W., Takagaki, K. & Lippert, M.T. Measurement, modeling, and prediction of\ntemperature rise due to optogenetic brain stimulation. Neurophotonics 3, 045007 (2016).\n72. Owen, S.F., Liu, M.H. & Kreitzer, A.C. Thermal constraints on in vivo optogenetic manipulations.\nNat Neurosci 22, 1061-1065 (2019).\n73. Sullivan, G.M. & Feinn, R. Using Effect Size-or Why the P Value Is Not Enough. J Grad Med Educ\n4, 279-282 (2012).\n74. Buzsáki, G. & Mizuseki, K. The log-dynamic brain: how skewed distributions affect network\noperations. Nature Reviews Neuroscience 15, 264-278 (2014).\n75. Vidaurre, D., Llera, A., Smith, S.M. & Woolrich, M.W. Behavioural relevance of spontaneous,\ntransient brain network interactions in fMRI. Neuroimage 229, 117713 (2021).\n76. Song, H., Park, B.Y., Park, H. & Shim, W.M. Cognitive and Neural State Dynamics of Narrative\nComprehension. J Neurosci 41, 8972-8990 (2021).\n77. Song, H., Shim, W.M. & Rosenberg, M.D. Large-scale neural dynamics in a shared low-dimensional\nstate space reflect cognitive and attentional dynamics. Elife 12(2023).\n78. Vidaurre, D., et al. Discovering dynamic brain networks from big data in rest and task. Neuroimage\n180, 646-656 (2018).\n79. Taghia, J., et al. Uncovering hidden brain state dynamics that regulate performance and decision-\nmaking during cognition. Nat Commun 9, 2505 (2018).\n80. Nghiem, T.E., et al. Space wandering in the rodent default mode network. Proc Natl Acad Sci U S A\n121, e2315167121 (2024).\n81. Swick, D., Ashley, V. & Turken, U. Are the neural correlates of stopping and not going identical?\nQuantitative meta-analysis of two response inhibition tasks. Neuroimage 56, 1655-1665 (2011).\n82. Chen, T., Cai, W., Ryali, S., Supekar, K. & Menon, V. Distinct Global Brain Dynamics and\nSpatiotemporal Organization of the Salience Network. PLoS Biol 14, e1002469 (2016).\n83. Lee, S.H. & Dan, Y. Neuromodulation of brain states. Neuron 76, 209-222 (2012).\n84. Allen, P.J., Josephs, O. & Turner, R. A method for removing imaging artifact from continuous EEG\nrecorded during functional MRI. Neuroimage 12, 230-239 (2000).\n85. Jung, W.B., Jiang, H., Lee, S. & Kim, S.G. Dissection of brain-wide resting-state and functional\nsomatosensory circuits by fMRI with optogenetic silencing. Proc Natl Acad Sci U S A 119(2022).\n86. Shepherd, G.M.G. & Yamawaki, N. Untangling the cortico-thalamo-cortical loop: cellular pieces of a\nknotty circuit puzzle. Nat Rev Neurosci 22, 389-406 (2021).\n87. Theyel, B.B., Llano, D.A. & Sherman, S.M. The corticothalamocortical circuit drives higher-order\ncortex in the mouse. Nat Neurosci 13, 84-88 (2010).\n88. Tian, J., et al. Dissection of the long-range projections of specific neurons at the synaptic level in the\nwhole mouse brain. Proc Natl Acad Sci U S A 119, e2202536119 (2022).\n89. Turchi, J., et al. The Basal Forebrain Regulates Global Resting-State fMRI Fluctuations. Neuron 97,\n940-952 e944 (2018).\n90. Parker, P.R.L., Brown, M.A., Smear, M.C. & Niell, C.M. Movement-Related Signals in Sensory\nAreas: Roles in Natural Behavior. Trends Neurosci 43, 581-595 (2020).\n91. Cohen, D., Nakai, T. & Nishimoto, S. Brain networks are decoupled from external stimuli during\ninternal cognition. Neuroimage 256, 119230 (2022).\n92. Anticevic, A., et al. The role of default network deactivation in cognition and disease. Trends Cogn\nSci 16, 584-592 (2012).\n93. Zhang, Q., et al. Brain-wide ongoing activity is responsible for significant cross-trial BOLD\nvariability. Cereb Cortex 32, 5311-5329 (2022)."}]

Decision letter