Long-range projections coordinate distributed brain-wide neural activity with a specific spatiotemporal profile

Significance What makes the brain tick? A simple yet challenging question that has captivated our minds for centuries. This sentiment was fittingly reflected in the launch of The BRAIN Initiative 3 years ago, spurred by the rapid advancement of noninvasive brain imaging and neuronal mapping technologies that have advanced our understanding of neural networks, which are central to brain functions and behavior. Here, we study the patterns of large-scale brain-wide interactions mediated by thalamo-cortical networks through optogenetics and functional MRI. We found that the thalamus can recruit long-range cortical and subcortical networks and initiate their interactions in a spatiotemporally specific manner. This finding provides a fresh impetus to study the mysteries of the brain. One challenge in contemporary neuroscience is to achieve an integrated understanding of the large-scale brain-wide interactions, particularly the spatiotemporal patterns of neural activity that give rise to functions and behavior. At present, little is known about the spatiotemporal properties of long-range neuronal networks. We examined brain-wide neural activity patterns elicited by stimulating ventral posteromedial (VPM) thalamo-cortical excitatory neurons through combined optogenetic stimulation and functional MRI (fMRI). We detected robust optogenetically evoked fMRI activation bilaterally in primary visual, somatosensory, and auditory cortices at low (1 Hz) but not high frequencies (5–40 Hz). Subsequent electrophysiological recordings indicated interactions over long temporal windows across thalamo-cortical, cortico-cortical, and interhemispheric callosal projections at low frequencies. We further observed enhanced visually evoked fMRI activation during and after VPM stimulation in the superior colliculus, indicating that visual processing was subcortically modulated by low-frequency activity originating from VPM. Stimulating posteromedial complex thalamo-cortical excitatory neurons also evoked brain-wide blood-oxygenation-level–dependent activation, although with a distinct spatiotemporal profile. Our results directly demonstrate that low-frequency activity governs large-scale, brain-wide connectivity and interactions through long-range excitatory projections to coordinate the functional integration of remote brain regions. This low-frequency phenomenon contributes to the neural basis of long-range functional connectivity as measured by resting-state fMRI.

One challenge in contemporary neuroscience is to achieve an integrated understanding of the large-scale brain-wide interactions, particularly the spatiotemporal patterns of neural activity that give rise to functions and behavior. At present, little is known about the spatiotemporal properties of long-range neuronal networks. We examined brain-wide neural activity patterns elicited by stimulating ventral posteromedial (VPM) thalamo-cortical excitatory neurons through combined optogenetic stimulation and functional MRI (fMRI). We detected robust optogenetically evoked fMRI activation bilaterally in primary visual, somatosensory, and auditory cortices at low (1 Hz) but not high frequencies . Subsequent electrophysiological recordings indicated interactions over long temporal windows across thalamo-cortical, cortico-cortical, and interhemispheric callosal projections at low frequencies. We further observed enhanced visually evoked fMRI activation during and after VPM stimulation in the superior colliculus, indicating that visual processing was subcortically modulated by low-frequency activity originating from VPM. Stimulating posteromedial complex thalamo-cortical excitatory neurons also evoked brain-wide bloodoxygenation-level-dependent activation, although with a distinct spatiotemporal profile. Our results directly demonstrate that lowfrequency activity governs large-scale, brain-wide connectivity and interactions through long-range excitatory projections to coordinate the functional integration of remote brain regions. This low-frequency phenomenon contributes to the neural basis of long-range functional connectivity as measured by resting-state fMRI. fMRI | optogenetic | brain connectivity | low frequency | thalamus T he brain is a highly complex, interconnected structure with parallel and hierarchical networks distributed within and between neural systems (1,2). This integrative architecture dictates the underlying principles of how brain-wide neural connectivity supports and organizes sensory, behavioral, and cognitive processes (3). Recent advances in structural (2,4) and functional connectivity (5)(6)(7)(8)(9)(10)(11)(12) mapping, as well as neural circuit modulatory tools, such as optogenetics (13,14), permit detailed, high-resolution neural examinations at unprecedented scales. In particular, functional connectivity mapping using resting-state functional MRI (rsfMRI) produces noninvasive visualization of slow and spontaneous hemodynamic fluctuations in humans (5)(6)(7)(8)(9) and animals (10)(11)(12). Coherent low-frequency (<1 Hz) fluctuations across functionally coupled, large-scale networks is an intriguing property, although the exact underlying neural bases and functional significance remain unclear. Previous studies integrating large-scale electrical recordings, voltage-sensitive dye (VSD), and Ca 2+ -imaging techniques (15)(16)(17)(18) support the hypothesis that slow, oscillating neural activity constrains and elicits these hemodynamic fluctuations. Furthermore, low-frequency activity can temporally synchronize remote brain regions for extended periods to enable information processing (19)(20)(21). Although studies demonstrated that rsfMRI connectivity is constrained by well-defined structural networks (3,17), connectivity itself can reorganize during learning-related tasks (22). Given these multifaceted findings, it is now imperative to directly investigate the spatiotemporal properties of large-scale neural interactions. Specifically, how do long-range projection networks facilitate and govern large-scale neural activity?
The cerebral cortex contains long-range axonal projections of pyramidal neurons reciprocally connected to the thalamus by cortico-thalamic and thalamo-cortical projections and within the cortex via cortico-cortical and interhemispheric projections (23). These long-range projections are predominantly glutamatergic and generate most excitatory inputs to their target regions (24). Thalamo-cortical circuits, particularly the somatosensory thalamo-cortical circuits in rodents, are excellent models of longrange projections often studied for their straightforward circuit architecture with monosynaptic excitatory cortical input (25,26). In addition, the thalamus can control cortical states (27)(28)(29) and modulate spontaneous cortical activity based on different cortical states (30,31). These studies suggest that the thalamo-corticothalamic networks are integral to facilitating diverse functional neural integration across the entire brain.
The neuronal architecture within the brain is an intricate web of interconnected excitatory neurons and inhibitory interneurons,

