Focal infrared neural stimulation with high-field functional MRI: A rapid way to map mesoscale brain connectomes

See allHide authors and affiliations

Science Advances  24 Apr 2019:
Vol. 5, no. 4, eaau7046
DOI: 10.1126/sciadv.aau7046


We have developed a way to map brain-wide networks using focal pulsed infrared neural stimulation in ultrahigh-field magnetic resonance imaging (MRI). The patterns of connections revealed are similar to those of connections previously mapped with anatomical tract tracing methods. These include connections between cortex and subcortical locations and long-range cortico-cortical connections. Studies of local cortical connections reveal columnar-sized laminar activation, consistent with feed-forward and feedback projection signatures. This method is broadly applicable and can be applied to multiple areas of the brain in different species and across different MRI platforms. Systematic point-by-point application of this method may lead to fundamental advances in our understanding of brain connectomes.


A new method for mapping circuits with mesoscale resolution

Mapping connection patterns in the brain is essential for understanding brain networks and its relationship to behavior and disease. There are many approaches to mapping connections in the brain. However, currently, there is no method to systematically map connections in vivo rapidly with millimeter-scale resolution at a brain-wide scale. Anatomical mapping methods based on tracer injections 1 to 5 mm in size are limited to a small number of injection sites (typically three to five tracers) and require 2 to 3 weeks for tracer transport, animal sacrifice, and time-consuming reconstruction. Mapping of resting state connectivity (1), based on low-frequency correlations between brain areas, is commonly used for the study of brain networks; however, the relationship of these correlations to anatomical connectivity is still uncertain. Diffusion methods (such as diffusion tensor imaging) can achieve high spatial resolution but are limited to larger fiber tracks when used at a whole-brain scale. There have been recent developments with electrical stimulation mapping with functional magnetic resonance imaging (fMRI) in animals [e.g., (2, 3)] and in humans [e.g., (4)]. However, issues of current spread, distortion, and signal dropout near the electrode in the magnet make it difficult to achieve focal millimeter-scale stimulation and to map local connections. fMRI in combination with optogenetic stimulation is a novel and focal method of brain-wide mapping (5); however, in animals such as primates [e.g., (6, 7)], this relies on injection of virus, which takes several weeks for expression and is limited to a small number of sites. Here, we have developed an alternative method for mapping brain connections that combines pulsed near-infrared neural stimulation (INS) with ultrahigh-field fMRI.


INS is an optical stimulation method that induces excitatory and inhibitory neuronal responses, as assessed with intrinsic signal optical imaging, calcium imaging, and electrophysiology (811). INS is mediated via rapid heat changes conferred by pulsed near-infrared wavelength light (12). The mechanism by which a transient temperature rise induces neural depolarization is unknown. Shapiro et al. (13) found that infrared light induces temperature changes characterized by a rapid rise followed by a slow decay (100 ms) similar to the heat relaxation time of water and elicits intensity-dependent depolarizing currents, suggesting that neuronal activation may be mediated by changes in membrane capacitance induced by the transient heat change profile (14,15). Other mechanisms that have been proposed include nanoporation, photomechanical effects, and activation via heat-sensitive transient receptor potential channels [for a review, see (16, 17)]. A wavelength of 1875 nm is selected because, in neural tissue, it achieves thermal confinement (is away from absorption minima where the increased penetration depth leads to greater tissue scattering) and achieves a tissue penetration depth of a few hundred micrometers (12); in cortical tissue, this reaches the superficial cortical layer neurons and apical dendrites of deep layer neurons. In general, INS leads to excitatory effects on neurons and neuronal populations and produces intensity-dependent effects. However, high-intensity stimulation of the cerebral cortex can result in suppressive effects due to activation of inhibitory cortical networks [(8, 9); for a review, see (16)]. Thus, in the brain, INS affects both excitatory and inhibitory circuits, the balance of which brings about a resulting net response. This net population response can also be monitored via imaging of hemodynamic response (8, 9, 16). Within certain parameters (dependent on optical fiber size, frequency, and radiant energy), this method can be applied repeatedly without tissue damage or detriment to cortical function or health of the animal [for a review, see (17)]. Note that, with different stimulation parameters (lower energies and longer duration), INS can be used to directly inhibit neural response (10, 18, 19).

