Distinct local and brain-wide networks are activated by layer-specic optogenetic stimulations of motor cortex

Primary motor cortex consists of a stack of interconnected but distinct layers, and plays a prominent role in motor control through large-scale networks. However, differential effects of M1 layer-specic functional pathways remain elusive, especially at the macroscopic and mesoscopic scales. Here, we combined layer-specic Cre-driver mouse lines, optogenetics, and fMRI with electrophysiological recordings to identify distinct M1 layer-specic networks. Neuronal activities initiated in L2/3 were mainly conned within M1, while stimulation of L4, L5, and L6 evoked distinct responses in M1 and motor-related subcortical regions, including the striatum and motor thalamus. Although motor cortex has long been considered agranular (without L4), our results structurally, functionally, and neurovascularly conrm the presence of L4. We also nd that layer-specic fMRI responses closely couple with laminar electrophysiological recordings. Overall, our results elucidate distinct brain-wide neural archetypes of M1 layer-specic cortical circuits that provide important insights in uncovering the motor system architecture. advanced Parkinson's disease.


Introduction
Primary motor cortex (M1) plays a prominent role in motor behavior and the cortical control of movements through large-scale networks. Long-range inputs converge onto M1 through cortico-cortical, thalamo-cortical and neuromodulatory projections, while M1 outputs project to spinal/bulbar motor centers, striatum, thalamus, subthalamus, red nucleus, and pons 1,2 . Like the sensory cortices, M1 consists of a stack of interconnected but distinct layers. Traditionally, layers 2/3 (L2/3) are considered to be a prominent source of intracortical excitation of L5 3,4 . Additionally, L5 has output connections to the ventral anterior thalamus (VA), striatum and spinal cord, while L6 has output connections to the ventral lateral thalamus (VL) 3 . Recently, L4 has been proposed to have prototypical synaptic circuit connectivity 5 .
Despite this seemingly well-de ned structural organization, the brain-wide functional in uences of M1 layer-speci c pathways have yet to be elucidated at the whole-brain scale.
Historically, it has been di cult to disentangle the functional properties of different cortical layers as they are highly anatomically intermingled. However, some neurons that constitute layer-speci c pathways share relatively distinct neurochemical identities, suggesting that the activation of such discretized layers may result in highly dissimilar brain-wide responses. In particular, some neurons in L2/3, L4, L5, and L6 of the neocortex express the dopamine receptor D3 (Drd3), the sodium channel non-voltage-gated 1 alpha (Scnn1a), the RNA polymerase II subunit B32 (Rpb4) and the neurotensin receptor type 1 (Nstr1), respectively. Advances in molecular genetics, such as the advent of Cre-recombinase driver lines 6 and optogenetic tools 7 , have thus made it possible to selectively express transgenes in the various neocortical layers. Several studies have exploited this capability to selectively excite speci c layers and measure the downstream effects using in vivo electrophysiology. For example, in the somatosensory cortex, L2/3 was shown to suppress L5 to enhance the selectivity and output range of downstream activity for nely-tuned cortical coding 8 . L4 activity of the somatosensory cortex was also shown to directly inhibit L5 to sharpen the spatial representation of L5 neurons 9 . In addition, L6 corticothalamic neurons in the visual cortex have been shown to activate L5 cortical outputs 10 and to regulate the strength of cortical responses 11 .
These studies have mainly characterized the downstream activities within the cortex, while only a few L5 studies have investigated long-range pathways such as the corticotectal projections in the visual system 12 and the corticofugal projection in the auditory system 13 . Furthermore, these studies have focused on the sensory cortices (auditory, somatosensory and visual cortices), which are structurally and functionally distinct from M1. Since layer-speci c M1 activities distinctly contribute to cortical control of movements, and since different layers in M1 have distinct contributions to motor-related dysfunctions, such as Parkinson's disease (PD), it is essential to characterize the downstream effects of layer-speci c M1 excitation. However, to date, there has been no direct evidence of differential effects of M1 layerspeci c pathways on brain-wide circuit function, especially at the macroscopic and mesoscopic scale.
Conventional functional MRI (fMRI) is performed at a macroscopic scale. However, recent developments in high-spatial resolution fMRI provide opportunities for in vivo measurements of mesoscale layerspeci c cortical responses 14 . Laminar fMRI can be used to tackle key cognitive neuroscience questions by distinguishing bottom-up and top-down cortical responses and examining the interactions between the two 15 . It can also provide a robust testbed for examining in uential theories of brain-wide circuit function by assigning speci c computational roles to different cortical layers 16 . However, the mesoscale layerspeci c fMRI representation of neuronal activity has not been established.
The integration of optogenetic stimulation with fMRI (ofMRI) has enabled the causal in uences of genetically de ned neuronal populations on downstream regions to be measured directly [17][18][19][20][21][22] . Here, we combined targeted optogenetic stimulation of M1 L2/3, L4, L5, and L6 with large eld-of-view highresolution fMRI to reveal the causal in uences of activities originating in each layer on brain-wide regions, including different M1 layers, the thalamus, and the caudate putamen (CPu). Subsequently, we characterized these downstream region activities locally and remotely using in vivo electrophysiology to delineate the neuronal underpinnings of the fMRI responses. Last, we explored neuronal origins of the laminar ofMRI responses.