Significance
What makes the brain tick? A simple yet challenging question that has captivated our minds for centuries. This sentiment was fittingly reflected in the launch of The BRAIN Initiative 3 years ago, spurred by the rapid advancement of noninvasive brain imaging and neuronal mapping technologies that have advanced our understanding of neural networks, which are central to brain functions and behavior. Here, we study the patterns of large-scale brain-wide interactions mediated by thalamo-cortical networks through optogenetics and functional MRI. We found that the thalamus can recruit long-range cortical and subcortical networks and initiate their interactions in a spatiotemporally specific manner. This finding provides a fresh impetus to study the mysteries of the brain. which creates highly dynamic yet stable excitation-inhibition networks. Normal operation of this architecture or brain-wide networks requires precise and spatiotemporally structured activity propagation and interaction patterns that support the functional properties of these networks, including downstream activity propagation characteristics within long-range projections. Intuitively, simple perturbations using local electrical microstimulation will not likely permit a reliable probe of the spatiotemporal activity propagation dynamics and their functional output. Previous studies examining the spatiotemporal response properties within thalamocortical circuits were limited by nonspecific electrical stimulation. Electrical stimulation cannot target cell-type-specific neurons, and it excites circuit components in the stimulation region that may not be involved in the particular behavior under physiological conditions. Indeed, electrical stimulation can disrupt both cortico-cortical and cortico-thalamic signal propagation (32) and result in enlarged regions of spatial activation in somatosensory cortex (33). These studies indicate the inadequacy of local electrical stimulation to evoke and characterize downstream activity propagation and interaction within long-range and polysynaptic projection networks. Furthermore, current electrode recording technology limits both the recording site density and spatial coverage. Therefore, thalamo-cortical circuits have been characterized only between thalamic nuclei and directly innervated cortical regions (34,35). This characterization yields a limited understanding into widespread neural interactions beyond thalamo-cortico-thalamic networks. Hence, a disconnection exists at the systems level describing how thalamo-cortical circuits mediate large-scale brain-wide neural activity.
In this study, we sought to characterize large-scale brain-wide neural activity propagation and interaction dynamics and examine their functional roles. We determined the spatiotemporal response properties of long-range excitatory thalamic projections to cortical and subcortical regions. We overcame the limitations of electrical stimulation nonspecificity (34,35) and electrode recording spatial coverage by using a multimodal approach. We used optogenetic functional MRI (fMRI) in combination with electrophysiological recordings to interrogate the vibrissae (whisker) somatosensory thalamo-cortical system. Integration of fMRI and optogenetics permits large-scale view and detection of brainwide neural activity and enables reversible, millisecond precision and cell-type-specific modulation of thalamo-cortical neurons. We selectively targeted thalamo-cortical excitatory neurons in ventral posteromedial thalamus (VPM), a primary lemniscal somatosensory pathway, and stimulated them at different frequencies . The posteromedial complex of the thalamus (POm), a nonlemniscal somatosensory pathway, was also examined. We investigated the spatiotemporal properties of downstream excitatory signal propagation beyond the defined thalamo-cortical circuit, identified the requisite long-range projections, and tested whether VPM stimulation modulates other sensory processing in remote brain regions.

