Early fMRI responses to somatosensory and optogenetic stimulation reflect neural information flow.

Blood oxygenation level-dependent (BOLD) functional magnetic resonance imaging (fMRI) has been widely used to localize brain functions. To further advance understanding of brain functions, it is critical to understand the direction of information flow, such as thalamocortical versus corticothalamic projections. For this work, we performed ultrahigh spatiotemporal resolution fMRI at 15.2 T of the mouse somatosensory network during forepaw somatosensory stimulation and optogenetic stimulation of the primary motor cortex (M1). Somatosensory stimulation induced the earliest BOLD response in the ventral posterolateral nucleus (VPL), followed by the primary somatosensory cortex (S1) and then M1 and posterior thalamic nucleus. Optogenetic stimulation of excitatory neurons in M1 induced the earliest BOLD response in M1, followed by S1 and then VPL. Within S1, the middle cortical layers responded to somatosensory stimulation earlier than the upper or lower layers, whereas the upper cortical layers responded earlier than the other two layers to optogenetic stimulation in M1. The order of early BOLD responses was consistent with the canonical understanding of somatosensory network connections and cannot be explained by regional variabilities in the hemodynamic response functions measured using hypercapnic stimulation. Our data demonstrate that early BOLD responses reflect the information flow in the mouse somatosensory network, suggesting that high-field fMRI can be used for systems-level network analyses.

Blood oxygenation level–dependent (BOLD) functional magnetic resonance imaging (fMRI) has been widely used to localize brain regions and networks associated with sensation, perception, and behavior (1–3). Different functional brain regions are connected through feedforward and feedback, ascending and descending, or bottom-up and top-down projections within the network (4). Therefore, it is critical to determine the direction of information flow to better understand brain functions. However, BOLD fMRI response, which is sensitive to vascular density and baseline physiological parameters (5), varies among brain regions and subjects (6, 7). For example, a region with large draining veins has a BOLD fMRI response that is delayed by a few seconds compared with a region that contains only capillaries within the parenchyma (8–10). Consequently, it has been argued that the order of neural events that occur within a few to tens of milliseconds of one another can be biased in fMRI dynamics due to regionally variable hemodynamic response functions (HRFs) (11).

One ultimate goal of fMRI research is to demonstrate the causality and temporal sequences of neural events in humans (12). Because the contribution of capillaries to BOLD fMRI increases with the magnetic field strength, we hypothesized that early hemodynamic responses at ultrahigh fields would reflect the timing of neural activation. During forepaw somatosensory stimulation in rats, the earliest fMRI response was observed in the thalamocortical (TC) input layer 4 (L4) within the primary somatosensory cortex (S1) (13–16), suggesting that early fMRI signals reflect synaptic input. However, differences in the BOLD onset responses between layers or regions could be related to differences in their HRFs. Thus, a systematic study of different neural processing orders is crucial to determine whether the onset times of BOLD fMRI responses indeed follow the order of neural events in functional interconnected regions that include thalamic nuclei with potentially different HRFs.

To investigate whether early BOLD response timing reflects the direction of neural information flow, we performed high spatiotemporal–resolution fMRI at an ultrahigh magnetic field of 15.2 T while conducting somatosensory stimulation, optogenetic stimulation, and a vascular challenge in lightly anesthetized mice. Forepaw somatosensory stimulation induced significant BOLD fMRI responses across multiple interconnected brain regions, including the ventral posterolateral nucleus (VPL), the posterior complex of the thalamic nucleus (PO), S1, and the primary motor cortex (M1) (17), among which the expected information flow is VPL → S1 → M1 (18, 19). To investigate whether the early BOLD responses were reversed when the activation sequence was reversed, we performed optogenetic stimulation of excitatory neurons in M1. In addition, to investigate whether the differences in onset times among brain regions were driven by different HRFs, we used hypercapnic challenge–induced vasodilation as a vascular control condition (20, 21). The dynamic characteristics of the BOLD responses were compared among the active somatosensory regions and among the cortical layers of S1. We found that the order of onset times among the active regions and layers clearly coincided with the known sequence of neural activation, indicating that early BOLD fMRI responses can be used to identify the direction of neural information flow.

Results

Order of Early BOLD Responses in Somatosensory Regions Reflects the Known Sequence of Neural Activation.