One advantage of INS is that light delivered through an optical fiber produces a spatially restricted stimulation site (20, 21). On the basis of light scattering modeling studies (21), direct measurement in neural tissue (22), and histological evidence (23), 1875-nm light distribution in tissue is similar in size to the optical fiber diameter and reaches several hundred micrometers in tissue depth. This focal effect is consistent with in vivo observations of INS stimulation of fascicles within peripheral nerves (24) and modulation of auditory thresholds in cat cochlea (11). For this study, we have demonstrated, by mapping hemodynamic signals using intrinsic signal imaging, that INS delivered through a 200–μm optical fiber can selectively target and modulate single ocular dominance columns in monkey visual cortex to achieve enhanced response of matching eye columns and relative suppression of nonmatching eye columns (9). This suggests that such an approach can be used to stimulate single submillimeter domains and their associated networks. Here, we examine whether INS-induced hemodynamics signals can be used in ultrahigh-field MRI to map mesoscale resolution networks at both global and local scales.


Global scale mapping

We chose to study the connectivity of cat primary visual cortex because its cortico-cortical circuitry is well understood and is based on extensive neuroanatomical tracing studies. We targeted the border of area 17 and 18 [the vertical meridian near the horizontal meridian (25, 26)] in hopes of visualizing both ipsilateral and contralateral connections. Since INS applied to the surface of the cortex reaches a penetration depth of a few hundred micrometers, it is expected to activate cortico-cortical projection neurons within the superficial layers and apical dendrites of subcortical projection neurons in the deep layers. Thus, based on known anatomy (Fig. 1B and fig. S1), we expected to see activation in the thalamus [lateral geniculate nucleus (LGN) and lateral posterior nucleus and pulvinar (LP-pulvinar complex)], primary cortical targets (areas 17, 18, and 19), and higher cortical areas [areas 20 and 21, AMLS (anterior medial lateral suprasylvian areas), and PMLS (posterior medial lateral suprasylvian areas)]. We also wondered whether polysynaptic circuits would be revealed.

Fig. 1 Global connections revealed by INS stimulation in cat visual cortex.

(A) Stimulation paradigm. Each block of INS comprises 20 trials of pulse trains at a single energy. Train duration, 0.5 s; intertrain interval, 35.5 s; pulses, 250 μs; 200 Hz; 0.1 to 1.0 J/cm2 per pulse. (B) Connections of cat area 17/18. Red arrows indicate cat visual projections from area 17 [adapted from (49)]. The color bar shows T statistics for significant voxels in (C) and (D), thresholded at P < 0.01. (C to E) Data from cat #4. (C) Thalamic activations. Left: Stimulation of area 17/18 at 0.3 J/cm2 [red asterisk in (D)] reveals ipsilateral LGN activation. Right: At 0.7 J/cm2, stimulation produces larger ipsilateral activation in LGN and pulvinar (LGN/Pul). In addition, strong activations are seen in contralateral pulvinar and lateral posterior nucleus (LP). (D) Cortical activations in three slices (dorsal, intermediate, and ventral). The red asterisk denotes the stimulation of area 17/18 at 0.7 J/cm2. Significant voxels seen in areas 17, 18, and 19; AMLS; anterior ectosylvian area (AES); and cingulate visual area (CVA). (E) Intensity dependence (average of 20 trials; error bars, SEM). Top: Time courses (left) and amplitudes (right) at the laser tip site in (D). The yellow bar indicates INS onset. Bottom: Response amplitudes at connected ipsilateral 18 and contralateral 17/18 sites in (D).

In these experiments, a 200-μm optical fiber delivered INS pulse trains (250 μs pulse width, 200 Hz, 0.5 s duration, 36 s per trial, and 20 trials per intensity; see Fig. 1A) to the cortical surface. Nondamaging radiant energies (0.1, 0.3, 0.5, 0.7, or 1.0 J/cm2 per pulse) were used. Blocks of different radiant energies were presented in a randomly interleaved fashion. All images were acquired at 1- or 1.5-mm voxel size using a single-shot echo-planar imaging (EPI) on a 7-T MRI system. Responding voxels were identified with a P value threshold at a false discovery rate (FDR) of 10 or 20%.

Activation sites. The result of INS stimulation in the visual cortex of the right hemisphere at 0.7 J/cm2 is shown in Fig. 1 (C and D) (all slices are shown in Fig. 2). On the basis of previously published studies of the cat visual system and a digital atlas (2529), the stimulation site in Fig. 1D (red asterisk, 0.7 J/cm2) corresponds to a visual cortical location near the area 17/18 border (25, 26) that extends into both areas 17 and 18 (Fig. 1C). This activation spanned several slices (area 17/18). Activations from this site will be interpreted as arising from either area 17 or area 18 and will be referred to as area 17/18.