Results
Brain-Wide fMRI Mapping of Downstream Signal Propagation from VPM.
We expressed Ca 2+ /calmodulin-dependent protein kinase II alpha (CaMKIIα)-dependent ChR2(H134R) fused with mCherry in VPM thalamo-cortical excitatory neurons in normal adult Sprague-Dawley rats, as the majority of thalamic nuclei, except for the reticular nucleus, are predominantly populated by CaMKIIα-expressing excitatory neurons (36) (Fig. 1A and Fig. S1A). We confirmed specific expression in VPM excitatory neurons through colocalization of mCherry with CaMKIIα staining. We performed optogenetic stimulation in lightly anesthetized rats (1% isoflurane). Blue light pulses at five frequencies (1 Hz: 10% pulse width duty cycle; 5, 10, 20, and 40 Hz: 30% duty cycle; light intensity: 40 mW/mm 2 ) were  (1,5,10,20, and 40 Hz) were presented in a pseudorandomized order. All frequencies were presented at 30% duty cycle except for 1 Hz that was presented at 10% duty cycle. delivered to VPM neurons in an interleaving block paradigm ( Fig. 1  B and C). These frequencies cover natural whisking frequencies (5 Hz and 10 Hz), below (1 Hz) and above (40 Hz). We chose a reduced duty cycle of 10% for l Hz of stimulation to avoid excessively long stimulation pulse widths that may not be physiological. We made the exposed optical fiber cannula (used to deliver blue light pulses) opaque using heat-shrinkable sleeves to prevent light leakage that may cause undesired visual stimulation to animals.
We detected robust large-scale blood-oxygenation-leveldependent (BOLD) fMRI activation within and beyond the somatosensory thalamo-cortical circuit at 1 Hz of ipsilateral VPM stimulation ( Fig. 2 and Fig. S2). We identified robust and diffuse BOLD responses in a contralateral primary somatosensory barrel field (S1BF) and upper lip area (S1ULp), secondary somatosensory cortex (S2), bilateral primary visual (V1) and auditory (A1) cortices, and cingulate cortex (Cg), in addition to ipsilateral SIBF, S1ULp, and S2 (see Fig. 2C and Fig. S2C for the BOLD signal profiles). Despite the long interstimulus duration and low duty cycle, we reliably detected BOLD activation beyond thalamo-cortical networks, which revealed recruitment of brain-wide and population-wide neural activity by low-frequency stimulation. We also performed 1 Hz of stimulation at a lower light intensity (20 mW/mm 2 ; Fig. S3A) and observed no apparent change in the overall spatial extent of broad cortical BOLD responses although BOLD signal amplitudes decreased as expected. At 5 and 10 Hz, BOLD activation was no longer present in contralateral somatosensory cortices and bilateral V1 and A1 ( Fig.  2 and Fig. S2). Instead, robust and diffuse BOLD activation was detected in ipsilateral primary motor cortex (M1), in addition to ipsilateral S1BF, S2, and S1ULp. At 20 Hz and 40 Hz, we detected only BOLD activation in ipsilateral S1BF, S2, and S1ULp. We observed no BOLD activation in naive animals under the same light stimulation paradigm. The absolute BOLD signal amplitude within the ipsilateral S1BF at 1 Hz was generally lower than signal amplitudes at 5 Hz and 10 Hz mainly due to the reduced-pulsewidth duty cycle of 10% used at 1 Hz. VPM BOLD activation was not apparent at 1 Hz because the lower duty cycle elicited a weaker BOLD signal in VPM.
The above fMRI visualization of large-scale thalamo-cortical downstream signal propagation demonstrates that evoked spatiotemporal dynamics are not restricted to monosynaptic thalamocortical projections from VPM to somatosensory cortices ( Fig. 2 and Figs. S1B and S2). We found frequency-specific response properties that transition activation from somatosensory-visual and somatosensory-auditory regions at low frequency (1 Hz) to sensorimotor areas at whisking frequency (5 Hz and 10 Hz). Together, these fMRI results indicate that excitatory signals evoked by optogenetic stimulation of VPM propagate to remote cortical regions polysynaptically and robustly at low frequencies.
Extracellular Recordings Underlying Brain-Wide BOLD Activations.
Guided by our fMRI results, we performed local field potential (LFP) recordings using the same stimulation paradigm as shown in Fig. 1C at five sites: ipsilateral VPM, bilateral S1, and bilateral V1 (Fig. 3A). Averaged 30-s traces of LFP recordings largely concurred with the magnitude and temporal profiles of BOLD activation at ipsilateral S1BF and bilateral V1 (Fig. 3B vs. Fig. 2C). The temporal evolution of BOLD signals was generally slower due to the inherent nature of the hemodynamic response to neural activity (37). Similar levels of LFP responses were elicited at all frequencies (1-40 Hz) for the first optogenetic pulse at each recording site. LFP response levels to optogenetic pulses were largely constant at 1 Hz. However, they decreased quickly during consecutive pulses at 5 Hz and 10 Hz and sharply decayed after the first pulse at 20 Hz and 40 Hz. LFP recordings performed in naive animals showed no evoked LFPs at these five sites upon simple VPM light stimulation (Fig. S4). In addition, laser light pulses can generate high photocurrents in neurons (38), which may produce unphysiological states. Therefore, we examined LFP recordings during 1 Hz of stimulation using different light intensities with the default 40 mW/mm 2 vs. 20 mW/mm 2 and 15 mW/mm 2 (Fig. S3B). We observed no apparent changes in LFP characteristics. LGN: lateral geniculate nucleus; V2: secondary visual cortex; Cg: cingulate cortex. (B) Averaged BOLD activation maps at 1, 5, 10, 20, and 40 Hz of optogenetic stimulation (n = 6; P < 0.001; asterisk: stimulation site). Robust and diffuse positive BOLD activations were observed in bilateral S1BF, V1/V2, A1, S2, and Cg cortices upon 1 Hz of stimulation. When stimulation frequency was increased to 5 and 10 Hz, robust and diffuse positive BOLD activations were observed in ipsilateral S1BF, M1, S1ULp, and S2. Further increase in stimulation frequency to 20 and 40 Hz activated only the ipsilateral S1BF, S1ULp, and S2. (C) BOLD signal profiles extracted from ROIs defined in A (error bar indicates ± SEM).
These LFP recordings showed that this frequency-dependent adaptation occurs first in VPM. We measured the LFP peak amplitude decay constant, τ, for 5 and 10 Hz stimulation ( Fig. 3 C and D). The measurements indicated further adaptation in ipsilateral S1, as τ S1 was larger than τ VPM (2.10 ± 0.15 s vs. 1.50 ± 0.17 s, P < 0.01 for 5 Hz; 1.05 ± 0.13 s vs. 0.73 ± 0.11 s, P < 0.05 for 10 Hz). Hence, simple failure of VPM excitatory neurons due to sustained membrane depolarization at high stimulation frequencies does not underlie this adaptation phenomenon. Both cortico-thalamic and intracortical interactions can contribute to these dynamic properties.
Multiunit activity (MUA) recordings in VPM using an optrode showed that spikes were optogenetically evoked in thalamo-cortical neurons at all stimulation frequencies (Fig. 4A). Spike frequency analysis demonstrated that optogenetic pulses successfully evoked action potentials in a stimulus-locked manner, even at high frequencies including 20 Hz and 40 Hz, the range at which noncell-type-specific electrical stimulation fails (39). We observed an onset latency of 3.3 ± 1.4 ms, consistent with VPM thalamo-cortical neuron characteristics (40). Surprisingly, MUA recordings revealed different spiking patterns that transitioned from thalamic burst activity to tonic activity with increasing stimulation frequency. As long stimulation pulse durations can confound evoked thalamic burst activity, we repeated these recordings using a 10-ms stimulation pulse and observed similar spiking characteristics (Fig.  4B). These findings indicate that the occurrence of evoked burst activity in VPM is determined through interstimulus intervals, not stimulus duration.