The mouse somatosensory network, which contains feedforward TC and feedback corticothalamic (CT) and corticocortical (CC) projections, has been investigated extensively using electrophysiology (19, 22), optical imaging and anatomical circuit tracing studies (23–25). Forepaw or whisker stimulation induces neuronal activity in the relay thalamic nucleus (26) that then projects to L4 in S1 (26, 27). From there, the activity spreads to layer 2/3 (L2/3) and layer 5 (L5) and then projects to M1 (23, 28). To examine whether fMRI can detect that information flow, we acquired gradient-echo (GE) BOLD fMRI data from lightly anesthetized mice during somatosensory stimulation; we used an ultrahigh magnetic field of 15.2 T, a temporal resolution of 250 ms, and a spatial resolution of 156 × 156 × 500 μm3 (Fig. 1). Three coronal slices containing the primary somatosensory forelimb area (S1FL), M1, the secondary somatosensory cortex (S2), and thalamic areas were selected based on scout whole-brain fMRI studies (). Somatosensory stimulation was applied to a single forepaw in the form of repeated electric pulses with a current intensity of 0.5 mA, a pulse width of 0.5 ms, and a frequency of 4 Hz (29). Each fMRI trial consisted of 40 s prestimulus, 20 s stimulus, 60 s interstimulus, 20 s stimulus, and 60 s post-stimulus periods.

Fig. 1.

High spatiotemporal BOLD fMRI responses during forepaw stimulation in anesthetized mice reflect bottom-up processing in the somatosensory network. Well-known forepaw somatosensory stimulation was used to activate regions in the somatosensory network and investigate whether the temporal order of BOLD responses reflects the neural processing sequence. (A) Schematic of high temporal–resolution somatosensory fMRI. Three coronal EPI slices containing cortical and thalamic areas of the somatosensory network were acquired at 15.2 T with a TR of 250 ms during 20 s forepaw stimulation to measure bottom-up somatosensory processing (VPL → S1FL → M1). (B) Functional maps of one representative animal (i) and group-averaged data from seven animals (ii) were overlaid on the original EPI images. The somatosensory network contralateral to the stimulated side was activated (green arrow: S1FL; yellow arrow: PO; purple arrow: VPL). (C) Based on the somatosensory fMRI map, five different ROIs were defined using the Allen mouse brain atlas. (D) fMRI time courses of group-averaged data were obtained from the ROIs contralateral to the somatosensory stimulation in the three cortical (i; S1FL, blue; M1, red; S2, green) and two thalamic regions (ii; VPL, purple; PO, yellow). The BOLD responses follow the two 20 s stimulation periods. (E) Mean signal changes were obtained from averages between 6 and 20 s after the onset of stimulation. Clearly, the S1FL response is the highest during forepaw stimulation. (F) The fMRI time courses during two stimulus blocks were averaged (black) and fitted by double gamma variate functions (S1FL, blue; M1, red; S2, green; VPL, purple; PO, yellow) to determine the HRFs (i for the representative animal shown in B, i; ii for group-averaged data from the seven animals shown in B, ii). The fitting curves were normalized to the peak (iii; showing from 0 to 5 s after stimulus onset), and the times to reach 5% and 30% of the peak were then measured in all ROIs (iv). The order of early BOLD responses to somatosensory stimulation is VPL → S1FL → PO/M1/S2. Gray horizontal bar, 20 s stimulus duration; error bars in E and F, iv and gray shades in F, ii: SEM; colored circles in E and F, iv: individual animal data; *P < 0.05 and **P < 0.01 in E and F, iv (n = 7; one-way ANOVA with repeated measures followed by Tukey’s post hoc test).

High-quality single-shot echo planar images (EPI) of mouse brains were used for coregistration among animals and as the underlay for the fMRI maps (Fig. 1). Radial dark lines within the cortex in the high-resolution GE–EPI images indicate venous vessels. Somatosensory stimulation induced significant BOLD responses in the somatosensory network contralateral to the stimulated forepaw (Fig. 1 , i for one animal; Fig. 1 , ii for group data from seven mice). Only voxels showing significant responses to the stimulation (uncorrected P < 0.001 in Fig. 1 , i; corrected P < 0.05 in Fig. 1 , ii) were color coded. The highest change was observed in S1FL (green arrow in Fig. 1 , i), with smaller changes visible in S2. Two distinct foci were detected in the thalamus; the ventral and lateral focus (purple arrow in Fig. 1 , i) is VPL, and the dorsal and medial site (yellow arrow in Fig. 1 , i) is PO. During the somatosensory stimulation, BOLD responses peaked around the end of the 20 s stimulation (Fig. 1 , i and ii), and the average signal change in S1FL was significantly higher than those in other areas (Fig. 1 and for individual animal data).