Fig. 2 All slices of voxels activated by laser stimulation of the area 17/18 border (0.7 J/cm2).

Significant voxels (P < 0.01, FDR of 20%) overlie locations with known anatomical connections to laser stimulation site in area 17/18. These include LGN; LP/pulvinar; areas 17, 18, 19, 20, and 21; and sylvian areas (AMLS and PMLS). Polysynaptic areas include contralateral LP/pulvinar, AES, CVA, and areas 4, 5, and 7 (see the main text). Note that LGN, which has very strong connections with area 17/18, contains some of the most significant voxels. The color bar shows T statistics for significant voxels. See Fig. 1A for the expanded forms of the abbreviations.

We observed clusters of activation in locations known to be directly connected to area 17/18, which appeared focal and specific. In the thalamus, this included the ipsilateral LGN (Fig. 1, C and D), where strong activation was observed, as well as ipsilateral LP and pulvinar (Fig. 1C, right). Laser stimulation delivered at a lower intensity (0.3 J/cm2) produced activation primarily in the ipsilateral LGN (Fig. 1C, left), suggesting an intensity dependence.

In the cortex, activation was observed in ipsilateral areas 18, 19, 20, and 21, as well as in AMLS and PMLS (Figs. 1D and 2). Clusters of voxels were also observed in contralateral areas 18, 19, 20, and 21 (Figs. 1, C and D, and 2). We also noted that the callosal projection to contralateral area 18 is not mirror symmetric, as previously described (30).

We then examined whether INS stimulation of area 17/18 would lead to activation at sites that can only be polysynaptically connected. Unexpectedly, stimulation at 0.7 J/cm2 produced activation in the contralateral LP and pulvinar (Fig. 1C), a result that can only come about from activation of callosally connected cortex (such as contralateral area 18 or 19; Figs. 1C and 2) (31). This latter finding suggests that focal cortical INS stimulation can lead to polysynaptic (ipsilateral visual cortex–to–contralateral visual cortex–to–contralateral thalamus) activation. Note that different thresholds can reveal slightly different activations. For example, a threshold of P < 0.01 (FDR of 20%) reveals a large thalamic activation (Fig. 2, fourth row, two rightmost slices), whereas a higher threshold of P < 0.0025 (FDR of 10%) of these same data reveals two distinct nuclei (LP and pulvinar) within this region (see fig. S2, fourth row, two rightmost slices).

Other cortical sites of indirect activation (Figs. 1, C and D, and 2) included (i) multisensory area anterior ectosylvian area (AES), which receives visual inputs via thalamic LP and cortical lateral sylvian areas (Fig. 1D) (32, 33); (ii) cingulate visual area (CVA; Fig. 1D, left) and splenial visual area (Fig. 2), which receives input from 17/18 via the LP/pulvinar complex (34, 35); (iii) area 5 and area 7 (Fig. 2), areas with strong inputs from pulvinar (36); and (iv) part of motor area 4 (Fig. 2), involved with visually guided motor movements (37, 38). In summary, the overall set of sites appears focal and specific and consistent with known anatomical connectivity from a focal site in areas 17 and 18 (see fig. S1).

Intensity dependence. We had previously shown, using intrinsic signal optical imaging, that INS induced a consistent hemodynamic response at the laser tip characterized by intensity dependence and a signal time course with faster rise times than sensory-induced intrinsic signals (8, 9). To further examine the reliability of INS-induced MR signal, we measured the intensity dependence of response at the stimulation site and connected sites.

As shown in Fig. 1E (top panels), intensity dependence was evident in the response amplitude at the laser tip site. Stimulation at the lowest intensity (0.1 J/cm2) produced the weakest response, whereas increasing intensity led to an increasing magnitude of response (response magnitude of the voxel at the laser tip; time courses are plotted in the top left panel and amplitudes are plotted in top right panel). In addition, note that, consistent with the optical imaging studies, time courses exhibited faster rise times than sensory-induced MR signals (Fig. 1E, top left).