Long-Range Excitatory Projections Support Signal Propagation at
Low Frequencies and Interactions with Long Temporal Windows. To identify downstream signal propagation pathways at 1 Hz, we analyzed the latency of evoked LFPs by using the first LFP peak location to represent the initial neuronal population activity (Fig. 5 A-C). Negative LFP peaks were observed first and used for ipsilateral VPM and ipsilateral and contralateral S1, whereas positive ones were observed for ipsilateral and contralateral V1. Analysis revealed a response latency of 3.4 ± 0.3 ms for VPM. As shown in Fig. 5D, a delay of 9.6 ± 0.6 ms occurred between VPM and ipsilateral S1, consistent with the reported 8.5-ms latency value for monosynaptic VPM to S1 projection (41). Signal then propagated from ipsilateral S1 to ipsilateral V1 through cortico-cortical circuit projections in 21.5 ± 1.6 ms, in line with the reported value of 25 ms (17), and to contralateral S1 through interhemispheric monosynaptic callosal projection in 6.7 ± 0.5 ms, consistent with the 6.0 ms estimated from reported axonal conduction velocity (42). Finally, the signal propagated in 6.2 ± 0.9 ms from ipsilateral V1 to contralateral V1 through interhemispheric callosal projection. However, we cannot preclude that propagation from contralateral S1 to contralateral V1 occurs through cortico-cortical circuit projections in 21.0 ± 1.8 ms because it closely matches the ipsilateral S1 to ipsilateral V1 delay. Taken together, these analyses indicate direct participation of long-range excitatory projections as demonstrated through their initial population activity propagation delays.
As shown in Fig. 5A and in contrast to Fig. 5E, the averaged LFP evoked at 1 Hz typically persisted for 200-300 ms, which included strong secondary peaks long after the primary peaks, Fig. 3. LFP recordings in sensory cortices and VPM confirm neuronal activity underlying BOLD activation and indicate adaptation at the thalamic and cortical levels. (A) Illustration of the location of recording electrodes in bilateral S1BF, bilateral V1, and optrode in VPM. (B) Averaged LFP from a representative animal during VPM optogenetic stimulation (eight blocks; 1 Hz, 10% duty cycle; 5, 10, 20, and 40 Hz, 30% duty cycle). Response levels of evoked LFPs remained fairly constant at 1 Hz. However, evoked LFPs decreased quickly during the initial consecutive pulses at 5-40 Hz. These LFP recordings were consistent with the observed BOLD responses. (C) Adaptation characterization at ipsilateral VPM during 5 Hz of optogenetic stimulation. Adaptation was quantified by normalization of the amplitudes of positive peaks at each of the evoked potentials to the amplitude of the first positive peak, which was then fitted with a monoexponential decay curve to quantify its decay constant, τ. (D) There was additional adaptation in ipsilateral S1 at 5 Hz (τ VPM = 2.10 ± 0.15 s vs. τ S1 = 1.50 ± 0.17 s; error bar indicates SD; paired t test with post hoc Holm-Bonferroni correction; **P < 0.01) and 10 Hz (τ VPM = 1.00 ± 0.13 s vs. τ S1 = 0.73 ± 0.11 s; *P < 0.05). This finding indicates that the observed adaptation involved interactions within the thalamo-cortico-thalamic network and was not solely from sustained VPM membrane depolarization.
particularly in contralateral and ipsilateral V1 and contralateral S1. This finding indicates that a long interstimulus duration (>200 ms) is required to allow and enable recurrent, long-lasting interactions between and within long-range cortico-cortical and thalamo-cortico-thalamic networks. We also performed LFP recordings at 1 Hz using a shorter pulse duration (10 ms instead of 100 ms). We ob-served no significant changes in overall latencies and LFP temporal characteristics (Fig. 5F).

Visual Subcortical BOLD Activation Enhanced by Low-Frequency VPM
Stimulation. To test whether low-frequency activity over longrange projections influences cortical and subcortical processing of visual stimuli, we used a dual stimulation paradigm with a continuous 1-Hz VPM optogenetic stimulation and binocular visual stimulation (Fig. 6A). Visual stimulation was applied before (baseline), during (OG-On), and after (OG-Off) continuous VPM stimulation.
Baseline visual stimulation evoked positive BOLD activation along the visual pathway, including the lateral geniculate nucleus (LGN), superior colliculus (SC), and visual cortices (V1 and V2) (Fig. 6B). This finding was expected and consistent with the previous studies (43)(44)(45). In the presence of continuous optogenetic stimulation, visually evoked BOLD responses significantly increased in the SC (Fig. 6 C and D; ipsilateral SC, P < 0.01; contralateral SC, P < 0.01; one-way ANOVA). Moreover, visual-stimulation-evoked BOLD responses remained elevated in SC after optogenetic stimulation during the ∼2.5-h experimental window (Fig. 6 C and D). These findings indicate that excitability in the visual pathway is modulated by low-frequency activity over long-range projections originating from the primary somatosensory thalamic nucleus, VPM.
fMRI Mapping of Downstream Signal Propagation from POm. We also used optogenetic fMRI to examine brain-wide BOLD activity evoked by excitatory neuronal stimulation in the nonlemniscal somatosensory thalamus, POm (Fig. 7). We observed activation of thalamo-cortical and thalamo-tectal projections upon 1 Hz of stimulation, which evoked robust BOLD activation in ipsilateral S1BF, bilateral SC, and LGN. At 5 Hz, BOLD activation was detected in ipsilateral S1BF, bilateral V1, LGN, and SC. At 10 Hz, BOLD activation was observed only in ipsilateral S1BF (Fig. 7 B and C). Neither ipsilateral S1BF nor SC was activated at 40 Hz of stimulation. These findings again demonstrate the brain-wide and population-wide neural activity evoked by low-frequency stimulation of POm. Furthermore, these spatiotemporal response characteristics are distinct from those observed by VPM stimulation, highlighting the differential roles of VPM and POm in coordinating brain-wide neural activity. Histological examination of POm animals confirmed that local viral expression was restricted to the POm nucleus with no observable expression in VPM (Fig. S5).