To measure the dynamic properties of BOLD responses, we averaged the time courses of 320 data points from two stimulus blocks (Fig. 1 , i and ii for fMRI time courses from the individual shown in Fig. 1 , i and group-averaged data shown in Fig. 1 , ii, respectively) and fitted them using two gamma variate functions widely used to determine HRFs from noisy time series data (7, 14, 30–32). Individual animal data with the best-fitted curves are presented in , and HRF parameters, such as the peak time and full width at half maximum (FWHM), are reported with the goodness of fit in . Since the peak time and FWHM are sensitive to all sizes of venous vessels, early dynamic characteristics were determined before the contribution of large draining veins is dominant. Thus, the fitted HRFs were normalized to the peak to determine the times to reach 5% and 30% of the peak in each subject (Fig. 1 , iii and iv). The time to 5% of the peak was used to capture the onset time of hemodynamic responses without significant contribution from the drainage of downstream vessels (gray dashed arrow in Fig. 1 , iii). During somatosensory stimulation, the onset time of the response in VPL was the fastest, followed by that in S1FL, and the responses of M1, PO, and S2 were slower (Fig. 1 , iv and for individual data). When the averaged time courses of seven mice were used for fitting, the onset times were 0.62, 0.95, 1.28, 1.29, and 1.30 s for VPL, S1FL, M1, PO, and S2, respectively. Similar trends were consistently observed in time to 30% and time to 50% calculations (). The order of early BOLD responses among interconnected regions followed the known order of neural activation within the somatosensory network: VPL → S1FL → M1/PO.

Order of Early BOLD Responses Can Decode the Sequence of Neural Activation during Optogenetic Stimulation.

To ensure that our previous finding that the order of BOLD responses matched the known neural activity sequence was not driven by innate differences in the HRFs of the somatosensory network, we initially performed hypercapnic stimulation to induce BOLD responses without evoking neural activity. The BOLD responses evoked by a 20 s hypercapnic challenge were widespread, and the onset times of the three cortical regions of interest (ROIs) (2.25 to 2.59 s) were significantly faster than those of the two thalamic nuclei (3.03 to 3.17 s) (). Although the hypercapnic responses were slower than that of forepaw activation (), one interesting observation is that the hypercapnic responses of M1 (2.25 s) and PO (3.03 s) are slightly faster than those of nearby S1FL (2.59 s) and VPL (3.17 s), respectively, whereas the responses of somatosensory-evoked S1FL and VPL were faster than those of M1 and PO. This suggests that the difference in neural activity–induced BOLD onset times cannot be explained by differential blood delivery.

Next, we performed optogenetic stimulation to reverse the order of information flow. Optogenetic stimulation in M1 of transgenic Thy1-ChR2 mice was used to elicit the M1 → S1FL → VPL sequence within the somatosensory network. The location of the optic fiber in M1 was identified in the standard anatomical GE image (fast low angle shot imaging, FLASH) and the functional EPI (Fig. 2 for one representative animal). Because image artifacts caused by the optic fiber implant are minimized in EPI, BOLD fMRI was performed while channelrhodopsin-expressing excitatory neurons in M1 were activated by a 20 s photostimulation (Fig. 2; 473 nm, 20 Hz, 10 ms duration, 2 mW at the fiber tip, and n = 6 Thy1-ChR2 mice).

Fig. 2.