At connected sites, in both the ipsilateral and the contralateral hemisphere, we also observed a general increase of response amplitude with increasing stimulation intensity (Fig. 1E, bottom panels). Consistent with activation due to synaptic transmission, the amplitudes of response activation at these sites were weaker and exhibited greater variability than that at the laser tip (compare error bars in Fig. 1E). Some of these sites increased in area at higher activation intensities (e.g., compare LGN activation in both panels of Fig. 1C). The intensity dependence of distant site activation to the INS stimulation intensity suggests that these activations were connected to the laser stimulation site.

Local scale mapping

We also examined whether INS could be used to reveal local intra-areal and inter-areal connections at high resolution. As previously shown (3941), squirrel monkey primary somatosensory cortex (SI; comprising areas 3a, 3b, 1, and 2) contains a somatotopic map of the digits (D1 lateral to D5 medial; Fig. 3A). The digit networks are organized in two orthogonal axes: an intra-areal axis that links digits within a single area (Fig. 3E, red ovals in the mediolateral axis), proposed to underlie sensation during manual power grasp, and an inter-areal axis that links single-digit information across areas (Fig. 3E, red ovals in the anteroposterior axis), proposed to underlie sensation during individuated digit movements. SI is also characterized by specific inter-areal feed-forward (middle layer labeling) and feedback (superficial and deep) patterns of connectivity (Fig. 3F) (42, 43).

Fig. 3 Local connections revealed by INS stimulation in squirrel monkey somatosensory cortex.