Discussion
Long-range projections are essential for neural communication between remote brain regions. The brain has an intrinsic capability to recruit various projections within interconnected networks to serve various functional demands without requiring individual long-range projections to perform a single function. Despite the wealth of positive evidence, it remains a challenge to establish how long-range projection recruitment and modulation occur across the broad range of neural targets. Furthermore, although long-range projections are prevalent in all mammalian brains, most studies focus on large-scale cortico-cortical interactions with minimal emphasis on modulation through subcortical interactions. Here, we examined both issues using optogenetics and fMRI to monitor brain-wide neural responses during thalamic stimulation. Our fMRI and subsequent extracellular recording results show that long-range cortico-cortical and interhemispheric callosal projections are recruited polysynaptically during low-frequency stimulation of VPM excitatory neurons. Similarly, POm stimulation leads to brain-wide neural activity but with a distinct spatiotemporal profile. Our results indicate that this low-frequency phenomenon likely is causally related to reciprocal interactions over thalamo-cortico-thalamic, cortico-cortical, and interhemispheric callosal projections to facilitate activity propagation with extended time windows. Furthermore, visual subcortical SC responses to visual stimulation are enhanced during and after sustained activation of long-range projections via VPM stimulation, which suggests that long-range projections and interactions shape multisensory integration processes in the brain. Together, these results indicate that low-frequency activity mediates brain-wide interactions through long-range excitatory projections. Such flexible characteristics of low-frequency, brain-wide neural activity can contribute to the neural basis of long-range functional organization and connectivity as measured by rsfMRI (5-12), extracellular recording (15,16), VSD (17), and Ca 2+ imaging (18).
Surprisingly, we demonstrated that low-frequency stimulation of excitatory neurons in VPM evokes BOLD responses in regions far beyond somatosensory cortices (Fig. 2). Reciprocal projections between the thalamus and recurrent cortical circuits distinguish the mammalian sensory cortex. The laminar and columnar architecture of the sensory cortex offers the most intuitive structural design with these specific functional characteristics. Sensory cortical pyramidal neurons receive thalamic input at specific cortical layers and distribute excitatory signals within and outside the cortex, whereas inhibitory interneurons are interspersed throughout cortical layers to provide the requisite excitation-inhibition balance. Despite exclusive stimulation of excitatory thalamo-cortical neurons in VPM, downstream signal propagation recruits both S1 excitatory pyramidal neurons and inhibitory interneurons at thalamocortical synapses (46). Most pyramidal neurons have axon collaterals that project back to and synapse in superficial layers, whereas others distribute excitation horizontally. These collaterals form a strong, recurrent excitatory network both within the local microcircuit and among long-range thalamo-cortico-thalamic and cortico-cortical networks (47). The strong input amplification elicited through this positive feedback loop is tightly regulated by various GABAergic interneurons (48) interposed within this pyramidal network to target different neuronal subdomains (23). In our study, evoked LFPs by 1 Hz of VPM stimulation persisted over long periods and exhibited multiple peaks in the ipsilateral VPM, bilateral S1, and V1 (Fig. 5 A  and E). This observation suggests that feedforward and feedback mechanisms governing cortico-cortical and thalamo-cortico-thalamic networks play a causal role to modulate these LFP responses and directly contribute to robust brain-wide BOLD activations detected during 1 Hz of stimulation.
The question remains how and why cortical long-range projections are recruited. We found that low-frequency VPM stimulation elicited thalamic burst activity, whereas high-frequency stimulation elicited tonic activity (Fig. 4A). Thalamo-cortical excitatory neurons have intrinsic membrane potential properties that allow them to fire in two modes, namely burst and tonic/ single spike (49), which require cortico-thalamic activity (50, 51). Optogenetic stimulation of VPM thalamo-cortical neurons elicited action potentials that excited ipsilateral S1 through thalamo-cortical projections. Subsequently, lowfrequency activity propagated cortico-cortically from ipsilateral S1 to ipsilateral V1, and interhemispheric callosal projections between ipsilateral and contralateral S1 were recruited. Signal then was propagated from ipsilateral V1 to contralateral V1 via interhemispheric callosal projections. (E) Averaged evoked LFPs during 5 Hz of VPM optogenetic stimulation (n = 3; 30% duty cycle; error bar indicates ± SD). Comparision with 1 Hz reveals the absence of secondary peaks in ipsilateral and contralateral V1. This indicates that long-lasting and recurring interactions between and within cortico-cortical and thalamo-corticothalamic networks require long interstimulus durations. (F) Evoked potentials from a representative animal during 1 Hz of VPM optogenetic stimulation using a short stimulation pulse (1% duty cycle) shows no significant differences in overall latencies and LFP temporal characteristics.
Thalamic burst activity, originally identified in the sleep/anesthetized state, is a proposed mechanism to decouple thalamocortical activity from the peripheral world (52). However, recent evidence suggests that thalamic bursting plays a fundamental role in sensory processing (53,54). Studies of visual (55) and somatosensory (30,53) pathways demonstrate the presence of thalamic bursting in awake and attentive states, implying that bursting activity modulates sensory perception. Intermittent activation of multiple long-range circuits to integrate major sensory cortices suggests that multisensory integration is essential to shape sensory perceptions that could not otherwise be achieved with only one sense (56,57). This premise raises an intriguing proposition that excitatory long-range projections and their dynamic interactions are recruited during low-frequency thalamic stimulation through burst activity propagation from VPM to S1 and, possibly, to other sensory cortices. In this study, we observed that thalamic burst activity was successfully elicited during VPM stimulation and sustained only when the interstimulus interval was more than 200 ms (Fig. 4A), when brain-wide BOLD activation was observed. Such temporal characteristics indicate that the brain-wide integration of spiking information from multiple thalamic neurons simultaneously, possibly through bursting activity, requires extended time windows. These interactions presumably synchronize remote brain regions that require large conduction delays (19)(20)(21).
Synaptic depression by anesthesia at the stimulation site may diminish polysynaptic activity propagation at high frequencies (39). However, our LFP recordings after cessation of light anesthesia show that evoked V1 responses did not exhibit apparent changes, despite increases in fidelity of signal transmission from VPM to S1 (Fig. S6). A previous study in awake rats showed that desynchronized activity occurred during whisking, as indicated by decreased slow oscillations (∼1 Hz) in somatosensory cortex, and that these effects could be mimicked using high-frequency optogenetic stimulation of thalamo-cortical excitatory neurons (28). These findings support our contention that high-frequency thalamic stimulation may not activate long-range circuits beyond its projected area, even during the awake state, when high-frequency signals (>20 Hz) are transmitted reliably to the cortex. Therefore, the recruitment of long-range projections by low frequency is not an effect of anesthesia.
In this study, we show that constant recruitment of dynamic network interactions through continuous low-frequency VPM stimulation increased SC BOLD activation to external visual stimulus and that such SC enhancement persisted after cessation of VPM optogenetic stimulation (Fig. 6). This observation indicates that sensory gain, a fundamental component of sensory information processing (58,59), can be modulated in the SC during and after continual activation of long-range circuits. Low-frequency stimulation coincided with VPM burst activity, which is modulated by corticothalamic input (50,51). Furthermore, thalamic bursting activity can increase due to enhancement of sensory feedback gain from cortical to thalamic regions (50). Hence, the enhanced SC activation observed here could result from burst information relayed to sensory cortices from VPM and facilitated through direct cortico-tectal projections from V1 and/or S1 (24,60). Neural plasticity mechanisms are also likely recruited. These mechanisms can readjust the neuronal components stimulating homeostatic regulation of neural circuit functions, such as neuronal or network excitability changes. Here, we define neural plasticity as the ability of the nervous system to adopt a new functional or structural state in response to extrinsic and intrinsic factors (61). Plasticity occurs at different neurodevelopmental stages spontaneously (62) and as driven by sensory experiences (63). Such sequential combination of spontaneously LGN: lateral geniculate nucleus; V1/V2: primary/secondary visual cortex. Averaged BOLD activation maps before, during, and after 1 Hz of optogenetic stimulation (n = 6; P < 0.001; asterisk: stimulation site) (Bottom). (C) BOLD signal profiles extracted from ROIs defined in B (error bar indicates ± SEM). BOLD activations in ipsilateral and contralateral SC were enhanced during optogenetic stimulation, which remained elevated after cessation of optogenetic stimulation. These findings indicate that excitability in the visual pathway is modulated by low-frequency activity over long-range projections originating from VPM. (D) Averaged β-values extracted from ROIs defined in B (n = 6; error bar indicates ± SEM, one-way ANOVA with post hoc Bonferroni correction; **P < 0.01; n.s, not significant). generated and experience-dependent neural activity endows the brain with the ability to accommodate changing, dynamic inputs throughout the life span. Studies demonstrate that both spontaneous and experience-dependent neural activity predominates at low frequencies (64)(65)(66), which further corroborates our findings that lowfrequency activity is intimately coupled with long-range projections and their dynamic interactions.
At present, functional connectivity is defined as the temporal coherence of neural activity patterns among spatially segregated brain regions at low and ultra-low frequencies, which has been most widely measured by rsfMRI (5)(6)(7)(8)(9)(10)(11)(12). However, interpretation of functional connectivity studies remains challenging because functional coupling changes dynamically. This dynamic property suggests that it is constrained by, but not fully dictated by, anatomical connectivity. Furthermore, the neural mechanisms underlying long-range, functional coupling remain unknown. Few studies have investigated this question directly (21,67). In this study, we also performed rsfMRI measurements (Fig. S7). Our preliminary rsfMRI results show an increase in long-range interhemispheric connectivity in S1BF, S2, V1, and A1 during and after low-frequency VPM stimulation. Our findings suggest that functional connectivity can be directly altered by low-frequency stimulation. Thus, brain-wide functional connectivity is intimately linked to low-frequency neural activity.
In this study, we also examined the nonlemniscal somatosensory thalamus, POm (Fig. 7). We observed brain-wide neural activity upon POm stimulation at low frequencies (1 and 5 Hz). Interestingly, its spatiotemporal response properties are different from those observed through VPM stimulation. The thalamus is a relay center for cortical sensory information. Axonal tracing (68), electrophysiology (41,69), and behavioral (70) studies show parallel, but functionally distinct, thalamo-cortical projections, including lemniscal VPM vs. paralemniscal POm. Previous studies also characterized the differences between VPM and POm thalamocortical projections in S1BF, in which VPM thalamo-cortical projections define an essential cortical activity pattern, whereas POm thalamo-cortical projections affect the effectiveness of transmission (69,71). Here, our fMRI results reveal that POm stimulation yields highly different spatiotemporal response characteristics from those by VPM stimulation. These findings suggest that POm is functionally distinct from VPM and that the paralemniscal pathway also plays an essential role in multisensory processing.
Optogenetic VPM or POm stimulation might antidromically activate somatosensory cortex. However, the experimental evidence from our fMRI and electrophysiological recordings strongly indicates that the observed spatiotemporally distinct brain-wide responses are likely not caused by antidromic S1 activation. Specifically, MUA recordings show the latency between optogenetic pulse onset and initial spikes as 3.3 ± 1.4 ms and 4.3 ± 2.3 ms for VPM and POm stimulation, respectively ( Fig. 4A and Fig. S5B). These latencies are characteristic of orthodromic activation of VPM or POm thalamo-cortical excitatory neurons (40,72), whereas antidromic activation exhibits a shorter latency with less variation (73). The VPM to ipsilateral S1 propagation time of 9.6 ± 0.6 ms (Fig. 5D) also corresponds well to published values for orthodromic propagation (41). These timing measurements support that our VPM or POm stimulation activates S1 orthodromically rather than antidromically, at least in a predominant manner. In addition, antidromic activation is unlikely because the anterograde and monosynaptic nature of our used viral transduction expresses the optogenetic construct only in VPM or POm and its ascending cortical targets (i.e., AAV5-CaMKIIα-ChR2) (Figs. S1 and S5), without further transfecting descending corticothalamic projections from somatosensory cortex (34,74).
In summary, our results directly demonstrate that low-frequency activity supports large-scale brain-wide interactions over long-range excitatory projections, such that low-frequency activity plays a pivotal role in functional integration of remote brain regions. Such spatiotemporal characteristics can contribute to the neural basis of brain functional connectivity measured by rsfMRI. In addition, this study presents the combined use of large-view fMRI, optogenetic neuromodulation, and external stimuli or tasks as a valuable means to dissect brain long-range interactions and their functions over a large scale.  LGN: lateral geniculate nucleus; SC: superior colliculus; S1ULp: primary somatosensory upper lip area; S2: secondary somatosensory cortex; M1: motor cortex; V1: primary visual cortex, V2: secondary visual cortex; A1: primary auditory cortex. (B) Averaged BOLD activation maps at 1, 5, 10, and 40 Hz of optogenetic stimulation (n = 6; P < 0.001; asterisk: stimulation site). A 1-Hz stimulation showed positive BOLD activation at visual subcortical regions, LGN and SC. Robust positive BOLD activation was observed in ipsilateral S1BF during 5 Hz and 10 Hz of stimulation, which then expanded to visual subcortical activations at 5 Hz. Stimulation at 40 Hz did not show any significant BOLD activation. These findings indicate a specific spatiotemporal response profile different from that by VPM stimulation, highlighting the differential roles of POm and VPM in coordinating brain-wide neural activity. (C) BOLD signal profiles extracted from ROIs defined in A (error bar indicates ± SEM). Viral titer in particles per milliliter was 4 × 10 12 for AAV5-CaMKIIα-ChR2 (H134R)-mCherry. Maps are available online from www.stanford.edu/group/ dlab/optogenetics. Stereotactic Surgery for Viral Injection. After the craniotomy was made, viral injections were performed at two depths in VPM [−3.6 mm posterior to Bregma, +3.0 mm medial-lateral (ML) right hemisphere, −5.8 mm and −6.2 mm from surface of dura] and POm (−3.6 mm posterior to Bregma, +2.1 mm ML right hemisphere, −5.0 mm and −5.4 mm from surface of dura) with subsequent verification by high-resolution anatomic MRI. Viral constructs were delivered at a 1.5-μL volume at each depth through a 5-μL syringe and 33-gauge beveled needle injected at a rate of 150 nL/min. Animals recovered for 4 wk before fMRI experiments were performed.
Animal Preparation for MRI Experiment. Surgery was performed under 2% isoflurane to implant an optical fiber cannula with a 450-μm diameter at the injection site with high-resolution anatomic MRI verification. The fiber-tip surface was beveled to facilitate the fiber insertion and minimize injury to brain tissue. Dental cement fixed the fiber on the skull.