High spatiotemporal BOLD fMRI responses to optogenetic excitation of excitatory neurons in M1 measure CC and CT projections. The optogenetic stimulation of M1 excitatory neurons was performed to reverse the neural information flow compared with somatosensory stimulation. (A) Anatomic (FLASH) and functional (EPI) MRI images (15.2 T) of the brain with an optical fiber targeting M1 in transgenic Thy1-ChR2 mice. The image artifacts caused by the optic fiber implant were minimal in the original EPI and FLASH. (B) Schematic of optogenetic fMRI for M1 excitation to induce top-down somatosensory processing (M1 → S1FL → VPL) in reverse order to that from forepaw stimulation. (C) Functional maps of one representative animal (i) and group-averaged data from six animals (ii) were overlaid on the original EPI images. The optogenetically driven brain activations were in the ipsilateral M1, S1FL, S2, and multiple thalamic areas as well as the ipsilateral dorsal striatum and the contralateral M1. (D) fMRI time courses of group-averaged data were obtained from ROIs ipsilateral to the M1 stimulation (S1FL, blue; M1, red; S2, green; VPL, purple; PO, yellow). The fMRI time courses follow the two 20 s stimulation periods. (E) Mean signal changes were obtained from averages between 6 and 20 s after the onset of stimulation. The responses of S1FL and M1 are significantly higher than those of the other ROIs during the M1 stimulation. (F) The fMRI time courses during two stimulus blocks were averaged (black) and fitted by triple gamma variate functions (S1FL, blue; M1, red; S2, green; VPL, purple; PO, yellow) to estimate the HRFs for an initial rise, subsequent steady increase to peak, and falling shapes (i for the representative animal from C, i; ii for group-averaged data from the six animals shown in C, ii). The fitting curves were normalized to the peak (iii; showing from 0 to 5 s after stimulus onset; see also for the ipsilateral striatum and the contralateral M1), and the times to reach 5% and 30% of the peak were then measured in all of the ROIs (iv). The order of early BOLD responses to the M1 stimulation is M1 → S1FL → S2/PO/VPL. Gray horizontal bar, 20 s stimulus duration; error bars in E and F, iv and gray shades in F, ii: SEM; colored circles in E and F, iv: individual animal data; *P < 0.05 and **P < 0.01 in E and F, iv (n = 6; one-way ANOVA with repeated measures followed by Tukey’s post hoc test).

The BOLD signal changes induced by M1 stimulation were observed in the ipsilateral M1, S1FL, S2, and multiple thalamic areas as well as in the ipsilateral dorsal striatum and the contralateral M1 (Fig. 2 , i for one animal and Fig. 2 , ii for group data from six mice). Because these ipsilateral active sites overlap with areas that are functional during somatosensory stimulation, we examined whether the onset time of early BOLD responses reflected the top-down somatosensory processing evoked by optogenetic stimulation of M1. In the five predefined ipsilateral ROIs of the somatosensory network, BOLD responses to the optogenetic M1 excitation peaked around the stimulus offset (Fig. 2 , i and ii), and the average BOLD signal changes in M1 and S1FL were higher than those in other areas (Fig. 2 and for individual data). Overall, BOLD responses were higher and faster than those to forepaw stimulation ().

The dynamic properties of BOLD responses were measured using a curve fitting process (Fig. 2 , i and ii show the fMRI time courses from the individual and group-averaged data in Fig. 2 , i and ii, respectively). In the fitted curves of animal-wise averaged time courses, the onset times were 0.57, 0.68, 0.97, 0.97, and 1.04 s for M1, S1FL, S2, PO, and VPL, respectively (Fig. 2 , iii). In general, the onset response of M1 was evoked first, followed by S1FL, with the responses of PO, S2, and VPL being delayed (Fig. 2 , iv and for individual data). The onset times of the ipsilateral striatum and contralateral M1 were between those of S1FL and S2/PO (). The order of BOLD onset responses to optogenetic M1 stimulation (M1 → S1FL → PO/VPL) was the reverse of the response to somatosensory stimulation (VPL → S1FL → M1/PO), indicating that the order of early BOLD responses faithfully depicts the neural activation sequence in the functional network.

Dynamics of Layer-Dependent BOLD in S1FL during Neural Stimulation versus Hypercapnia.

Because the temporal order of early BOLD responses reflected information flow among long-range somatosensory regions, our next question was whether the neural activation sequence could be identified by the layer-dependent onset times. In S1FL, neurons in L4 are known to respond first during forepaw stimulation due to TC inputs from VPL (33, 34), whereas neurons in L2/3 respond first during optogenetic stimulation in M1 due to CC inputs from M1 (33, 35). To measure spatiotemporal changes at a laminar resolution, we flattened the cortical section containing S1FL and MI (Fig. 3) and generated fMRI maps from group-averaged images. To quantify the dynamics and magnitude of BOLD changes, S1FL was subdivided into three evenly spaced ROIs by depth: an upper cortical ROI containing L2/3, the middle layers (L4/5), and the lower layer (L6) (Fig. 3).