Results
Layer-speci c M1 stimulations activate distinct brain-wide networks To selectively activate L2/3, L4, L5, and L6 of primary motor cortex (M1) in vivo, we used transgenic mouse lines expressing Cre-recombinase 9,23−25 under control of Drd3, Scnn1a, Rpb4 and Nstr1 receptor regulatory elements, respectively. An AAV5 virus was injected into the M1 to express the excitatory opsin ChR2 in Cre-positive neurons to enable selective layer-speci c optogenetic control of M1. Histological and immunohistochemical examination con rmed that ChR2-EYFP was localized to the neurons in their respective layers and their intra-cortical projections (Fig. 1). Speci cally, ChR2-EYFP expression was observed in L2/3 M1 neurons and their L5 projections for the Drd3 L2/3 Cre-line, in L4 M1 neurons and their L2/3 projections for the Scnn1a L4 Cre-line, in L5 M1 neurons and their L2/3 projections for the Rbp4 L5 Cre-line, and in L6 M1 neurons and their L4 projections for the Ntrs1 L6 Cre-line.
Whole-brain ofMRI was applied to determine layer-dependent spatiotemporal characteristics of brain-wide evoked responses driven by layer-speci c M1 stimulation. Optical pulses were delivered at 5 Hz (30% pulse width duty cycle; light intensity, 30-50 mW/mm 2 ) and a general linear model (GLM) was used to identify voxels signi cantly modulated during stimulation. fMRI activation maps ( Fig. 2A, p < 0.001, FDR corrected, n = 12) and blood oxygen level dependent (BOLD) signal pro les (Fig. 3A) show that layerspeci c M1 stimulations at 5 Hz activate distinct brain-wide networks. Note that BOLD signal pro les were extracted using atlas-based, anatomically-de ned regions of interest (ROIs; Fig. 2B  Since the CPu and VL were robustly activated during L5 and L6 stimulations, respectively, histology and immunohistochemistry of these regions were examined. Long-range M1 L5 and L6 projections were revealed in the CPu and VL ( Supplementary Fig. 1), respectively, supporting our fMRI results. Since stimulation frequency can evoke distinct local and brain-wide responses 18,21,22,26 , we also applied 10, 20, and 40 Hz layer-speci c M1 stimulations within the same animals to investigate frequency-dependent spatiotemporal characteristics. The responses evoked by these stimulation frequencies largely resemble responses evoked by the 5 Hz stimulation ( Supplementary Fig. 2). The extracted t-values had a decreasing trend with frequency in ipsilateral M1 during L4 stimulation, while both the extracted t-values and fraction of ROI positively modulated had an increasing trend with frequency in contralateral CPu and ipsilateral VL during L5 stimulation ( Fig. 2D; one-way ANOVA followed by trend analysis). For L6 stimulation, ipsilateral M1 responses had increasing trends with increasing stimulation frequency for tvalues and fraction of ROI positively modulated, while a decreasing trend with frequency was observed for fraction of ROI negatively modulated ( Fig. 2D; one-way ANOVA followed by trend analysis). For area under the curve of the BOLD signal, they showed an increasing trend in ipsilateral M1 during L5 and L6 stimulations, as well as ipsilateral CPu and VL during L5 stimulation ( Fig. 3C; one-way ANOVA followed by trend analysis). In addition, the BOLD signal in the ipsilateral M1 transitioned from negative to positive during 5 and 10 Hz L6 stimulations, but not 20 and 40 Hz stimulations ( Supplementary Fig. 3A). No evoked responses were observed in a naïve animal ( Supplementary Fig. 4), indicating that the observed responses were a direct consequence of ChR2 stimulation rather than heat induced artifacts or undesired light-induced activations 27,28 . Neuronal underpinnings of the brain-wide fMRI responses  Table 3).
No evoked responses were observed in the naïve animal ( Supplementary Fig. 7), indicating that the observed responses were direct consequences of ChR2 stimulation and not photovoltaic induced artifacts or undesired light-induced activations 29,30 . Layer-speci c fMRI responses and their neuronal origins To examine layer-speci c fMRI representation of neuronal activity, we rst analyzed the local fMRI responses during layer-speci c M1 stimulations. Local fMRI activation maps ( Fig. 6A; p < 0.001, FDR corrected, n = 12) and BOLD signal pro les extracted from atlas-based anatomically-de ned ROIs show that layer-speci c M1 stimulation activates distinct local responses (Fig. 6D). Interestingly, all fMRI responses exhibit an increasing trend along the cortical depth ( Fig. 6B; one-way ANOVA followed by trend analysis). Parts of the ROI positively modulated exhibit an increasing trend along the cortical depth during L4 and L5 stimulations, while parts of the ROI negatively modulated exhibit a decreasing trend along the cortical depth during L4 stimulation ( Fig. 6C; one-way ANOVA followed by trend analysis). BOLD signal (Fig. 6D) area under the curve exhibit an increasing trend along the cortical depth during L5 and L6 stimulations ( Fig. 6E; one-way ANOVA followed by trend analysis). Similar results were obtained during 10, 20, and 40 Hz layer-speci c stimulations ( Supplementary Fig. 8). Notably, during L6 stimulation, the BOLD signal shows initial negative response, which reduces with cortical depth (Fig. 6D) and higher stimulation frequency ( Supplementary Fig. 8D).
To understand the neuronal origins of such layer-speci c fMRI responses, we analyzed in vivo extracellular recordings obtained along the M1 cortical depth (Figs. 7 and 8). Surprisingly, the LFP amplitude decreased along the cortical depth during L2/3 and L5 stimulations, while the amplitude increased along the cortical depth during L6 stimulation ( Fig. 7B and 7C). The discrepancy between the fMRI and LFP response trend across layers during stimulation of each layers potentially point to layerspeci c neurovascular coupling mechanisms. Scatter plot shows the overall correlation between the LFP and the laminar BOLD-fMRI signal across all four cortical layer stimulations and resulting responses in all layers (Fig. 7D). The amplitude of the LFP responses generally decreased with stimulation frequency with all four layer-speci c stimulations ( Supplementary Fig. 9C). LFP responses at the beginning of the stimulation was similar across stimulation frequencies ( Supplementary Fig. 9A). We also examined neuronal spiking activity along the M1 cortical depth in response to M1 layer-speci c stimulations ( . The scatter plots shows the correlations between the spike rate, the laminar BOLD-fMRI signal, and the LFP along the cortical depth ( Fig. 8D and Supplementary Fig. 10).

Discussion
Here, we combined layer-speci c Cre-driver mouse lines, optogenetics, and large-view fMRI with subsequent electrophysiological recordings to reveal distinct brain-wide ( Fig. 2-5) and local ( Fig. 6-8) M1 layer-speci c networks. These techniques enable investigations of the mechanisms of both local and brain-wide layer-speci c cortical circuit functional changes during development, aging, diseases and pharmacological interventions in future studies. In particular, how different cortical layers contribute to the local networks, and how these local networks contribute to brain-wide functional networks in normal and diseased state should be addressed. For instance, our results show that laminar BOLD-fMRI (Fig. 6), LFP (Fig. 7), and spike (Fig. 8) responses were mainly con ned within the deeper layers (L5 and L6) of M1 during L2/3 stimulations while L5 stimulations evoked strong responses in all layers. Globally, our results show that BOLD (Figs. 2 and 3), LFP (Fig. 4), and spike (Fig. 5) responses were mainly con ned within M1 during L2/3 stimulations, while L5 stimulations evoked local (M1) and brain-wide (CPu and VL) responses. Therefore, where L5 stimulation appeared to generate widespread activation both within and outside of M1, L2/3 activation imposed local and global bounds on neural activation. Since L2/3 of the somatosensory cortex can suppress L5 8 , we speculate that L2/3 stimulation similarly limited activity in the deeper layers of M1 and inhibited downstream propagation to regions such as the CPu and VL. Future studies may address the role(s) of each layer in coordinating and gating large-scale motor networks.
Traditionally, M1 is an agranular cortical area which features a lack of cytologically distinct L4 31,32 . The existence of a M1 L4 has emerged from previous studies as one that is composed of a thin band of pyramidal cells 33,34 and with thalamic afferents innervated at the middle of M1 35,36 . These concepts are in contrast to the traditional view of M1 layer architecture and has only been recently challenged with the reporting of a functional L4 5,37 . Our orescence images show ChR2-EYFP expression in L4 neurons for the Scnn1a Cre-line mice in M1 (Fig. 1B -D), and the thickness of the expression in L4 was ~ 100 µm ( Fig. 1B -C) which is consistent with previous reports 5 . ChR2-EYFP expression was observed in L2/3 projections from L4 and is also consistent with the previous study 5 . From a functional perspective, our electrophysiological recordings demonstrate that M1 L4 responses can be evoked by L4 optogenetic stimulation, or by that of L2/3, L5 or L6 (Figs. 7 and 8). Activities initiated from M1 L4 can even propagate downstream to the CPu and VL (Figs. 4 and 5). Furthermore, these neuronal responses parallel the fMRI responses which are connected through neurovascular coupling (Figs. 2-5). Our histology, electrophysiological recordings, and ofMRI results support the presence of L4 in M1 structurally, functionally, and neurovascularly, which facilitates future studies to focus on what types of information M1 L4 neurons process, how they do so, and how this relates to motor behavior in healthy, development, aging, and disease states.
Parkinson's disease (PD) is a progressive neurodegenerative disorder that affects movement. Although striatal dopamine de ciency caused by neuronal loss in the substantia nigra has been one of the hallmarks of PD [38][39][40] , recent evidence demonstrated that the motor cortex plays a prominent role [41][42][43] . Furthermore, some studies indicate that M1 is fundamental to the pathophysiology of PD which propagates along the cortico-striatal-thalamo-cortical network 38,41,43 . Since motor cortex consists of different layers, it is imperative to study how different layers in M1 contribute to PD. Previous postmortem studies have shown that PD patients have decreased M1 cortical thickness with ~ 70% reduction in L1 and L3, and with 25% reduction in L6 41,44 . However, the layer-speci c M1 local and brain-wide circuit function alterations in PD remain largely unknown. Our results demonstrate that L5 and L6 M1 activities propagate to the CPu and VL (Fig. 2-5), which indicates that these layers directly in uence the corticostriatal-thalamo-cortical network. Additionally, we show that neuronal activity initiated in different layers results in distinct local networks (Fig. 6-8), yet how these local networks affect cortico-striatal-thalamocortical network remain unknown. Future studies may utilize layer-speci c ofMRI to study different PD animal models to understand and characterize local and brain-wide layer-speci c circuit dysfunction, and the interactions between the local and global networks in PD. Subsequent electrophysiological recordings can explore the neural origins of these networks, and ultimately establish foundations for early diagnosis and treatment of Parkinson's disease.
A better understanding of the contribution of the motor cortex to PD symptoms will facilitate the development of novel therapeutic approaches. Subthalamic nucleus (STN) deep brain stimulation (DBS) has been recognized as a therapeutic method for PD 45 , yet many patients are reluctant to undergo invasive procedures or are not eligible. Recent studies have explored non-invasive methods using transcranial magnetic stimulation [46][47][48][49] and pharmacological modulation 50 of the motor cortex. However, these methods have not been proven to be reliable, partly due to the non-speci c stimulation and modulation of the whole motor cortex. Nevertheless, these studies suggest that TMS and pharmacological modulation of the motor cortex can induce striatal DA release which attenuate PD symptoms. Moreover, 120 Hz stimulation of the M1 L5 cells attenuated PD symptoms in parkinsonian rodent model 51 . Our results demonstrate that L5 stimulations evoked striatal responses (Figs. [2][3][4][5]. Taken together with the previous studies, our study results show that L5 layer-speci c stimulations can potentially be used as a neuromodulation technique for PD treatment.
Motor cortex stimulation is also used as a treatment for chronic neuropathic pain [52][53][54] or for rehabilitation after stroke 55-57 , but with mixed success. One limitation of motor cortex stimulation is that it cannot separately modulate different layers or different neuronal elements within a region. This makes it di cult to identify the exact mechanisms of different motor cortex stimulation paradigms. By using optogenetics to selectively stimulate different layers of M1, we found that each layer drives unique local and brain-wide responses. Future studies interested in therapeutic applications may wish to investigate particular layers within this circuit. We also show that the stimulation of layer-speci c M1 can drive distinct effects depending on the precise frequency used. Frequency can be a key parameter in optimizing the e cacy of motor cortex stimulation, and our current and past studies illustrate this effect.
The sensory cortex is heterogeneous with differential functions that include visual, auditory and somatosensory processing, and similarly consists of a stack of six interconnected but distinct layers 31 . Although the current study mapped the propagation of layer-speci c M1 activities, the propagation of sensory cortical activities is yet to be characterized. Future studies may determine the speci c spatiotemporal propagation of layer-speci c sensory cortical activities using layer-speci c ofMRI. Furthermore, the mesoscale fMRI representation of neuronal activity can also be established in the sensory cortices, which can provide foundations for studying the alterations of local and brain-wide layerspeci c cortical networks in different neurodegenerative diseases. In addition, cognitive functions, such as memory acquisition, memory consolidation, and memory retrieval, involves in complex brain functional networks 58,59 . However, how different cortical layers are involved in these processes remains largely unknown. Layer-speci c ofMRI may be applied to characterize and explore the spatiotemporal dynamics and functional roles of these cognitive networks.
The nature of neurovascular coupling remains an active area of research 60,61 with various studies suggesting that positive BOLD signals re ect increases in synaptic input 62 and others linking the BOLD signal to spiking activity 17,18,63 . In the current study we examined the neurovascular coupling in response to layer-speci c M1 stimulation in core motor-related regions such as the cortex, striatum and thalamus. Overall, the temporally resolved coupling between neuronal spiking/LFP and the BOLD response was region speci c. We showed that positive BOLD responses in the ipsilateral M1 correspond to negative peaks of the LFP responses ( Supplementary Fig. 11) and the increase in spiking for all layer-speci c stimulations, which is consistent with various neurovascular coupling studies 17,18,60,62,63 . For the ipsilateral CPu, we found that a negative and positive BOLD response were respectively linked to positive and negative peaks of the LFP during L4 and L5 stimulations ( Supplementary Fig. 11), respectively. Speci cally, L6 stimulations evoked BOLD responses that transitioned from negative to positive which corresponded to a similar transition of the LFP response from positive peaks to negative peaks ( Supplementary Fig. 11). Furthermore, all layer-speci c stimulations resulted in a spiking increase in the ipsilateral CPu, even though L2/3 and L4 stimulation evoked no and negative BOLD responses, respectively. These ipsilateral CPu results suggest that the BOLD and LFP responses inversely covary, while changes in spiking might not be highly correlated to the BOLD response. These CPu results are supported by previous striatal studies which show that negative BOLD is associated with large increases in neuronal activities [64][65][66] . For the ipsilateral VL, on the other hand, we found that both positive and negative BOLD responses were linked to positive LFP peaks for L4, L5 and L6 stimulations (Supplementary Fig. 11; L2/3 stimulation resulted in no BOLD nor LFP responses). Unchanged or negative BOLD responses corresponded to a spiking increase in the ipsilateral VL during L2/3 and L4 stimulations, while positive BOLD corresponded to a mix of spiking increase and decrease in different units during L5 and L6 stimulations. While the neurovascular coupling properties varied by region in the current work, additional studies using simultaneous fMRI and electrophysiological recordings should be performed to elucidate region-and/or layer-speci c relationships.
More recently, neurovascular coupling research has expanded to using mesoscale laminar fMRI methods 14,67 . Among the M1 cortical layers in mice, L4 is considered the thinnest at ~ 100 µm and L5 and L6 the thickest at ~ 300 µm 5 . While resolving these thicknesses approaches the limitations of the spatial resolution of fMRI, we attempted to establish the mesoscale layer-speci c representation or neuron activity. Although layer-speci c fMRI responses were distinct during L2/3, L4, L5, and L6 stimulations, all fMRI responses increased along the cortical depth. This phenomenon, however, was not observed in the LFP and spike recordings. Despite the distinct M1 recordings that we observed during layer-speci c stimulations, LFP responses and increased spiking were detected across all cortical layers.
Since draining veins carry deoxygenated blood from deeper cortical layers to the surface 68 , L6 neurons would be the rst to consume the oxygen through the supply of nearby capillaries. This would leave deoxyhemoglobin in the blood stream to be transported upwards to the super cial layers, and may cause an initial negative or decrease in BOLD at the super cial layers. In particular, the BOLD responses transitioned from negative to positive in the super cial layers during low frequency (5 Hz -20 Hz) L6 stimulation, and the amplitude of this initial negative BOLD response decreased along the cortical depth. This observation supports the hypothesis that draining veins transports deoxyhemoglobin from the deeper layers to the super cial layers causing an initial negative or decrease in BOLD. Our results provide insight into the feasibility of studying neurovascular coupling at the mesoscale using ofMRI. Future studies may take into account the neurovasculature to elucidate the exact mechanisms of mesoscale layer-speci c neurovascular coupling. Carolina viral vector core (titer of 4 × 10 12 particles/mL). All stereotactic surgeries were performed with mice at 8 -10 weeks of age. Animals were anesthetized with iso urane (induction 5%, maintenance 1.5% -2%; Sigma-Aldrich) and secured in a stereotactic frame with nonrupturing ear bars (Kopf Instruments). A heating pad was used to maintain body temperature, and sterile ocular lubricant was applied to the eyes to prevent desiccation during surgery. Buprenorphine (1 mg/kg) was injected subcutaneously for analgesia. After a midline incision along the scalp, a small craniotomy and viral injection/cannula implantation were performed at the primary motor cortex (+ 0.75 mm AP [anterior-posterior], − 1.5 mm ML [medial-lateral], injection at + 0.15 mm, + 0.35 mm, + 0.6 mm, or + 0.8 mm DV [dorsal-ventral] for L2/3, L4, L5, and L6 respectively). 1.5 μL of the AAV5/DIO-ChR2-EYFP virus was delivered using a 10 μL syringe with a 34G metal needle (World Precision Instruments Inc.) at a 75 nL/min ow rate driven by a microsyringe pump controller. Following injection, the injection needle was held in place for 10 min before slowly retracting it from the brain. A custom-designed ber-optic cannula was mounted and secured on the skull using light-cured dental cement (Kuraray Inc.), with the optical ber extending from the cannula's base to the desired depth (0.1 mm above the injection site). Following surgery, mice were kept on a heating pad until recovery from anesthesia and were given carprofen (5 mg/kg, subcutaneously [s.c.]) daily for 2 days to minimize post-operative discomfort. All experiments were conducted at least 4 weeks following virus injection to ensure optimal ChR2 expression. Probe locations were validated in all animals used for ofMRI experiments with T2-weighed structural MRI images.

Histology and immunohistochemistry to validate expression of ChR2
To con rm the precise targeting of ChR2 to L2/3, L4, L5 or L6, a cohort of mice injected with the DIOrecombinant virus was deeply anesthetized with pentobarbital and transcardially perfused with 0.1 M PBS followed by ice-cold 4% paraformaldehyde (PFA) in PBS. Brains were extracted and xed in 4 % PFA overnight at 4 °C. The brains were equilibrated in 20 % sucrose in PBS at 4 °C overnight. Axial sections of 50 μm thickness were prepared on a freezing microtome (Micron HM550 Microtome, Thermo Scienti c Inc.). For immunohistochemistry, free oating sections were washed with 0.1 M PBS for 20 mins at room temperature. Sections were then exposed to 200 ng/ml of DAPI in PBS at room temperature for 20 mins. Slices were then washed and mounted using Fluoromount-G (SouthernBiotech, Birmingham, AL). Immuno-uorescence was assessed with a laser confocal microscope (LSM 880 inverted confocal with Airyscan, Zeiss Inc.) at the Cell Sciences Imaging Facility at Beckman Center for Molecular and Genetics Medicine. The images were then analyzed to validate the expression of ChR2. Since the primary motor cortex is 1000 µm thick, and since the AAV5 spread is roughly 500 µm, region of interest ( Figure 1A) was de ned as a 500 µm × 1000 µm (width × height) rectangular box centered at the injection site. Anatomical landmarks including the midline, corpus callosum, cortex and striatum were identi ed and the ROI was de ned based on the mouse brain atlas and injection site location to ensure consistency. For sensitivity, 12%, 15%, 23% and 15% of cells identi ed within the injection area were ChR2-EYFP-positive for Drd3-Cre (L2/3), Scnn1-Cre (L4), Rpb4-Cre (L5), or Ntsr1-Cre (L6) mouse lines, respectively. All mouse lines have 100% speci city.

ofMRI Experiments
All MRI experiments were carried out on a 7 Tesla Bruker-converted Biospec small animal MRI system at the SCI3 facility at Clark Center using a transmit-only birdcage coil in combination with a customdesigned actively decoupled receive-only surface coil. Animals were initially anesthetized in an induction chamber with 5% iso urane before placement into the restraining apparatus. The restrained animal was then placed onto the MRI-compatible cradle with ears, teeth, and head secured. To maximize signal-tonoise ratio, the receiver coil was rst placed on top of the head and centered over the brain and then placed into the iso-center of the magnet. The animals were lightly anesthetized and sedated using a combination of iso urane and dexmedetomidine (0.25% iso urane mixed with O 2 and medical air; an initial bolus of 0.1 mg/kg dexmedetomidine; continuous infusion at 0.1 mg/kg/hr during scanning). During all fMRI experiments, continuous physiological monitoring was performed using an MRIcompatible system (SA Instruments). Vital signs were within normal physiological ranges (rectal temperature: 36.5 -37.5 °C, heart rate: 260 -420 beat/min, breathing: 80 -120 breath/min, oxygen saturation: > 90%) throughout the duration of the experiments.
Fourteen contiguous 1.0-mm slices were positioned in the transverse orientation according to the mouse brain atlas to cover the whole brain. T2-weighted high-resolution anatomical images were acquired prior to fMRI to check for brain damage and con rm accurate probe location. These anatomical images were acquired using a rapid acquisition with relaxation enhancement (RARE) sequence with eld of view (FOV) = 20 × 20 mm 2 , matrix = 256 × 256, RARE factor = 8, echo time (TE) = 33 ms, repetition time (TR) = 4,000 ms. ofMRI measurements were obtained for the same slices using a single-shot Gradient-Echo Echo-Planar-Imaging (GE-EPI) sequence with FOV = 20 × 20 mm 2 , matrix = 75 × 75, ip angle = 40°, TE = 14 ms, TR = 1,000 ms.
The MRI scanner and laser for optogenetic stimulation were synchronized using a Master pulse stimulator system (A.M.P.I.). The light delivery system was kept outside the magnet and long optical cables (~5 m) were used to deliver light into the scanner bore. Blue light was delivered using a 473-nm DPSS laser measured before scanning as 3 mW at the ber tip (250 μm) corresponding to a light intensity of 60 mW/mm 2 for optogenetic stimulation. To determine the frequency-dependent spatiotemporal characteristics of evoked M1 layer-speci c responses, four frequencies were used (5 Hz, 10 Hz, 20 Hz, and 40 Hz with 30% duty cycle). M1 layer-speci c neurons were stimulated with a standard block design paradigm that consisted of 10-s light-on and 30-s light-off periods. Three to four trials were recorded for each frequency in a pseudo-random manner in each animal.

fMRI Data Analysis
For each animal, all fMRI images were corrected for slice timing differences, realigned to the rst image of the rst fMRI session, and spatially smoothed with a Gaussian kernel (FWHM = 2 pixels) using SPM12 (Wellcome Department of Imaging Neuroscience, University College, London). The fMRI data were then temporally linear detrended and temporally band-pass ltered (0.01-0.25 Hz) to eliminate baseline drift caused by system instability and physiological noise. For detrending, linear trend was rst calculated using linear regression of the global temporal signal obtained from the whole brain. Then, detrending was applied to the temporal signal of each voxel using the obtained linear trend. T2-weighted images from each animal were registered to a custom-made brain template acquired with the same settings.
Registration was performed by a ne transformation and Gaussian smoothing to maximize normalized mutual information (SPM12). Only then were the fMRI images resliced correspondingly. A general linear model (GLM) was applied. Student's t test was performed to identify activated voxels using the threshold p < 0.001. The mapped responses were then compared between each animal to assess result quality before group level analysis. At the group level, realigned, registered, and resliced images corresponding to the same fMRI session were averaged across all animals. To quantify the activations from the activation maps, t-values within the ROI were extracted and averaged, while fraction of ROI positively and negatively modulated were calculated with a voxel-wise threshold of p < 0.05 as positively or negatively modulated. BOLD temporal pro les for each stimulation frequency were extracted from identical ROIs delineated from the mouse brain atlas by regional averaging across the whole ROI.
Area under the curve of the BOLD pro les was calculated to quantify the strength of the BOLD responses for each individual animal, and then averaged across animals within the same group.

In Vivo Electrophysiology
In vivo electrophysiology was performed to directly measure the neuronal activity of the ipsilateral M1, CPu and VL. Similar to ofMRI experiments, the animals were lightly anesthetized and sedated using a combination of iso urane and dexmedetomidine (0.25% iso urane mixed with O 2 and medical air; an initial bolus of 0.1 mg/kg dexmedetomidine; continuous infusion at 0.1 mg/kg/hr during experiment).

Electrophysiology Data Analysis
For LFP recordings, raw data were initially imported into custom-written Matlab software and a notch lter centered at 60 Hz was applied to remove power-line noise. Individual animal LFPs were averaged across stimulation blocks to generate a single LFP trace for each stimulation frequency, which was then averaged across animals. For spike recordings, two-tailed paired t-tests were used to identify signi cant changes in ring rate within each unit between pre-stimulation and stimulation periods (10 seconds each). Only units with signi cant change (either increase or decrease) were included in the following analysis. Peri-event time histograms were averaged across stimulation blocks and units to generate an average peri-event time histograms for each stimulation frequency. One-way ANOVA followed by Bonferroni's post hoc test was applied to compare pre-stimulation, stimulation and post stimulation periods.
Although BOLD-fMRI and LFP responses re ects the overall signal within a region, which all recorded units should be considered when exploring the association between spiking, BOLD-fMRI and LFP, we still  Table 7), between these spiking parameters and long-range LFP recordings (Supplementary Table 8), between these spiking parameters and laminar BOLD-fMRI (Supplementary Table 9), and between these spiking parameters and laminar LFP recordings (Supplementary Table 10) generally showed that the spiking parameter, total absolute value of percentage spike rate change of all neurons, tends to yield the highest correlation.

AUTHOR CONTRIBUTIONS
RWC participated in planning the study, performed surgeries, fMRI experiments, electrophysiological recordings, immunohistochemistry, data analysis, and wrote the paper. MA, HJL, and BE helped with electrode design, and electrophysiological recordings. HA and DF participated in the original planning of the study and provided the transgenic mouse lines. JHL planned the study, advised all personnel involved in the study, and wrote the paper.