MRI-Synchronized Optogenetic and Visual Stimulation.
For optogenetic stimulation, blue light was delivered to VPM or POm nucleus using a 473-nm diodepumped solid-state laser with calibrated power of 8 mW (40 mW/mm 2 ) at the fiber tip. Optical fiber cannulas were made opaque using heat-shrinkable sleeves to prevent light leakage that may cause undesired visual stimulation. Both VPM and POm are large thalamic nuclei (Figs. S1 and S5). The optical fiber tip was typically situated in the center of the VPM or POm nucleus through stereotaxic implantation. The spatial spread from the fiber tip for 473 nm of 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 direction (74), confining the optogenetic stimulation within the VPM or POm nucleus. For both VPM and POm fMRI experiments, five frequencies were used (1 Hz with 10% duty cycle; 5, 10, 20, and 40 Hz with 30% duty cycle). Here, we chose a reduced duty cycle of 10% for 1 Hz of stimulation (in contrast to 30% for 5-40 Hz) to avoid the use of very long stimulation pulse widths (300 ms). Long ChR2 optogenetic pulse width can cause neuronal output silencing and adversely affect the neural activity efficacy (75). VPM or POm was stimulated with a block design paradigm that consisted of 20-s-on and 60-s-off periods. A total of nine fMRI sessions were performed on each animal. Each session included 10 blocks, with each frequency occupying two blocks. Different frequencies were presented in a pseudorandom manner. For dual optogenetic/ visual stimulation, additional blue light stimulation was presented via a separate optical patch cable at 20 mm in front of both eyes, with power calibrated at 0.5 mW. Visual stimulation (1 Hz and 10% duty cycle) was presented in a block design paradigm of 20-s-on and 40-s-off periods. Each fMRI session included four blocks of visual stimulation. Fifteen fMRI sessions were performed on each animal, 5 before (baseline), 5 during (OG-On), and 5 after (OG-Off) optogenetic stimulation. For rsfMRI experiments, 15 fMRI sessions were performed on each animal, 5 baseline, 5 OG-On, and 5 OG-Off.
MRI Acquisitions. All MRI experiments were carried out on a 7T MRI scanner. Scout and anatomical images were first acquired for accurate positioning and reference. fMRI and rsfMRI measurements were made using a single-shot gradient-echo echo-planar imaging (GE-EPI) sequence with 32 mm 2 × 32 mm 2 field of view, 64 × 64 matrix, flip angle 56°, and echo time/repetition time (TE/TR) = 20/1,000 ms. During MRI, anesthesia was maintained with 1.0% isoflurane. Animal heart rate, respiration rate, oxygen saturation, and rectal temperature were continuously monitored (10,43,76). Note that the fiber outside brain was made opaque using heat-shrinkable sleeves to avoid artificial visual stimulation.
fMRI and rsfMRI Data Analysis. fMRI analysis was performed using the procedure established in our previous studies (43,76). In brief, all GE-EPI images were first coregistered for each animal using SPM8 and registered to a custom-made brain template. For individual animals, data from fMRI sessions were averaged, smoothed, and high-pass-filtered. Then, a general linear was applied to calculate the response coefficient (β) maps for each stimulus. Student's t test was performed to identify activated voxels using the threshold P < 0.001 and cluster size ≥4. BOLD temporal profiles were extracted from regions of interest (ROIs) delineated from the rat brain atlas. rsfMRI analysis was performed according to a procedure reported in our previous study (10).
LFP and MUA Recordings and Analyses. Simultaneous LFP recordings at five brain regions were performed on animals using the same anesthesia protocol (1.0% isoflurane) as in MRI experiments. Craniotomies (1-mm diameter) were made over bilateral S1 (+0.1 mm anterior to Bregma, ±5.0 mm ML), bilateral V1 (−6.0 mm posterior to Bregma, ±4.0 mm ML), and a 5-mm-diameter craniotomy over the stimulation site VPM (−3.6 mm posterior to Bregma, +3.0 mm ML right hemisphere). Single tungsten microelectrodes (1 MΩ and 10-μm tip diameter) were placed in bilateral S1 and V1 approximately between layers IV and V. For VPM, a homemade optrode was used. It consisted of a single tungsten microelectrode (2 MΩ and 10-μm tip diameter) coupled to the fiber using protocols as described previously (77). LFP and MUA recordings were acquired using a 32-channel OpenEphys system with a 30-kHz sampling frequency. The optogenetic stimulation paradigm was the same as used in fMRI experiments. Again, these optical fiber cannulas were made opaque using heat-shrinkable sleeves to prevent light leakage that causes inadvertent direct visual stimulation.
Additional optogenetic stimulation paradigms were used to compare the effects of short-and long-light pulse durations (10/100, 10/60, and 10/30 ms for 1, 5, and 10 Hz, respectively) and different light intensities (40,20, and 15 mW/mm 2 ) on the evoked LFP responses. Additional anesthesia protocols were used to determine the effects of different anesthesia levels (0.8, 0.5, and 0.0% isoflurane) on evoked LFP responses. Recordings under 0.0% isoflurane were performed at 5 min and 15 min after cessation of isoflurane.
For LFP analysis, raw data were band-pass-filtered (1-300 Hz) and notchfiltered for 50 Hz using Matlab. Individual LFP and group-averaged LFP traces were generated. The LFP latency was measured from the onset of light pulse to the first peak location. MUA analysis was performed using Offline Sorter software (Plexon Inc.) for spike sorting and Neuroexplorer software (Nex Technologies) and Matlab for calculating spike frequencies and generating averaged traces. For spike identification, spikes were first sorted using an automated shape-based spike-sorting method, in which a negative and positive threshold crossing at twice the noise SD was imposed, and then this sorting was verified by visual inspection before spike frequency analysis. To quantify the optogenetically evoked spike frequency during the 20-s optogenetic stimulation block, spikes were identified during the optogenetic stimulation pulse period and verified to be stimulus-locked. That is, only the first spike within a predefined 2-to 5-ms window after the onset of each optogenetic stimulation pulse was counted as a spiking event, and the subsequent spikes within the stimulation pulse period were discarded. To quantify the spontaneous spike frequency before and after the 20-s optogenetic stimulation block, the sorted spikes were counted by imposing an interspike interval of 100 ms that was determined based on the temporal characteristics of the thalamic burst activity predominantly observed in the absence of stimuli. The above procedure ensured that both thalamic burst (typically, consecutive multiple spikes within the 20-to 80-ms window) and tonic activity (single spikes lasting less than 5 ms) were characterized as one single spiking event to measure optogenetically evoked spike frequency. The procedure was sensitive to detecting the failure of optogenetic pulses in evoking spikes.
Histology, Immunohistochemistry, and Confocal Imaging. Upon completion of all fMRI and electrophysiological recording experiments, a selected number of animals (n = 3 for VPM; n = 2 for POm) were transcardially perfused, and their brains sectioned and prepared. Consecutive sections (120 μm apart) were prepared according to a previously published procedure (78). Double or triple immunofluorescence was assessed using a confocal scanning laser microscope.
Data Availability. The data that support the findings of this study are available from the corresponding author on request.