Fig. 3.

Cortical depth-dependent BOLD responses in S1FL during neural stimulation and hypercapnia. Layer-dependent dynamic responses were examined to determine whether the order of onset times of the BOLD responses followed expected neural circuits. (A) The cortical areas (S1FL and M1) were linearized using radially projecting lines perpendicular to the cortical edges. Note that the fMRI data were spatially normalized in a common brain space (underlay, study-specific EPI template; overlay, Allen mouse brain atlas). S1FL was then divided into three evenly spaced depths to quantify layer-specific responses (upper S1FL, red; middle, green; lower, blue). (B) Dynamic BOLD percent change maps were calculated with an average window of 2 s duration at an interval of 1 s from group-averaged data (n = 7 for forepaw stimulation; n = 6 for M1 excitation; n = 7 for hypercapnia) (Movies S1–S3). During forepaw stimulation, BOLD responses first appeared in the middle cortical layer of S1FL at 2 s after stimulus onset (yellow arrow), gradually propagated to all the layers of S1FL at 5 s, and later peaked in the upper layer. On the other hand, the early responses during M1 excitation occurred in the upper areas of M1 and S1FL (white arrows at 2 s) and spread downward over time. Hypercapnic-induced responses started in the upper layer of M1 (red arrow at 2 s), moved to the upper layer of S1FL at 5 s, and then propagated to the deeper cortical layers. (C) Mean signal changes during forepaw stimulation, optogenetic M1 excitation, and hypercapnia monotonically decreased with depth (i). However, the laminar BOLD responses normalized by the response at deeper cortical ROIs were consistent under all three experimental conditions (ii). (D) The group-averaged fMRI time courses in the S1FL layer-specific ROIs were fitted using double (or triple) gamma variate functions (upper S1FL, red; middle, green; lower, blue) to determine the HRFs. (E and F) The fitting curves were normalized to the peak (showing from 0 to 3 s after stimulus onset), and the times to reach 5% and 30% of the peak were then measured in the S1FL layer-specific ROIs. The order of early BOLD responses was middle → upper → lower for forepaw stimulation and upper → middle → lower for M1 excitation and hypercapnia. Gray horizontal bar, 20 s stimulus duration; error bars in C and F, SEM; colored circles in C and F, individual animal data; *P < 0.05 and **P < 0.01 in C, i and F (one-way ANOVA with repeated measures followed by Tukey’s post hoc test).

The cortical depth–dependent responses in S1FL and M1 over time were visualized using time-dependent percent change maps with a moving average window of 2 s duration at an interval of 1 s (Fig. 3 and Movies S1–S3). During forepaw stimulation, BOLD response within S1FL was evoked first at the middle layer of the cortex at 2 s (yellow arrow in Fig. 3) and then gradually propagated to the surface of the cortex and the deeper cortical layers. Later, the highest response was observed in the upper layer of S1FL, where large draining veins are concentrated (Fig. 3 and Movie S1). The fMRI response evoked by M1 optogenetic excitation was observed first in the upper cortical areas of M1 and S1FL at 2 s (white arrows in Fig. 3 and Movie S2), and it propagated downward over time, although the largest signal was maintained in the upper cortical area. Similarly, the hypercapnia-induced BOLD response began at the cortical surface in M1 (red arrow in Fig. 3) and then shifted to the upper cortical layer of S1FL before penetrating into the deeper cortical layers (Fig. 3 and Movie S3), which was consistent with the expected penetrating blood flow direction from the cortical surface to the deeper cortical layers (36–38). These results suggest that the laminar BOLD responses to neural stimulation and the hypercapnic challenge have different origins.

The average BOLD signal changes during stimulation were highest in the upper layers and decreased monotonically with cortical depth (Fig. 3 , i and 3 and for individual data). To detect whether relative BOLD responses across layers differed between neural stimulation and the hypercapnic challenge, the laminar BOLD responses were normalized by the response at a deep cortical depth, and the normalized BOLD responses were consistent across all experimental conditions (Fig. 3 , ii). Similar depth-dependent BOLD magnitude changes were previously observed in rats (13, 39), cats (40, 41), monkeys (42), and humans (43–45). Overall, the peak BOLD magnitude changes cannot identify the layers receiving synaptic input.