(A) Somatotopy in areas 3b and 1 (38). (B) Stimulation paradigm. ISI, inter-stimulus interval. (C) A 400-μm optic fiber on top of an artificial dura in an optical chamber over SI. Surface coil, 3 cm. Dark vessel, central sulcus. (D) Two planes of MRI imaging, with the chamber seen on top of the cortex. Scale bar, 1 cm. (E) Schematic of digit connectivity in SI (33). (F) Schematic of feed-forward and feedback inter-areal connectivity in SI (35). (G) Cortical responses to stimulation at D3. Slice in tangential plane. The red dot indicates the optical fiber tip in area 3b. See Fig. 4 (A and C) for determination of areal borders (dashed white lines). CS, central sulcus; M, medial; P, posterior. Scale bar, 5 mm (monkey #2). (H) Cortical response to stimulation of area 2 (red arrowhead, fiber tip; dark spot, signal dropout). Slice in orthogonal plane. Green arrowheads indicate feed-forward activations. White and blue arrowheads denote feedback. See Fig. 4 (B and D) for determination of areal boundaries. The white rectangle indicates the silicone plug in the chamber. D, dorsal; LS, lateral sulcus. Scale bar, 5 mm. The color bar shows T statistics for significant voxels in (G) and (H) (monkey #3). (Photo credit: Gang Chen and Anna Wang Roe.)

To examine whether INS could reveal these local patterns of connectivity, we first determined digit locations in areas 3b and 1 using intrinsic signal optical imaging and vibrotactile stimulation of digit tips (two cases shown in Fig. 4, A and B; note that digit representations, colored circles, are ~1 mm in size). We then applied INS stimulation (12 pulse trains per trial; pulse train duration, 0.5 s; intertrain interval, 2 s; and 100 pulses per pulse train; Fig. 3B) through a 400-μm INS optical fiber targeting a selected digit site (Fig. 3C). Using a custom 3-cm surface coil (Fig. 3C) in a 9.4-T MRI scanner and a contrast agent [monocrystalline iron oxide nanoparticle (MION)], we acquired imaging slices in planes tangential (Fig. 3D, red) or orthogonal (Fig. 3D, green) to the cortical surface.

Fig. 4 Determination of areal boundaries via optical imaging

(A) Monkey #2 (Fig. 3G, tangential case) and (B) monkey #3 (Fig. 3H, orthogonal case). Left: Optical chamber. Blue boxes denote imaged fields of view. Right: Optical images to stimulation of digit tips [cf. (42)]. White arrows in D3 of (A) correspond to those in (C) and (E). Field of view for D1 and D2 images is shifted relative to that for D3 and D4. (C and D) Summary of digit tip topography from (A) and (B). (C) Arrows indicate corresponding locations with (E). The red dot indicates the laser stimulation site. White dotted lines denote the borders between areas. (D) Colored arrowheads indicate corresponding locations with (F). The red arrowhead denotes the laser stimulation site. The horizontal purple line indicates the approximate location of slice in (F). Green arrowheads denote the feed-forward projections in (F). White (area 3b) and blue (area 1) arrowheads indicate feedback projections in (F). (A and C) Monkey #2. (B and D) Monkey #3. (E and F) Compare with (C) and (D) for alignment of optical imaging and fMRI maps. (E) Same as Fig. 3G. (F) Same as Fig. 3H. These figures are shown to easily see the correspondence between optical and MRI images. (Photo credit: Robert M. Friedman and Anna Wang Roe.)

Activation patterns. Figure 3G (tangential slice) illustrates an example in which the fiber tip was targeted to a D3 digit location [red dot; see Fig. 4, A and C, for digit maps; green and white arrows point to corresponding locations in Figs. 3G and 4C; see (44) for optical imaging/fMRI alignment methods]. INS stimulation (0.29 J/cm2) produced a ~2-mm activation at the laser tip site as well as 200- to 300-μm-wide focal activations away from the stimulation focus (white and green arrows). Activation sites were consistent in location with matching D3 digit representation in areas 3a and 1 (green arrows) and with anatomical interdigit connections within area 3b (small white arrows), similar to the two orthogonal axes previously described (Fig. 3E). Since the fiber tip targeted a digit tip location, the activations fell largely near the area 3a/3b border and not near the area 3b/1 border, topographically consistent with tip-to-tip representation at the area 3a/3b border and palm-to-palm representation at the area 3b/1 border (41).

In a separate case (Fig. 3H; see Fig. 4, B and D, for digit maps; colored arrowheads indicate corresponding areal locations), the laser fiber tip was placed in area 2 (red arrowhead) and slices were scanned in the orthogonal plane. Near the INS tip site, stimulation (0.50 J/cm2) resulted in a ~2-mm activation as well as anterior and posterior focal activations (colored arrowheads). On the basis of our optical imaging of digit maps (Fig. 4, B and D), we determined the locations that fell in areas M1, 3a, 3b, and 1 anteriorly and in area 5 posteriorly. We noted that the active voxels in areas 3a and M1 fell in the middle cortical layers (green arrowheads) and that those in areas 3b and 1 exhibited a bilaminar distribution (white and blue arrowheads, respectively). In area 5, an area known to receive feed-forward inputs from area 2, the laminar position of activation was less clear as it fell on the crown of the lateral sulcus where large vessels reside. Overall, these patterns are consistent with previously described feed-forward (to areas 3a and M1) and feedback (to areas 3b and 1) patterns of connectivity (Fig. 3F) (42). Activation sites were similar in size to anatomically labeled patches (~200 to 300 μm).

Intensity dependence and reliability. We also examined the intensity dependence of the MR signal with INS laser stimulation. As shown in Fig. 5A, MR signal was consistent across runs in terms of both timing and magnitude (see error bars). In general, the greater the stimulation intensity, the stronger the MR signal change (see color legend in Fig. 5A; amplitudes are plotted in Fig. 5B, one dot per trial). It is unclear why, at the 0.39 J/cm2 energy level (blue dots), the response reached a lower than expected response level. As shown in Fig. 5 (C and D, 0.5 and 0.29 J/cm2, respectively), spatial aspects of activation were also intensity dependent. The inter-areal connection patterns were more prominent at higher stimulation intensities (0.5 J/cm2) and weaker when intensity was reduced by roughly half (0.29 J/cm2) (compare Fig. 5, C and D), supporting the hypothesis that distant activations were related to INS stimulation.

Fig. 5 Intensity dependence.

(A) Time course of MR signal change at different intensities of INS stimulation (average of 13 trials). The same case as Figs. 3G and 4 (A, C, and E) (monkey #2). (B) Each dot denotes the peak amplitude of MR signal of individual trials (n = 13). Error bars, SEM. The dashed line represents the regression line. (C and D) Local activation at high [(C) 0.50 J/cm2] and low [(D) 0.29 J/cm2] intensities. The same case as in Figs. 3H and 4 (B and D) (monkey #3). Stimulation in area 2 (optical fiber, red arrow). Third panel: Same as in Fig. 3H. Adjacent 1.5-mm-thick slices from lateral (left) to medial (right). Central sulcus [seen in the two rightmost panels of (C)] aided in demarcation of area 3b (white arrowheads) and area 1 (blue arrowheads). (D) Same as (C) except for lower INS intensity. Scale bar, 5 mm. The color bar shows T statistics for significant voxels in (C) and (D).

To further test the reliability of our data, we compared the average of even versus odd trials. We find that the half-trial images are largely similar (Fig. 6). Since half as many trials will have a lower signal-to-noise ratio, we used a P value threshold of P < 0.0005 (for half-trial averages) rather than P < 0.0001 (for all trial images). As shown in Fig. 6, comparison of the half-trials reveals that the potential feedback activations (two white and two blue arrowheads) are (almost) all observed in both the even and odd trials. The potential feed-forward activations (green arrowheads) are observed in the even and odd trials in area 5 but appear only in the odd trials in areas 3a and M1. We suspect that this may be due to an issue of lower signal-to-noise ratio and that collection of a greater number of trials would improve the stability of these signals. Although the data are suggestive, the strength of this finding needs to be further examined in experiments with a larger number of trials.

Fig. 6 Reliability of activation.

Panels A, B, C and D represent adjacent slices. Top row: Original panels (Fig. 5 C, n = 13 trials) showing feed-forward (green arrowheads) and feedback (white and blue arrowheads) connections from point stimulation in area 2 (red arrow). Middle row: Sum of even trials (n = 6). Bottom row: Sum of odd trials (n = 7). Since half as many trials will have lower statistical power, we used P < 0.0005 for half-trial images (middle and bottom rows) instead of the P < 0.0001 used for all trials (top row). Comparison of the half-trials reveals that the images are largely similar (from monkey #3).


In summary, we demonstrate the feasibility of using focal INS stimulation in the MRI for in vivo mapping of brain networks. Activation locations were consistent with known long-range connections in cat visual cortex and local short-range connections in squirrel monkey somatosensory cortex. The specificity and focal nature of activated target sites make it unlikely that these patterns were due to random activations and suggest a true reflection of anatomical connectivity. Furthermore, these signals can be mapped under quite different conditions (two different sensory cortical areas, two different species, and two different ultrahigh-field MRI systems), indicating generality of methodology. Compared with anatomical mapping, INS-fMRI offers certain advantages. Connections can be studied in vivo, thereby reducing the number of animals needed. fMRI maps are rendered in three dimensions directly on anatomical brain scans (eliminating time-consuming reconstruction). Furthermore, it is rapid (can be acquired within a single 1- to 2-hour fMRI session). These advantages permit the study of multiple networks within a single animal, further opening possibilities of correlating multiple types of datasets (e.g., electrophysiology and behavioral data) from an individual animal. These results demonstrate that the use of ultrahigh-field fMRI provides sufficient signal-to-noise ratio for high spatial resolution mapping, allowing the study of brain-wide networks at mesoscale.


Cat protocol

Animal preparation. All procedures were in compliance with and approved by the Institutional Animal Care and Use Committee of Zhejiang University and followed the guidelines of the National Institutes of Health Guide for the Care and Use of Laboratory Animals. Four cats were anesthetized with ketamine hydrochloride (10 mg/kg)/atropine (0.03 mg/kg) and maintained with propofol (5 mg/kg per hour; induction, 5 mg/kg) or sufentanil (2 to 4 μg/kg per hour CRI (continuous rate infusion); induction, 3 μg/kg) supplemented with up to 1% isoflurane. Animals were intubated and artificially ventilated. End-tidal CO2, respiration rate, SpO2, heart rate, and electrocardiogram were monitored and maintained. Temperature (37.5° to 38.5°C) was monitored and maintained via a circulating water blanket. Laser stimulation (Laserglow Technologies, Toronto, Canada) was delivered via a 200-μm-diameter optic fiber through a 1-mm burr hole in the skull overlying visual cortex.

Stimulus protocol. For all INS experiments, stimuli consisted of a train of laser pulses (wavelength, 1875 ± 10 nm; radiant energy, 0.1 to 1.0 J/cm2 per pulse). Laser pulses were calibrated with a power meter and transmitted through a 200-μm-diameter optical fiber positioned normal to the brain surface. Each train consisted of 0.5-s duration burst (250 μs per pulse at 200 Hz) presented once every 36 s. In each scanning block, 20 pulse trains (trials) at a single intensity were delivered. Blocks of different intensities were randomly interleaved. For analysis, the first half (18 s) was used as the baseline period, and the second half was used as the stimulus response period (see Fig. 1A).

MRI methods. Cat data were acquired on a 7-T research scanner (Siemens Healthcare, Erlangen, Germany) equipped with a whole-body gradient set (70 mT/m and 200 T/m per second), with a single loop coil (RAPID MR International, Columbus, OH, USA) for transmission and signal reception. Structural images were obtained at a voxel size of 0.5-mm isotropic cube [echo time (TE), 20 ms; repetition time (TR), 3000 ms; matrix size, 128 × 128; flip angle, 120°; 30 slices]. Functional images were obtained with a single-shot EPI sequence at a voxel size of 1.0-mm isotropic cube (20 slices) or 1.5-mm isotropic cube (10 slices) (TE, 22 ms; TR, 2000 ms; matrix size, 64 × 50; flip angle, 90°).

MR data preprocessing. DICOM files were converted to AFNI (Analysis of Functional NeuroImages) data format with AFNI Dimon (45). EPI images were then processed with slice-timing correction, motion correction, and removal of baseline shift with a third-order polynomial function. EPI (1.0- and 1.5-mm isotropic cube) images were aligned to structural images (0.5-mm isotropic cube), with and 3dQwarp in AFNI. All EPI images were resampled to 0.5-mm resolution to match the grid space of the structural images. The transformation matrices acquired during EPI image registration were then used for presenting statistical test results on structural images (in Figs. 1 and 2 and fig. S2).

Detection of significant voxels. The detection of significant voxels followed the general linear model approach for fMRI analysis. The time course of each voxel was modeled as a weighted sum of several predictor variables (stimulus and motion predictors) plus an error term. The stimulus predictor was the convolution of INS onset with a slow rising function with a 10- to 20-s peak (46). The estimated weight of the stimulus predictor was subjected to a t test to examine its significance. Significant voxels were identified with a P value threshold at an FDR of 10 or 20% (Benjamini-Hochberg method).

Time course analysis. Individual voxel time courses were extracted from EPI data with AFNI 3dmaskdump. Time courses were then averaged over 20 repetitions and plotted. Peak amplitudes were estimated as the maximum percentage changes within 18 s after laser onset.

Atlas. Functional annotations of the response sites (in Figs. 1 and 2 and fig. S2) were determined according to the atlas called Catlas (29). To use this atlas, we aligned the atlas T1 images to the structural images of cats. We then applied the same transformation to the atlas segmentation to ensure good alignment with the anatomy. The coregistration was done with the AFNI program and 3dQwarp. The software ITK-SNAP was also used to visualize the atlas and the activation sites.

Monkey protocol

Animal preparation. All procedures were in compliance with and approved by the Institutional Animal Care and Use Committee of Vanderbilt University and followed the guidelines of the National Institutes of Health Guide for the Care and Use of Laboratory Animals. Three squirrel monkeys were anesthetized with ketamine hydrochloride (10 mg/kg)/atropine (0.03 mg/kg) and maintained with isoflurane (0.5 to 3%) delivered in a 70:30 O2/N2O mixture. Animals were intubated and artificially ventilated. End-tidal CO2, respiration rate, SpO2, heart rate, electrocardiogram, and rectal temperature were monitored by MR compatible capnography, fiber-optic pulse oximetry, and small animal monitoring system. The end-tidal CO2 was kept around 4% by adjusting the rate and volume of the ventilator, and the rectal temperature of 37.5° to 38.5°C was maintained by an animal warming system. Under general anesthesia (1 to 2% isoflurane), two monkeys were implanted with custom-made imaging chambers, with a transparent artificial dura over the primary somatosensory cortex, secured by ceramic screws. After placement in an MRI-compatible stereotaxic tube, laser stimulation (Lockheed Martin Aculight, Bothell, WA, USA) was delivered via a 400-μm-diameter optic fiber through the imaging chamber.

Stimulus protocol. Laser stimuli consisted of a train of laser pulses (wavelength, 1875 ± 10 nm; radiant energy, 0.1 to 1.0 J/cm2 per pulse, calibrated with a power meter) via a 400-μm-diameter optical fiber positioned normal to the brain surface. Each train consisted of 0.5-s duration bursts of 250 μs per pulse at 200 Hz. In a single trial, 12 pulse trains were repeated once every 2.5 s for 30 s (interstimulus interval between pulse trains was 2.0 s), followed by 30-s blank. Blocks of 10 to 15 trials (300 imaging volumes) were acquired within each scanning block of the same laser energy.

MRI methods. All MRI scans were performed on a 9.4-T (Varian Medical Systems, Palo Alto, CA, USA) horizontal bore magnet with an inner bore size of 21 cm and an actively shielded gradient (400 mT/m and 3000 T/m per second). A custom 3-cm-diameter transmit-receive surface coil was used in all MR scans, with injection of MIONs at 10 to 20 mg of Fe per kilogram of animal weight. The signal of all functional data with MION injection was inverted to facilitate comparison with BOLD data. T2-weighted oblique structural images at 68 μm by 68 μm in-plane resolution were acquired to identify venous structure on the cortical surface for locating the central sulcus and lateral sulcus and for coregistration of fMRI maps across imaging sessions. fMRI data were acquired from the same slices using a two-shot or four-shot EPI with a TE of 10 ms and a TR of 750 ms at voxel sizes of 270 μm by 270 μm by 2000 μm (tangential slices) or 270 μm by 270 μm by 1500 μm (orthogonal slices).

MR data processing. Functional EPI data were preprocessed for motion correction and then analyzed in MATLAB (Mathworks, Natick MA, USA) codes [for details, see (44)]. EPI data were collected as a 128 × 128 matrix. Raw EPIs were interpolated to a 512 × 512 matrix and overlaid on anatomical images. Time courses were drift-corrected using a linear model and temporally smoothed with a low-pass filter (cutoff frequency set above the second harmonic of task function). The correlation of each functional time course to a reference waveform was calculated, and functional maps were generated by identifying regions of clustered voxels whose correlation with the reference waveform was significant to P < 10−4. Choice of this threshold was guided by area, amplitude, and centroid of activation determined via optical imaging and electrophysiology in the same animals. Activations were overlaid on corresponding anatomical images and then compared with known topography.

Optical imaging/fMRI image registration. Somatosensory areal borders in squirrel monkeys were estimated by collecting optical imaging and by distinguishing receptive field preferences by electrophysiological single and multi-unit recording. Digit maps were obtained by vibrotactile stimulation of each digit. For coregistration, we first identified corresponding anatomical and blood vessel landmarks in surface blood vessel images and structural images, such as visible surface vessels and transcortical veins. These coordinates were then put into a point-based registration algorithm [implemented in MATLAB; for details, see (47, 48)]. We have used these mapping methods in a number of publications [e.g., (39, 40, 44)].


Supplementary material for this article is available at

Fig. S1. Cat anatomy.

Fig. S2. All slices of significant voxels controlled at an FDR of 10% under laser stimulation of the area 17/18 border (0.7 J/cm2).

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank John C. Gore and Li M. Chen for discussion and Fang Wang, Chaohui Tang, and Yimin Shen for technical assistance. Funding: This research was conducted at Zhejiang University and Vanderbilt University and was supported by National Natural Science Foundation grants 81430010 and 31627802 (to A.W.R.), 81701774 and 61771423 (to X.T.Z.), and 31471052 (to G.C.); National Hi-Tech Research and Development Program grant 2015AA020515 (to A.W.R.); Zhejiang Provincial Natural Science Foundation of China grant LR15C090001 (to G.C.); Fundamental Research Funds for the Central Universities grants 2015QN81007 (to G.C.) and 2016QN81018 (to X.T.Z.); NIH grants NS044375 and NS093998 (to A.W.R. and R.M.F.); the Vanderbilt Imaging Institute; and NIH grant R01 NS052407-01 (to D.J.). Author contributions: A.G.X. designed and conducted cat experiments, analyzed cat data, developed analysis pipeline, created figures, and wrote the manuscript. M.Q. developed cat methodologies and conducted cat experiments. F.T. conducted monkey data analysis. B.X. conducted additional analysis of monkey data. R.M.F. collected and analyzed monkey optical imaging data and helped in the writing of the manuscript. J.W. helped develop and implement cat data analysis methods. X.S. conducted cat experiments and helped develop cat laser-fMRI procedures. Y.S. is a Siemens on-site scientist for the 7-T facility and provided access to key pulse sequences. M.M.C. is an INS laser expert and played a critical role in calibration, design of laser pulse train delivery, and evaluating the safety of laser energies. J.M.C. provided laser expertise and helped with monkey data collection. E.D.J. and A.M.-J. provided Aculite laser and were involved in the early phases of this study. X.Z. is an MR physicist on our 7-T team; he designed the RF coil, optimized sequences, and conducted MR thermometry measurement. G.C. designed and conducted monkey experiments, analyzed monkey data, created figures, and wrote the manuscript. A.W.R. designed studies, obtained funding, conducted some of the cat and monkey fMRI and optical imaging experiments, supervised all phases of data collection and analysis, created figures, and wrote the manuscript. All authors contributed to the discussion of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article