The dynamic properties of BOLD responses in layer-dependent ROIs at three different depths were measured and analyzed using the gamma-variate fitting process to investigate whether the early BOLD response can identify the layers receiving synaptic input from TC and CC projections (Fig. 3 for averaged data; for individual ROIs). During forepaw stimulation, the onset response of the middle cortical layers tended to be slightly faster than that of the upper and lower layers of S1FL (Fig. 3 and for individual data). The onset times of the averaged time courses were 0.85, 0.98, and 1.06 s for the middle, upper, and lower layers, respectively (Fig. 3). These results were consistent with previous findings in the rat S1FL during forepaw stimulation (13, 14, 16).

Because dynamic differences between layers could have resulted from different neural inputs or different HRFs, we examined data from the optogenetic stimulation of M1 and hypercapnia conditions to separate the two potential contributions. The S1FL response driven by M1 excitation was first observed in the upper layer, followed by the middle and lower layers (Fig. 3 and for individual data). The onset times of those averaged time courses were 0.55, 0.64, and 0.76 s for the upper, middle, and lower layers, respectively (Fig. 3). Indeed, the CC projection from M1 produced the earliest BOLD response in the upper cortical layers in S1FL.

Hypercapnia-induced hemodynamic responses showed a tendency to begin in the upper layers and move to the middle and lower layers (Fig. 3 and for individual data). The onset times of the hypercapnic BOLD responses derived from the averaged time courses were 2.42, 2.59, and 2.76 s for the upper, middle, and lower layers, respectively (Fig. 3). This onset time difference can be explained by arterial blood delivery along the cortical depth (38). The upper cortical layers have a faster onset time for arterial blood delivery and shorter path lengths between the arterioles and venules than the deeper cortical layers. Consequently, the BOLD signals from the upper cortical layers have faster onset times than those from the deeper layers. With dynamic observations of neural activity–driven and hypercapnic BOLD responses, we can conclude that the synaptic input layers receiving TC and CC projections have the earliest BOLD responses.

Discussion

fMRI Onset Times of Somatosensory Brain Network Follow the Neural Information Flow.

In this work, we have demonstrated that the spatiotemporal dynamics of BOLD fMRI signals follow the neural activation flow, indicating that ultrahigh-resolution BOLD fMRI at ultrahigh fields can be used to determine the direction of information flow in a whole brain.

The putative direction of information flow for somatosensory processing was estimated from the BOLD onset times (Fig. 4): 1) The relay thalamic nucleus VPL responds; 2) L4 of S1FL receives the TC synaptic input (Δt = ∼230 ms) and propagates it to other layers (Δt = 130 to 210 ms); 3) L2/3 and L5 neurons in S1FL project to M1 (23, 28), and L5 neurons in S1FL project to the higher-order thalamic nucleus PO (18, 26, 46); and 4) S2 has a slower BOLD onset time than S1FL (Δt = ∼320 ms), but the vascular reaction times of S1FL and S2 during hypercapnia are quite similar. The neurons in S2 receive three different synaptic inputs: directly from VPL (47), directly from S1, and indirectly from S1 via PO (18, 46). The slow BOLD response in S2 indicates that the CC projections from S1FL might predominantly induce hemodynamic responses.

Fig. 4.

Putative order of functional activity in the somatosensory network as measured by BOLD fMRI. The orders of bottom-up and top-down processing in the somatosensory network were determined using the relative difference in the onset times of averaged BOLD responses between somatosensory regions. The arrows indicate information flow direction, and the difference in relative BOLD onset times between regions (Δt) is also reported. (A) Sensory-evoked signal from forepaw stimulation is first relayed through VPL into the middle layer of S1FL. The thalamocortical inputs next propagate to the upper and lower S1FL, and then the upper S1FL projects to M1 and S2 through the CC pathway, and the lower S1FL projects to the higher-order thalamic nucleus PO. (B) Neural activity in M1 following optogenetic stimulation is relayed into the upper layer of S1FL and then to S2. The signals propagated along the cortical layers are subsequently projected to PO and VPL. Solid line, expected functional pathway; gray dashed line, anatomically possible pathway.

Optogenetic stimulation of M1 induced activities in efferent pathways (Fig. 4). According to anterograde tracing studies from M1 (48, 49), neurons in L2/3 and L5a project mostly to the ipsilateral S1, S2, and striatum and to the contralateral M1, whereas neurons in L5b and L6 have extensive projections to the thalamus (49, 50). Within the somatosensory network, neurons in both S2 and PO can receive indirect input from M1 via S1, whereas neurons in VPL receive M1 inputs only indirectly via L6 of S1FL (51). Those previous findings are consistent with the results of our BOLD fMRI study here, wherein M1 stimulation induced activity first in L2/3 of S1FL, then in S2 and PO, and last in VPL.

Our fMRI results have shown that the sequential order of BOLD onset times (VPL → S1FL L4 → S1FL L2/3 → M1 for somatosensory stimulation; M1 → S1FL L2/3 → S1FL L6 → VPL for optogenetic stimulation) coincide with the order of expected neural circuit connections, indicating that the BOLD times capture local capillary responses near active neurons. BOLD sensitivity to vessel size depends on the magnetic field, increasing with magnetic field strength linearly for emerging venules and veins and supralinearly for capillaries and small venules (52). Therefore, higher magnetic fields increase the relative contribution that microvessels near sites of neural activity make to the BOLD signals (52–54). Due to the increased sensitivity to capillaries, the early BOLD response at ultrahigh fields could be more likely to originate from the capillaries before downstream draining becomes significant. Consequently, ultrahigh field fMRI is needed to measure the early BOLD responses with enough precision to detect the information processing order.

The Earliest Onset Times of BOLD fMRI within the Cortex Occur at Synaptic Input Layers.

Many studies using well-established animal models have investigated the feasibility of using dynamic BOLD fMRI to obtain information from the synaptic input layers (9, 13–16, 55–57). Dynamic BOLD changes in the S1FL in response to forepaw stimulation of rats anesthetized with α-chloralose were previously investigated at 7, 9.4, and 11.7 T (13, 14, 16, 55). Silva et al. found that the onset time was about 1.1 s in deep cortical areas and ∼1.5 s on the cortical surface (55). To improve the temporal resolution in previous single-slice fMRI studies (13, 14), a single k-space line was acquired in repeating fMRI runs with a 40 to 50 ms temporal resolution and phase-encoding steps. In that way, the onset times were found to be 0.59 s for L4/5, 1.28 s for L2/3, and 1.11 s for L6 (13). Yu et al. reported that layer-dependent onset times measured using a line scanning approach were 0.75, 0.96, and 0.92 s for L4, L2/3, and L5/6, respectively, and that the BOLD response in M1 was delayed by ∼200 ms relative to L4 of S1FL (14). Although the exact onset times depend heavily on BOLD sensitivity, temporal resolution, the strength of neural activity, and anesthesia, the common observation is that the BOLD response to somatosensory stimulation starts at L4 of S1FL, where baseline cerebral blood flow and blood volume are also the highest (37, 58, 59). In our study, forepaw somatosensory stimulation induced the earliest BOLD response at L4/5, whereas optogenetic stimulation of M1 produced the earliest BOLD responses in L2/3 of S1FL. Thus, the earliest onset of hemodynamic response indeed occurs at the layers receiving synaptic input.

Functional Response versus Vascular Reactivity to Hypercapnia.

The magnitude and dynamics in BOLD responses are often considered to be affected by different intrinsic vascular responses among regions. Thus, a hypercapnic challenge may be used to separate the dynamics of vascular-driven and neural-driven HRFs. In our study, neural activity–driven BOLD responses (0.6 to 1.5 s onset times) occurred much earlier than hypercapnic responses (2.3 to 3.2 s) (Figs. 1 , iii and iv and 2, iii and iv versus , iii and iv), indicating that the early BOLD response reflects neural activity. The difference in the onset times of BOLD responses induced by neural stimulation and hypercapnia can be explained in terms of underlying mechanisms. Neural activity–driven vascular responses occur first in the capillary bed via pericytes, precapillary sphincters, smooth muscle cells, and endothelial cells and then propagate to upstream (arterial network toward the pial surface) and downstream (venous) vessels (15, 60–64). Thus, neural activity–driven capillary action is the earliest within the vasculature. Meanwhile, vascular action, caused by an increased CO2 level, is likely mediated by an increase in extracellular pH that leads to the relaxation of vascular smooth muscles in the arteries (65, 66) but does not involve the endothelium (67) or neurons (68). Therefore, hypercapnia-induced vascular changes occur first in the arteries, followed by the capillaries and veins (69). The early BOLD response to vascular changes is likely related to the arrival of blood to the capillary bed in a downstream process. In summary, the onset time of the BOLD response to hypercapnia is highly sensitive to arterial blood delivery, whereas that to neural stimulation reflects the capillary responses near active neurons.

How the Order of Onset Times Reflects the Sequence of Neural Events.

Onset response differences of a few hundred milliseconds among different layers or regions appear to reflect the sequence of neural activation on the scale of a few milliseconds. Note that the difference of neural latency times, which are generally determined by the time of the first spike, is ∼5 ms for first-order thalamic relay nucleus to S1 (26), ∼8 ms for S1 to M1 (23, 28), and ∼10 ms for S1 to higher-order thalamic nucleus PO (26). The most critical issue is why the hemodynamic onset response is much slower than the neural latency time. The hemodynamic response measured by BOLD fMRI is closely related to local field potentials (70–73), which reflect the activity of synchronized synaptic afferents, not output spiking activity (74). Even though the first spike occurs within a few milliseconds, the synaptic input activity required to induce one spike might not induce a detectable hemodynamic response. The vascular response is expected to monotonically increase with the magnitude of local field potentials (70–73), which is closely related to the number of active neurons in a given volume and time (75). To be detected by fMRI (or other hemodynamic tools), a measurable vascular response that meets the “detectable vascular threshold” is required. Note that the vascular threshold depends on the sensitivity of the BOLD fMRI. The detectable vascular threshold can be reached faster by a higher number of active neurons and more synchronized activity, which effectively increases vasoactive signaling by releasing ions and vasoactive mediators such as adenosine, prostaglandin E2, and nitric oxide (63, 64).

Wherever neural projections occur, synaptic input activity at the downstream sites is likely to be less synchronized than that at the upstream site, causing it to require a longer time to reach the vascular threshold. Perhaps the hemodynamic delay tracks the sequence of the vascular threshold–reaching neural activity by propagating the onset response at every synaptic connection. Further systematic studies are needed to understand the temporal hemodynamics relative to the amount of neural activity in different layers or regions.

Limitations and Potential Pitfalls.

To utilize early dynamics of BOLD fMRI for systems-level fMRI research, some caution should be exercised. In fMRI, each voxel contains all sizes of vessels, including draining veins. In order to detect early capillary-sensitive BOLD responses with high spatiotemporal resolution, the use of ultrahigh fields is preferable. In our studies, anesthesia was used to minimize head motions and to increase experimental time for extensive averaging, consequently reducing the contribution of CC feedback circuits (76, 77). Under wakefulness, both feedforward and feedback networks contribute to BOLD HRFs, potentially complicating the interpretation of onset times. To interpret onset time differences as sequential neural information flow, a prior network knowledge is necessary, which can be obtained from anatomical and/or functional connectivity (24, 78–80). In addition, Granger causality and dynamic causal modeling can be supplemented to determine neural information flow.

Conclusion

In conclusion, the order of early BOLD responses to somatosensory and optogenetic stimulation was consistent with the sequence of neural information flow in a well-understood somatosensory network. Mesoscopic, population-level, synaptic input–based, ultrahigh-resolution BOLD fMRI at an ultrahigh field can thus be used to investigate whole-brain functional information flow. For example, PO can be involved in spinothalamic ascending input and CT descending input. Our data clearly indicate that the BOLD response of PO is due to the CT projection from S1, which is supported by electrophysiology studies (81). As ultrahigh magnetic fields of ≥7 T become increasingly available, fMRI with high spatiotemporal resolution could be used to explore sequential neural events in a large brain region, which could have a major impact on system neuroscience research.

Materials and Methods

All animal experiments were conducted with approval from the Institutional Animal Care and Use Committee of Sungkyunkwan University and under published standards for humane animal care. Naïve C57BL/6 mice were used to determine the timing of functional responses to somatosensory stimulation (n = 8) and hypercapnic challenge (n = 7), and transgenic Thy1-ChR2 mice (n = 6) were used in the optogenetic fMRI experiment. The details of animal preparation, fMRI experiments, and data analysis are provided in the .