Behaviors and brain disorders involve neural circuits that are widely distributed in the brain. The ability to map the functional connectivity of distributed circuits, and to assess how this connectivity evolves over time, will be facilitated by methods for characterizing the network impact of activating a specific subcircuit, cell type, or projection pathway. We describe here an approach using high-resolution blood oxygenation level-dependent (BOLD) functional MRI (fMRI) of the awake mouse brain-to measure the distributed BOLD response evoked by optical activation of a local, defined cell class expressing the light-gated ion channel channelrhodopsin-2 (ChR2). The utility of this opto-fMRI approach was explored by identifying known cortical and subcortical targets of pyramidal cells of the primary somatosensory cortex (SI) and by analyzing how the set of regions recruited by optogenetically driven SI activity differs between the awake and anesthetized states. Results showed positive BOLD responses in a distributed network that included secondary somatosensory cortex (SII), primary motor cortex (MI), caudoputamen (CP), and contralateral SI (c-SI). Measures in awake compared with anesthetized mice (0.7% isoflurane) showed significantly increased BOLD response in the local region (SI) and indirectly stimulated regions (SII, MI, CP, and c-SI), as well as increased BOLD signal temporal correlations between pairs of regions. These collective results suggest opto-fMRI can provide a controlled means for characterizing the distributed network downstream of a defined cell class in the awake brain. Opto-fMRI may find use in examining causal links between defined circuit elements in diverse behaviors and pathologies.
Functional MRI (fMRI) based on the blood oxygenation level-dependent (BOLD) signal (Kwong et al. 1992; Ogawa et al. 1990, 1992) is widely used to indirectly measure neural activity in distributed brain networks in humans and nonhuman primates. A recent study has shown that optogenetic strategies, using activation of channelrhodopsin-2 (ChR2)-expressing neurons (Boyden et al. 2005), can be used with high-field MRI imaging to evoke BOLD signals in the anesthetized rodent (Lee et al. 2010). We present data obtained using optogenetic activation of a specific cell class in the awake animal, mapping and characterizing the distributed network responses that result. These opto-fMRI studies collectively open up a wide array of opportunities for exploring the relation between BOLD responses and neural activity and to use rodent models to assess the impact of causal manipulations on distributed brain networks.
The recent paper of Lee et al. (2010) illustrates the potential of optical neural control and fMRI in the anesthetized state, showing a BOLD response to light-driven epochs of stimulation in MI in anesthetized rats transfected with ChR2 in Ca2+/calmodulin-dependent protein kinase IIα (CaMKIIα)-expressing neurons. Like prior analyses of sensory-evoked responses in the human (Boynton et al. 1996), the measured BOLD response began a few seconds after stimulation and decreased ∼6 s after stimulation offset. Lee et al. (2010) also showed downstream effects of optical stimulation by showing a robust BOLD response in the thalamus during local MI stimulation.
These results raise the possibility that opto-fMRI could serve as a tool to explore properties of distributed brain networks—for example, enabling researchers to focus their electrophysiological or anatomical experiments on opto-fMRI-identified sets of brain regions or to examine how large-scale brain network properties are affected by molecular genetic manipulations. By providing a bridge between the increasingly widespread use of optical methods to probe causal neural circuit functions in animals and fMRI, which in humans is typically performed in the awake state, opto-fMRI in the awake mouse may play an important role in the translation of neural circuit insights derived from animal experimentation toward basic and clinical neuroscience in the human.
In this paper, we describe an initial exploration of an awake mouse model of optogenetic stimulation performed during fMRI, thus enabling measurement of large-scale circuit dynamics in a behaviorally relevant state. Using awake mice expressing ChR2 in pyramidal cells in the primary somatosensory (SI) barrel cortex as a model system, we show the detection, in individual mice, of cortical and striatal structures known to be connected to SI (Aronoff et al. 2010; Carvell and Simons 1986, 1987; Chakrabarti and Alloway 2006; Diamond et al. 2008; Ferezou et al. 2007; Megevand et al. 2008; White and DeAmicis 1977). We also show the use of opto-fMRI to assess how the administration of the common anesthetic isoflurane modulates the set of recruited structures downstream of a cell type, as well as the connectivity between them.
Animals and surgery
All procedures were conducted in accordance with National Institutes of Health guidelines and with the approval of the MIT Committee on Animal Care. Eleven wild-type mice (C57BL/6, purchased from Charles River) and 11 transgenic ChR2 mice [line 18, stock 007612, strain B6.Cg-Tg (Thy1-COP4/EYFP)18Gfng/J from Jackson Labs, Bar Harbor, ME; bred in-house with wild-type mice] were used. A summary of animals used for the different experiments throughout this study is presented in Supplemental Table S1.1
For viral injection, eight wild-type mice were first anesthetized with isoflurane (∼1–2% mixed with oxygen), craniotomized (∼0.5 mm wide), and injected with lentivirus encoding for ChR2-GFP under the CaMKII promoter [FCK-ChR2-GFP, 1 μl, as used in Han et al. (2009), over a 30-min period, into the left SI barrel field]. Injections were performed using an injection pump (Quintessential Stereotaxic Injector, Stoelting Co., Wood Dale, IL) driving a 10 μl Hamilton syringe connected to a glass micropipette (100 μm tip) via polyethylene tubing. The system was filled with mineral oil. The coordinates of viral injection relative to bregma were as follows: 1.0 mm posterior, 3 mm lateral, 0.7 mm ventral.
For opto-fMRI, animals were surgically implanted with a headpost (Fig. 1A) atop their skulls using dental cement (CandB-Metabond, Parkell). The headpost (weighing ∼250 mg) and the headpost holder were custom made from Accura 55 plastic, an ABS-like plastic (Fig. 1A). For virally injected animals, the headpost surgery was done at least 3 wk after the surgery for viral injection. For Thy1-ChR2 mice, the kind of craniotomy performed depended on the experiment: in four Thy1-ChR2 mice used for controls for examining negative signal (Supplemental Fig. S2), a small craniotomy (∼0.5 mm) was also drilled at this time through the skull over SI to allow positioning of the optical fiber above the cortical surface; in five other Thy1-ChR2 mice used for near-microscopic resolution opto-fMRI (Fig. 5) and for three control experiments (a blood pressure control, a yellow light control, and a paw stimulation control), instead of a craniotomy over SI, the skull was thinned to one half of its thickness to reduce distortion even further for high-resolution imaging (Fig. 5, A and B). At the end of each surgery, exposed skull (except for that which was thinned) was covered with a thin and uniform layer of dental cement to minimize echo-planar image (EPI) distortion caused by susceptibility mismatch in fMRI. The reason for this final dental cement step was to reduce distortion: Supplemental Fig. S1B shows a pair of images (each maximum-intensity projected across 4 coronal slices) from a mouse, which shows the effect of a uniform layer of Metabond on the EPI distortion artifact. Specifically, this headposted wild-type mouse was imaged [spin-echo EPI (SE-EPI)] without any Metabond (Supplemental Fig. S1Bi) and reimaged after applying a thin and uniform layer of Metabond (Supplemental Fig. S1Bii), showing in the latter case a reduction in artifact.
Awake and anesthetized opto-fMRI of mouse: experimental setup and protocol
Opto-fMRI data were acquired on a 9.4-T (Bruker BioSpin MRI, Ettlingen, Germany), 20 cm ID, horizontal bore magnet. Custom-built radio frequency (RF) transmit-receive surface coils, specifically designed for mice, were used for imaging (Supplemental Fig. S1A). Functional data were acquired using an SE-EPI in the coronal orientation; we compared this protocol to gradient-echo EPI (GRE-EPI) in Supplemental Fig. S1C, in an acute (0.5% isoflurane anesthetized) experiment on a wild-type mouse, and found GRE-EPI to be more susceptible to distortion artifacts than SE-EPI.
Opto-fMRI experiments were performed on eight virally injected wild-type mice, nine Thy1-ChR2 mice, and two wild-type mice. Of the eight virally injected mice, five underwent awake imaging followed by 0.7% isoflurane anesthetized imaging, and three were imaged at 0.5% isoflurane or 1% isoflurane for laser power and negative BOLD controls (Fig. 2B; Supplemental Fig. S2B). Of the nine Thy1-ChR2 mice, three were imaged at 0.5% isoflurane for near-microscopic imaging (Fig. 5), four were imaged at 0.5% isoflurane for negative BOLD controls (Supplemental Fig. S2), and two were imaged at 0.7% used for the blood pressure control, a yellow light control, and a paw stimulation control.
For animals that underwent awake imaging, 3 consecutive days of acclimation to head restraint and scanner noise were first performed, followed by fMRI data collection starting on the fourth day. Past experiments in awake head-restrained mice have shown that, after a few days of such acclimation, overt stress is greatly reduced, as observed behaviorally (e.g., little or no struggling, eye secretions indicative of stress in mice, or excessive vocalizations or other signs of stress) (Boyden and Raymond 2003; Boyden et al. 2006). Each session was performed at approximately the same time of day, around 2:00 PM to minimize the influence of circadian rhythm. The headposted animal was inserted into the custom built G-10 fiberglass MRI cradle and restrained for ∼20 min to acclimate the animal to this head restrained position. The animal was given chocolate sprinkles while restrained on the MRI cradle. After 20 min of head restraint acclimation, the cradle was inserted into the 9.4-T MRI scanner. The cradle was locked to the MRI scanner stage to position the head of the animal in the center of the magnet bore. An SE-EPI, similar to the one used for fMRI data collection, was run for 40 min to acclimate the animal to scanner noise. The animal was taken out of the scanner and given two more chocolate sprinkles before finally being taken out of the body restraint tube and the headpost holder.
fMRI data collected during all three acclimation sessions were used to estimate head motion of the animal being scanned. Animal head motion was computed using the 3dvolreg program in AFNI. The 3dvolreg program assumes a rigid body transformation and minimizes a least squares difference between each source sub-brick volume and the base. We find the animal head motion reduced significantly [Supplemental Fig. S1D; F(1,48) = 261.82, P < 0.0001, main effect of acclimation in a 2-way ANOVA of acclimation × orientation of motion; F(5,48) = 59.38, P < 0.0001, main effect of orientation of motion in a 2-way ANOVA of acclimation × orientation of motion; F(5,48) = 34.122, P < 0.001, interaction between acclimation and orientation of motion] following 3 days of acclimation compared with unacclimated animal (1st day of acclimation). Changes in signal-to-noise (SNR) ratio for EPI scans were also computed for all 3 days of acclimation. Temporal SNR (tSNR) was computed from the scans without optical illumination from 100 voxels in the SI barrel field region (contralateral to the optical illumination). tSNR was determined from the mean voxel value across 200 time points, divided by temporal SD of that voxel's signal. The change in average tSNR after 3 days of acclimation compared with the first day of acclimation is plotted in Supplemental Fig. S3E. There was a trend toward an increase in tSNR following 3 days of acclimation (paired t-test, n = 5 mice; P = 0.12). During awake opto-fMRI, an RF transmit-receive coil (single copper loop, milled out from a copper plated 1/16-in epoxy material (FR-4, T-Tech, Norcross, GA), shown in Supplemental Fig. S1A) was positioned over the animal's head, surrounding the headpost. The animal was attached by its headpost to the headpost holder, which was in turn attached to the MRI cradle. The animal was positioned in a body restraint tube—a 4.5 cm diam, 2 mm wall thickness, plastic tube with an enclosed end at the tail end of the tube, padded with an absorbent material—that lightly fits around the animal's body and restricts motion during opto-fMRI experiments. The restraint tube was bolted to the MRI cradle. A 200 μm optical fiber (Ocean Optics, Dunedin, FL) attached to an adjustable optical fiber holder was also, at this time, positioned directly on the target (i.e., into the craniotomy, right above the brain), passing through the RF coil. A nose-cone for isoflurane delivery was positioned around the animals' snout. The MRI cradle was slid into the magnet bore. The cradle was locked to the MRI scanner stage to position the head of the animal in the center of the magnet bore. Animals were anesthesia induced with isoflurane right after the awake functional imaging experiment. The anesthesia level for these animals was maintained, as indicated by end-tidal isoflurane level at 0.7%.
For purely anesthetized imaging, animals were anesthetized with isoflurane (2–3% in oxygen) before positioning them on the MRI cradle. The MRI cradle was slid into the magnet bore. The cradle was locked to the MRI scanner stage to position the head of the animal in the center of the magnet bore. Breathing rate (Small Animal Monitoring 1025, SA Instruments, Stony Brook, NY) and end-tidal expired isoflurane (V9004 Capnograph Series, Surgivet, Waukesha, WI) were continuously monitored during awake and anesthetized imaging experiments. The anesthesia level was set at 0.5 or 1.0%, as indicated by end-tidal isoflurane level.
Functional images were collected at voxel resolution 150 × 150 × 500 μm (5 slices; mouse images in Figs. 1B, 2D, 3A, and 4A; Supplemental Figs. S3, A and B, S4, and S5A), 100 × 100 × 500 μm (5 slices; mouse images in Fig. 5, A and B; Supplemental Fig. S2A), or 200 × 200 × 500 μm (10 slices; mouse images in Fig. 2, A and C; Supplemental Figs. S1, B and C, and S3, B and C). For simplicity, a summary of fMRI parameters used for different experiments throughout this study is presented in Supplemental Table S1. For the high-resolution imaging sessions (that is, <150 × 150 × 500 μm in voxel size), we took advantage of the speed and strength of our gradient coils to frame the volume imaged with four saturation slices (e.g., Fig. 5A), thus avoiding wrap-around artifacts (Wang et al. 2004; Wilm et al. 2007). Functional data were acquired using a SE-EPI; 2.5 s repetition time (TR) and 25 ms echo time (TE). High-resolution T1-weighted anatomical images (78 × 78 × 500 μm) were acquired using a rapid acquisition process with relaxation enhancement (RARE) sequence in the coronal orientation, after physiological data acquisition was complete. Statistical maps of the correlation of BOLD percent signal change to the boxcar pattern of light delivery (aka boxcar correlation maps) were aligned to the high-resolution T1-weighted structural volumes using AFNI's align_epi_anat.py program.
The optical protocol used for opto-fMRI in this paper was a 10 s baseline period of darkness followed by 16 repetitions of 15 s on, 15 s off periods of 40 Hz trains of 8 ms laser pulses (the boxcar pattern shown in Fig. 1A, top right), applied to SI. Light was delivered with a 473 nm blue laser (Shanghai Dream Lasers, Shanghai, China), placed outside the magnet room and coupled to a 200 μm optical fiber (∼5 m in length), which was passed inside the magnet room through a small duct. A MATLAB program and a USB Data Acquisition Module (Cole-Parmer, Vernon Hills, IL) were used to control the laser to deliver the protocol boxcar pattern (Fig. 1A, top right).
For mapping experiments in the five awake mice (and the 0.7% isoflurane states) and for the near-microscopic opto-fMRI experiments, the laser power out the fiber was 5–10 mW (fiber tip irradiance, 150–300 mW/mm2). For control experiments, three virally injected mice, three Thy1-ChR2 mice, and two wild-type mice, we used 10–15 mW laser power; these last two wild-type mice were exposed to 25–30 mW as controls.
Three control experiments were performed on two 0.7% isoflurane-anesthetized Thy1-ChR2 mice to examine 1) blood pressure dependence on SI illumination, 2) fMRI response to yellow light, and 3) fMRI response to paw stimulation. For the first of these control mice, blood pressure, heart rate, body temperature, and respiration rate were continuously monitored (sampling rate of 5 s; Advisor, Vital Signs Monitor, Surgivet, Waukesha, WI) throughout the length of an opto-fMRI experiment, on one mouse, to observe the effect of laser stimulation on the animal's physiology. The femoral artery was cannulated with polyethylene tubing (ID 0.28 mm, OD 0.61 mm) to measure blood pressure and heart rate. The mean arterial blood pressure (mabp), mean heart rate (hr), mean temperature (t), and mean breathing rates (br) were 89.4 ± 8.7 mmHg, 577.4 ± 12.8 beats/min, 35.8 ± 1.5°C, and 98.7 ± 3.7 breaths/min, respectively, during the on periods (mean ± SE; taken across 64 on periods, in 4 scans of 16 on periods each) and 89.3 ± 8.5 mmHg, 576.5 ± 15.7 beats/min, 35.8 ± 1.0°C, and 98.9 ± 3.4 breaths/min during the off periods. On and off were not significantly different for any of these measures (paired t-test, n = 64 on periods, each followed by an off period; Pmabp > 0.90; Phr > 0.90; Pt > 0.55; Pbr > 0.85). For the second mouse, the effect of paw stimulation was measured (15 s, 40 Hz, 3 ms pulse duration, 3 mA-amplitude pulse current delivery) on the BOLD signal, finding the percent change in SI to be, at its temporal peak, a 2.2% increase (at 7.5 s after stimulus onset), with a shape that matches the hemodynamic response function (HRF) with P < 0.0005 (least squares linear regression analysis, comparing the canonical HRF with each of the 4 scans within the session; R2 = 0.9627, F = 32.93). In this same mouse, we measured the effect of yellow light (593 nm, 5 mW) on the fMRI signal and found no induced BOLD signal changes (positive or negative).
Control experiments investigating putative negative signals were conducted using an MRI phantom of 3% agarose in saline containing a LEGO brick (Fig. 2E).
Opto-fMRI of mouse: data analysis
The analysis pipeline is outlined in the bottom right of Fig. 1A. Boxcar correlation maps for EPI functional data were generated using AFNI (Cox 1996; Nelson et al. 2006) (National Institutes of Health, http://afni.nimh.nih.gov/afni, Bethesda, MD) and MATLAB (The Mathworks, Natick, MA). EPI functional data were motion corrected in all three spatial dimensions, but no spatial smoothing or undistortion was performed (to preserve full resolution). The images were slice-time corrected and detrended, as is standard for EPI scanning (Lindquist 2008; Smith et al. 1999). Percent signal change was computed by subtracting, from each voxel's BOLD time series, the temporal average of all off periods for all scans in a session and dividing by the average of the signal across all off periods. We performed four scans per session, and percent signal change was averaged across scans in a session unless otherwise indicated.
To determine which voxels had significant increases or decreases in BOLD signal, a voxelwise time series correlation of the percent signal change was performed, as determined above, with the protocol boxcar (Fig. 1A, top right), delayed by one repetition time (TR; 2.5 s) to compensate for the BOLD response time lag (Cox 1996; Nelson et al. 2006); we label this analysis “boxcar correlation” in the text. Activation in a region was deemed significant if a cluster (here used in the fMRI sense of the word) of at least six contiguous voxels had correlation P values at an uncorrected threshold of P < 0.005. This cluster size of 6 and this uncorrected P value threshold of 0.005 were objectively chosen via Monte Carlo simulations performed in AFNI (family-wise error correction using AlphaSim program), to result in a multiple comparisons corrected type I error of 5%, appropriate for the performance of statistics for individual voxels taken throughout the entire imaging volume (Forman et al. 1995; Xiong et al. 1995). The parameters used for the four Monte Carlo simulations (1 for each scanning resolution used) were as follows: number of voxels, 100 × 100 × number of slices; voxel size (functional resolution); size of smoothing filter (same as voxel size); number of Monte Carlo iterations, 10,000.
Region of interest identification and cross-registration
To identify regions showing positive BOLD responses, boxcar correlation maps (that is, the correlation of BOLD percent signal change to the boxcar model, as described above) were aligned to the high-resolution T1-weighted structural volumes using AFNI's align_epi_anat.py program. High-resolution T1-weighted structural images were cross-registered between mice using AFNI's 3dAllineate program; this alignment was used to bring the boxcar correlation maps into a standard coordinate space.
For each imaging session, the coordinate location of the voxel with peak boxcar correlation, in each contiguous set of significantly activated voxels, was identified. Using k-means clustering we clustered (note: the term cluster is now being used in the k-means sense, not in the fMRI preprocessing sense used in the previous section of the methods) the coordinate locations of the peak-correlation voxels thus identified. We adapted an algorithm (Duda et al. 2000) to perform the k-means analysis in an unsupervised fashion, following the steps listed below, starting with C = 1:
performing k-means clustering with C clusters;
computing for each resultant cluster the within-cluster sum-of-squares position variance – that is, the sum, over all peak-correlation voxels in a cluster, of the squared distance |(coordinate of the peak-correlation voxel) – (mean coordinate of all peak-correlation voxels within the cluster)2|;
comparing the C sum-of-squares variances thus computed, using a t-test, to those obtained when k-means was run with C + 1 clusters.
We iterated this process, increasing C by 1 each time until the P value of this t-test was no longer significant (P < 0.05); then, the number of clusters used for k-means was set to C.
Key regions of interest (ROIs) were named overlaying T1 anatomical images onto corresponding atlas plates from the Paxinos and Franklin (2001) atlas, by matching their anteroposterior coordinates and scaling atlas plates along the x- and y-axes until the borders of the T1 anatomical images and the atlas plate corresponded.
As a test of robustness, the k-means clustering analysis was repeated on correlation maps obtained for two different values of boxcar correlation P value threshold (5 times higher and 5 times lower than the 1 objectively chosen via Monte Carlo simulation to result in a multiple comparisons corrected type I error of 5%; Supplemental Fig. S5).
ROI time series extraction and analysis
Extracted time series from each key ROI were signal averaged across the contiguous significantly activated voxels associated with each key ROI and averaged across pulse trains within a session (consisting of 4 scans, of 16 pulse trains each) to generate an average time series of percent change in the BOLD signal for each key ROI (Fig. 4B). The peak percent change in the BOLD signal, across time, was extracted from these thus-averaged traces (Fig. 4C).
We assessed the extent to which BOLD signal time series could be fit with a human-derived canonical HRF. We produced a goodness-of-fit (R2) measure and probability that the fit is statistically significant, by applying a least squares linear regression for each key ROI, in each of the awake and anesthetized (0.7% isoflurane level) states, each time taking into account data from each of the five animals to compare against the canonical HRF. In a separate analysis, we sought to test whether the goodness-of-fits themselves differed between the awake and anesthetized states or between regions. For this analysis, the goodness-of-fit (R2) for the default value of delay of response (6 s, relative to onset) was first computed on each animal's individual time series (averaged across scans within a session) and then through Fisher's z-transform computed the z-score for the goodness-of-fit to the HRF. The resultant set of z-scores was subjected to a two-way ANOVA with factors of brain state and ROI. In a second analysis, we varied the delay of response parameter from 1 to 12 s. The delay of response for which the R2 value was maximum was computed for each animal independently, and the resultant set of delay of responses that resulted in peak R2 values was subjected to a two-way ANOVA with factors of brain state and ROI.
Normalized cross-correlations of BOLD percent signal change time series were computed between pairs of key ROIs (Fig. 4D) (Haralick and Shapiro 1992). We computed, for each key ROI, the average of all the normalized cross correlations between that ROI and all the other key ROIs; this average was termed the “total connectivity” of that region.
BOLD signal responses as a function of cortical depth were computed by drawing a line perpendicular to the cortical surface and choosing sets of three voxels (the 1 closest to the line and the 2 horizontally flanking ones) for each cortical depth. Peak percent BOLD signal change was averaged across each of these sets of three voxels and across animals (n = 3) for each depth below the surface of the cortex (Fig. 5C).
We used t-test and ANOVA (StatView, SAS Institute, Cary, NC) throughout, e.g., to analyze the dependence of BOLD signal change or significantly activated brain volume on anesthesia and brain region (Fig. 4C, etc.).
Local field potential and multiunit activity recordings
Two Thy1-ChR2 mice were anesthetized with isoflurane (0.5%) and headposted for electrophysiology, followed by the opening of a small craniotomy over the barrel cortex. A laminar silicon linear electrode array (Neuronexus Technologies) was connected to a Cheetah32 data acquisition system (Neuralynx). The array consisted of 16 contacts along the probe shank, each with a diameter of 15 μm and spaced 100 μm apart. The probe shank was inserted such that ≥12 contacts covered the cortical laminae from immediately below the surface to a depth of ∼1,200 μm continuously (depicted in Fig. 5D). Recordings were performed with minimal filtering (band-pass 0.1–9,000 Hz) and separated into local field potential (LFP) and multiunit activity (MUA) signals in software using MATLAB. Specifically, MUA was determined by filtering the original signal between 600 and 6,000 Hz and thresholding at 4 SD above background noise; LFP was determined by filtering the original signal between 7 and 200 Hz. MUA and LFP plots in Fig. 5, E and F, were averaged across two animals, for 10 trials of 1 s of light stimulation (40 Hz, 8 ms pulses) each time, and plotted as a function of cortical depth. LFP power was normalized to the 1 s baseline before light onset.
After opto-fMRI experiments were concluded, we verified the extent of ChR2 expression. Animals were transcardially perfused with 100 mM PBS followed by 4% formaldehyde in PBS. Brains were postfixed overnight at 4°C and cryoprotected in 30% sucrose for 24–48 h before sectioning. Free-floating sections (50 μm) were cut using a cryostat (Leica 3050S, Leica Microsystems, Wetzlar, Germany), mounted on glass slides with Vectashield mounting medium with DAPI (Vector Laboratories, Burlingame, CA), and coverslipped. Spread and labeling efficiency were estimated by examination of 50 μm coronal sections near the site of injection for the presence of GFP. For the five animals that underwent the core awake opto-fMRI experiments, four were histologically examined; the fifth mouse (CaMKII-ChR2-D) was not examined because of potential health issues unrelated to the opto-fMRI experiment. Examination of the entire brain between V1 and frontal cortex showed no GFP beyond the site of injection. To estimate the volume of the cortex containing ChR2-GFP expressing cells, we measured the area in each slice by measuring image intensity of GFP (field of view 1.9 × 1.9 mm, 4× objective, BX51, Olympus, Tokyo, Japan). Each injection of ∼1 μl labeled an ellipsoid volume with a diameter of 700–1,000 μm along the anterior-posterior aspect, 700–1,000 along the medial-lateral aspect, and 800–1,100 mm along the dorsal-ventral aspect. Next, we collected confocal image stacks (21 × 21 × 45 μm; FV10i, Olympus) sampling two regions along the dorsal-ventral axis of each GFP-positive area, on each of three equally spaced slices sampled from each animal. For each of these regions, we counted the number of GFP-positive and DAPI-positive cell bodies. By this measure, the local density of cells expressing ChR2 would be estimated as 314,938 ± 110,485/mm3 (SE) over the regions identified (n = 4 mice).
Optimization of opto-fMRI methods for awake mice
A diagram of the opto-fMRI hardware for awake mouse imaging is shown in Fig. 1A. This equipment consists of a blue laser coupled to a long optical fiber, a custom-engineered MRI cradle for holding an awake mouse (CAD design available on request), a body restraint tube within the cradle for minimizing body motion, an RF coil (Supplemental Fig. S1A), and an optical fiber and holder (which passes through the open loop of the RF coil) over the desired neural target. To facilitate placement of the animal in the MRI cradle, the animal is surgically prepared in advance with a plastic head post that is bolted into the cradle at the beginning of an experiment. The animal is also prepared with a craniotomy or thinned skull over the area of interest (in this case, the SI barrel cortex, which has been transgenically or virally labeled to express ChR2 in pyramidal neurons). The experimental opto-fMRI protocol (Fig. 1A, right) used EPI in conjunction with a boxcar protocol of light delivery (15 s periods of 40 Hz delivery of 8 ms laser pulses, separated by 15 s periods of no illumination). Significantly activated voxels in images were identified by correlating the percent signal change of the BOLD response with a boxcar function that corresponded to the periods of light delivery and selecting the voxels whose correlation coefficients were significant at a familywise multiple comparisons-corrected P value level of <0.05 (equivalent to an uncorrected P value of 0.005) and had five or more significant neighbors (i.e., a cluster threshold of 6). The uncorrected P value of 0.005 and the cluster threshold of 6 were objectively determined using the Monte Carlo simulation program of AFNI (Forman et al. 1995; Xiong et al. 1995).
Our method was optimized along a number of experimental axes. First, to enable fMRI in awake mice, it was critical to acclimate mice to head restraint and scanner noise over a multi-day period under conditions similar to actual opto-fMRI experimentation. In addition, the body restraint tube was important for reducing body motion. Although the body was not in the fMRI field of view, motion can cause susceptibility artifacts in brain images by changing the magnetic environment. Several other methodological choices were important for preventing magnetic susceptibility artifacts. First, although GRE-EPI is commonly used because of its excellent signal-to-noise properties, we found that SE-EPI minimized distortion caused by magnetic susceptibility mismatch between tissue and air (Bandettini et al. 1994). This sequence approach (SE-EPI) also provides for higher-resolution imaging (Kim et al. 2000; Lee et al. 1999; Moon et al. 2007) (Supplemental Fig. S1C). Second, adding a thin layer of dental cement homogeneously across the surface of the skull reduced magnetic susceptibility artifacts caused by the mismatch between intact tissue and air (Supplemental Fig. S1B). Third, we adapted a strategy from diffusion tensor imaging (Wang et al. 2004; Wilm et al. 2007), namely, the use of saturation slices, to eliminate signals from nearby tissue outside the field of view.
BOLD signal changes downstream of SI pyramidal cell activation
Opto-fMRI was performed on awake mice that had undergone viral infusion of the excitatory neuron-targeting lentivirus FCK-ChR2-GFP (Han et al. 2009) into the SI barrel cortex. These mice expressed ChR2-GFP in pyramidal neurons within a spherical volume between 700 and 1,000 μm in diameter, containing neurons distributed in layers 2/3 and layer 5 (Fig. 1Bi). On illumination of SI pyramidal cells, we observed robust, distributed BOLD responses in the brain ipsilateral and contralateral to the fiber illumination site, time locked to the illumination period. Figure 1Bii shows the results from a representative awake mouse. The time series of evoked BOLD signals in SI is plotted for this mouse in Fig. 1Biii (red trace). When the laser was not on, there was no visible change in the SI BOLD signal (Fig. 1Biii, black trace). We did not detect positive BOLD responses in control mice not expressing ChR2 undergoing opto-fMRI (Fig. 2A; n = 2 mice; representative mouse shown), even with laser powers increased to levels threefold to sixfold higher than used in the experiments here described (Fig. 2B). Few voxels were significantly activated in control even at an uncorrected P value of 0.15 (Fig. 2B). This control experiment mitigates the possibility that positive BOLD responses were produced by non-ChR2-related effects of cortical laser illumination, such as any potential temperature changes induced by illumination.
A putative negative BOLD response was observed in the brain during opto-fMRI in the SI barrel cortex immediately under the optical fiber tip (Fig. 2D). This negative BOLD response was similar in ChR2-negative mice (Fig. 2C). To probe the mechanism of negative signal, we carried out a specific set of experiments to compare the negative signal under the fiber tip in three ChR2-virally transduced mice, two non-ChR2-expressing mice, and four Thy1-ChR2 mice and found negative responses in all three sets of mice of comparable magnitude [Fig. 2F; P > 0.80, F(2,6) = 0.225 for peak reduction in signal, factor of mouse type, 1-way ANOVA]. The fact that lack of ChR2 did not alter the signal reduction, despite its elimination of the signal increases observed (e.g., Fig. 2, A and B) led us to hypothesize that the signal reduction observed under illumination was not related to ChR2 but instead reflected a direct and local effect of light on the sample that was observable using MRI. In support of this hypothesis, we observed small but significant signal reductions even in plain agarose (Fig. 2E), on illumination. This local effect was at least in part unrelated to optical activation of ChR2-expressing neurons and may therefore reflect temperature-related modulation of the fMRI signal (Yablonskiy et al. 2000) or other nonspecific local effects of illumination (see Supplemental Results and Supplemental Fig. S2 for more analysis of the negative signal, i.e., analysis for Thy1-ChR2 mice). We confine the rest of our analysis to observation of positive BOLD signals, because negative signal was not consistently observed except for this nonspecific local response.
To characterize the consistency of the opto-fMRI map across different mice, we performed opto-fMRI on five awake mice, each injected with 1 μl FCK-ChR2-GFP lentivirus in SI of the left hemisphere. To facilitate this comparison, the resultant functional scans were aligned for each mouse to its respective T1 anatomical scan, and these datasets were cross-registered across all mice. All mice exhibited BOLD activations in the illuminated SI and also displayed activations in nearby cortical regions, subcortical regions, and the contralateral cortex. Voxels that were significant in any of the five mice are color-coded according to the median (taken across mice) boxcar correlation value in Fig. 3A (raw datasets for all 5 mice are shown in Supplemental Fig. S4).
Neural targets recruited by SI pyramidal cell activation can be automatically identified
To identify, in an unbiased manner, the brain regions downstream of SI activation (Figs. 1Bii and 3A), we developed an unsupervised algorithm that first clusters sets of contiguous significantly activated voxels into ROIs and then localizes the centroid of each ROI to a reference mouse atlas. The clustering algorithm is based on an automated variant of the k-means clustering algorithm (see methods for details), which adjusts the number of clusters according to objective statistical criteria, and showed 11 ROIs for the five opto-fMRI experiments conducted (Fig. 3Bi). Of these 11 ROIs, we designated 5 as key ROIs, which were robustly activated in all five mice (Fig. 3Bii). The centroid-atlas alignment showed these key ROIs as ipsilateral SI, contralateral SI (c-SI), ipsilateral primary motor cortex (MI), ipsilateral secondary somatosensory cortex (SII), and the caudoputamen (CP), which are regions known to be targets of SI pyramidal neurons (Aronoff et al. 2010; Carvell and Simons 1986, 1987; Chakrabarti and Alloway 2006; Diamond et al. 2008; Ferezou et al. 2007; Megevand et al. 2008; White and DeAmicis 1977). The algorithm also reported ROIs that were less consistent across individual mice (Fig. 3Bii). These ROIs might represent additional activated circuits, for which the reliability of observation was not as high as the selected five key ROIs. For example, in a subset of the five opto-fMRI experiments, clusters emerged in regions that were atlas-aligned to thalamus and contralateral SII (Fig. 3Bi; these activated regions were also visible in the raw datasets shown in 11 ROIs; Supplemental Fig. S4). This variability could result from natural variations in the strength of connectivity of circuits from one brain to another, from subtle variations in the location of virally labeled neurons or of the placement of the optical fiber across mice, or from the fact that, in our algorithm, the choice of a boxcar correlation statistical threshold could admit significant voxels for some mice but reject others, leading to a perceived heterogeneity across mice.
To assess statistical robustness, we repeated the opto-fMRI algorithm on boxcar correlation maps obtained by varying the uncorrected boxcar correlation P value to values 5 times higher and 5 times lower than that used in Fig. 3. We found that varying the uncorrected boxcar correlation P value to 0.001 or 0.025 resulted in few qualitative changes in the appearance of the correlation maps (representative mouse shown in Supplemental Fig. S5Ai; median boxcar correlation map shown in Supplemental Fig. S5Aii). When we ran the unsupervised neural target identification algorithm on the opto-fMRI data, few differences emerged in the cluster map, with key ROIs preserved (Supplemental Figs. S5, B and C) and with most clusters remaining centered at the same location but expanding or contracting as the P value was relaxed or made stricter. Our algorithm thus exhibits robustness to variations in the specific P value chosen for the initial boxcar correlation threshold, indicating its utility in comparing neural maps across animals and across experiments.
Opto-fMRI evokes a greater BOLD response in awake than anesthetized animals
We anticipate a major use of opto-fMRI will be in awake animals, given the potential for comparing animal opto-fMRI data to human fMRI data (which is predominantly collected in the awake state) and for understanding behavior in head-posted rodents, an approach that has recently been optimized for a variety of imaging modalities (Andermann et al. 2010). To show the use of opto-fMRI in analyzing how light-activated neural networks can be modulated, we compared, for the five mice analyzed in Fig. 3, how administration of 0.7% isoflurane varied the network impact of SI illumination. Isoflurane was chosen because it is commonly used in small animal fMRI to prevent motion artifacts and immobilization stress (King et al. 2005; Lahti et al. 1998; Sicard et al. 2003), and it was also the anesthetic used in an earlier opto-fMRI study (Lee et al. 2010). The role that isoflurane anesthesia (or any anesthesia) has on local versus global network dynamics, and on the coupling of these neural signals to blood flow, is controversial (Corfield et al. 2001; Matta et al. 1999; Tsurugizawa et al. 2010).
To visualize the effects of anesthesia on the set of regions causally recruited by SI, we plot in Fig. 4A the clusters or ROIs identified by the algorithm described above, for the awake (top) and anesthetized (bottom) states. Raw datasets leading up to these ROI plots, akin to those shown for the awake state in Fig. 3, are provided in Supplemental Fig. S3. In the anesthetized state, only 6 ROIs were activated compared with the 11 ROIs in the awake state. Of the five key ROIs identified in the awake state (SI, c-SI, SII, MI, and CP), four were represented in the anesthetized state, but the striatal response was less prominent. This finding is consistent with previous studies using electrophysiology (Detsch et al. 2002; West 1998) reporting that anesthesia suppresses somatosensory-evoked information transfer to subcortical regions. In the isoflurane anesthetized condition, there was also a decrease in the BOLD response for the voxels in all five key ROIs identified in the awake state [Fig. 4, B and C; F(1,40) = 172.63, P < 0.0001, main effect of anesthesia level in a 2-way ANOVA of anesthesia level × region; F(4,40) = 7.830, P < 0.0001, main effect of region in a 2-way ANOVA of anesthesia level × region; no interaction, F(4,40) = 0.558, P = 0.6943]. These contrasts show the potential utility of opto-fMRI for probing dynamics in network activation.
To address the question of how the time course of the BOLD response is altered by anesthesia administration, we examined the goodness-of-fit of BOLD responses in different regions under different levels of anesthesia to the canonical HRF. The BOLD response in SI was well fit with the canonical HRF in the awake state [Fig. 4B, dotted line superimposed on red trace in left plot; R2 = 0.94, F(1,4) = 21.37, P = 0.0002, n = 5 mice; least square linear regression, allowing only the amplitude scaling of the HRF to vary, as well as the amplitude offset, but not allowing for scaling or shifting of the HRF along the time axis]. Indeed, all of the awake responses were well fit with the canonical HRF [R2 ranging from 0.78 to 0.94, F(1,4) ranging from 7.28 to 21.37, P values ranging from 0.0002 to 0.009; dotted lines superimposed on red trace in each plot of Fig. 4B]. The anesthetized responses were less well fit by the HRF [R2 ranging from 0.15 to 0.90, F(1,4) ranging from 1.57 to 17.57, P value ranging from 0.0008 to 0.27].
To apply a direct statistical comparison of the goodness-of-fit to the HRF for the awake and anesthetized BOLD signal changes, we computed an R2 for every opto-fMRI session BOLD signal versus the HRF. In our first analysis, we computed this R2 using the default value of delay in HRF (6 s relative to onset) and then through Fisher's z-transform computed the z-score for the goodness-of-fit to the HRF. We ran an ANOVA to test how these z-scores varied across anesthesia level and region. We found that anesthesia significantly lowered the z-score for the goodness-of-fit to the HRF of the BOLD signal relative to the awake state [F(1,40) = 167.62, P < 0.0001, main effect of anesthesia level in a 2-way ANOVA of anesthesia level × region], and regions also differed in their match to the HRF [F(4,40) = 34.42, P < 0.0001, main effect of region in a 2-way ANOVA of anesthesia level × region]. There was a significant interaction between anesthesia level and region [F(4,40) = 7.28, P = 0.0002]. Specifically, local SI activity was least affected in its BOLD signal response shape by anesthesia (z-score for the goodness-of-fit to the HRF decreased from 2.07 to 1.59), whereas the striatum was most affected (z-score for the goodness-of-fit to the HRF decreased from 1.48 to 0.35). In a second analysis, we varied the delay parameter for the HRF from 1 to 12 s and recomputed the R2 for each fit to determine the delay that would result in the highest R2 value. We found that anesthesia significantly increased the delay in HRF that resulted in the peak R2, relative to the best-fit delay seen for the awake state [F(1,40) = 98.33, P < 0.0001, main effect of anesthesia level in a 2-way ANOVA of anesthesia level × region], and regions also differed in their delay in HRF that would result in peak R2 [F(4,40) = 11.04, P < 0.0001, main effect of region in a 2-way ANOVA of anesthesia level × region]. There was a significant interaction between anesthesia level and region [F(4,40) = 18.96, P < 0.0001] in this delay computation.
Anesthesia affects BOLD functional connectivity in the network recruited by pyramidal cell activation
We applied opto-fMRI to probe the functional connectivity of the brain network downstream of SI as a function of brain state, again using isoflurane anesthesia as a model. We performed normalized cross-correlation analyses of the BOLD signal changes in different ROIs during opto-fMRI. This technique provides a measure of temporal correlation across regions while nullifying the effect of absolute changes in signal amplitude (as reported above), allowing focus on how anesthesia level modulates the coordination between regions. We found that the anesthetized state had generally lower normalized cross-correlations between pairs of regions than did pairs of regions in the awake state [Fig. 4D; F(1,80) = 71.67, P < 0.0001, main effect of anesthesia level in a 2-way ANOVA of anesthesia level × region pair]. This measure of correlation also varied between different pairs of regions [F(9,80) = 4.76, P < 0.0001, main effect of region in a 2-way ANOVA of anesthesia level × region pair, and no interaction between the two factors, F(9,80) = 1.22, P = 0.29].
The striatum showed lower correlations with other regions than did other pairings. The 17 post hoc tests of differences between the awake and anesthetized state that showed P values using Fisher's protected least significant difference (PLSD) test significant with respect to a 0.05 significance level, all involved the striatum. When we computed, for each region, the average normalized cross-correlation between the striatum and all the other regions, we found that this “total connectivity” was lowered by anesthesia [F(1,40) = 86.50, P < 0.0001, main effect of anesthesia level in a 2-way ANOVA of anesthesia level × region pair]. Different regions also exhibited different total levels of connectivity with the rest of the network [F(4,40) = 3.39, P < 0.02, main effect of region in 2-way ANOVA of anesthesia level × region pair; no interaction, F(4,40) = 0.76, P > 0.55], with the striatum exhibiting less total connectivity than the others (P < 0.02, Fisher's PLSD, comparison with all other regions).
Opto-fMRI applied to transgenic mice expressing ChR2 in defined cortical layers
Opto-fMRI can be applied to other models of potentially significant interest, and we explored these signals in transgenic mice. These preparations provide disease models and have been important in a variety of basic science studies. We applied opto-fMRI to Thy1-ChR2 transgenic mice that express ChR2-YFP in the neocortex in layer 5 pyramidal neurons (Arenkiel et al. 2007) (Fig. 5Di). In this mouse strain, multiunit neural activity appeared predominantly in layer 5 (Fig. 5, Dii and E); LFP power was broadly distributed in these mice, perhaps because of the extensive ChR2-positive dendrites (Fig. 5F). When Thy1-ChR2 mice received light isoflurane (0.5%) and underwent opto-fMRI, positive responses could be observed in SI sensory cortex (Fig. 5B), including prominent activity in layer 5. The ChR2-independent negative signal in Thy1-ChR2 mice (Supplemental Fig. S2) also was observed in the most superficial layers of the cortex (Fig. 5C). These findings extend the application of opto-fMRI to transgenic mice, suggesting that application to other lines will be successful.
Optogenetic strategies using ChR2 expression have previously been shown to evoke a BOLD response in the anesthetized rat (Lee et al. 2010). Here we extend this work is three ways. First, robust positive BOLD responses were shown in the awake mouse. The ability to perform opto-fMRI in the awake mouse opens up translational links between causal circuit studies in animals and behaviorally and clinically important human studies. Second, using photoactivation of cortical pyramidal cells of the mouse SI barrel cortex, which has well-characterized targets throughout the brain (Aronoff et al. 2010; Carvell and Simons 1986, 1987; Chakrabarti and Alloway 2006; Diamond et al. 2008; Ferezou et al. 2007; Megevand et al. 2008; White and DeAmicis 1977), we showed that distributed regions of the network could be evoked by light stimulation and measured using an unbiased, automated analysis. The network response included regions in the striatum: the MI, c-SI, and SII. Third, functional connectivity measures showed the BOLD response was correlated between regions within the network consistent with observations from evoked and intrinsic BOLD fluctuations in the human (Fox and Raichle 2007; Friston 1994; Hampson et al. 2006; Pezawas et al. 2005; Ramnani et al. 2004; Stephan et al. 2008; Van Dijk et al. 2010). The absolute magnitude of BOLD response in each region of the network and the functional correlation strengths between regions were attenuated by anesthesia. These collective results suggest that opto-fMRI provides an approach to map circuits in the awake mouse and thus complements fMRI approaches used in vivo in humans.
This study also provides a detailed account, intended to facilitate the implementation of this method by future researchers. A number of technical challenges were overcome to achieve the combination of optical neural activation and fMRI in the awake mouse. Specialized hardware for awake mouse imaging was implemented (Fig. 1) that enables positioning of optical fibers and immobilization without incurring significant motion artifacts, using a combination of head posts and restraint tubes. Surgical and pulse sequence strategies were devised for reducing magnetic susceptibility artifacts while preserving high resolution. For the majority of the experiments in this study, we chose to use 150 mW/mm2 irradiance, well within the range used for repeated stimulation (Ayling et al. 2009; Cardin et al. 2009, 2010; Kravitz et al. 2010; Zhang et al. 2010) and at the lower limit of what is needed to evoke EMG responses when targeting MI (Ayling et al. 2009); however, it will generally be useful to select the least dosage of light necessary for achieving a given scientific goal. We also performed a number of critical control experiments, for example, showing that negative signal changes under the optical fiber are nonspecific effects of opto-fMRI, unrelated to the activity of ChR2-expressing neurons. Future studies will be required to differentiate these nonspecific signals from negative BOLD responses.
Opto-fMRI provides a strategy for characterizing the network downstream of a target in a fashion that is amenable to repeated assessment, potentially over extended time periods suitable to longitudinal study of learning and plasticity, and disease progression. Thus this work complements recent studies by Lee et al. (2010) in the rat discussed above and Logothetis et al. (2010) in the monkey that used electrical stimulation to measure the BOLD response in distributed regions across the visual system. In the future, the addition of new technologies such as optical neural silencing reagents (Chow et al. 2010; Han and Boyden 2007; Zhang et al. 2007; Zhao et al. 2008) or multisite illuminators, the targeting of different neuron types, and further improvements of fMRI, will augment the performance of opto-fMRI. Awake mouse opto-fMRI also opens up the possibility of implementing a psychophysical paradigm and assessing brain-wide effects of changing the psychometric threshold curves through optogenetically manipulating a specific subcircuit (Cardin et al. 2009).
There are also clear limitations to opto-fMRI. Weak, modulatory functional connections or transient connections may not be detectable because of the slow temporal response of the BOLD signal and the nature of neural-to-hemodynamic coupling. This limitation of poor temporal resolution is in part offset by the utility of this approach for exploring widespread distributed regions and the specificity of the target cell types and location that can be achieved. Another limitation of opto-fMRI method is the need for a head-fixed animal. Head-fixed animals have been shown to have sparse firing rate in the sensory cortex compared with freely moving animals (de Kock and Sakmann 2009; Sakata and Harris 2009).The higher firing rate is thought be important for information processing in various ways, in a behaviorally relevant state (Vijayan et al. 2010).
Other approaches for mapping neural circuits in a causal fashion using fMRI exist. Electrical microstimulation has proven to be an important tool in combination with fMRI (Canals et al. 2008; Ekstrom et al. 2009; Logothetis 2003; Logothetis et al. 2010; Sultan et al. 2007; Tolias et al. 2005) for driving local circuits in the context of whole-brain mapping, but electrical activation can recruit a diversity of cell types within a microcircuit, as well as fibers of passage, smooth muscle, and other brain circuit elements (Ranck 1975). However, the passive spread of current can confound the analysis of functional connections (local and distal), as determined by fMRI, requiring indirect methods to characterize the current spread such as behavioral controls (Tehovnik et al. 2006). Future work that directly compares stimulation methods will be valuable for understanding tradeoffs between methods and resolve the uses of these tools in different spheres of neuroscience.
E. S. Boyden acknowledges funding by the National Institutes of Health Director's New Innovator Award (DP2 OD-002002-01), National Institutes of Health Challenge Grant 1RC1 MH-088182-01, National Institutes of Health Grand Opportunities Grant 1RC2 DE-020919-01, National Institutes of Health Grant 1R01 NS-067199-01, National Science Foundation (NSF) Grants 0835878 and 0848804, the McGovern Institute Neurotechnology Award Program, the Department of Defense Congressionally Directed Medical Research Program Post-Traumatic Stress Disorder Program, National Alliance for Research on Schizophrenia and Depression, the Alfred P. Sloan Foundation, Jerry and Marge Burnett, the Society for Neuroscience Research Award for Innovation in Neuroscience, the Massachusetts Institute of Technology Media Laboratory, the Benesse Foundation, and the Wallace H. Coulter Foundation. A. M. Graybiel acknowledges funding by National Institutes of Health Grants MH-060379 and NS-025529 and European Union Grant 201716. C. I. Moore acknowledges funding by the NSF Small Grant for Exploratory Research (SGER) program and McGovern Institute for Brain Research. N. Kopell acknowledges the NSF SGER program. I. Kahn and R. L. Buckner acknowledge the Howard Hughes Medical Institute.
No conflicts of interest, financial or otherwise, are declared by the authors.
We thank J. Siegle for designing the mouse head post, J. Cardin for technical assistance with histology, and the Harvard Center for Brain Science Optical Imaging Core and S. Turney for assistance with fluorescence imaging.
Present address of I. Kahn: Department of Physiology and Biophysics, Faculty of Medicine, Technion, Haifa, Israel.
↵1 The online version of this article contains supplemental data.
- Copyright © 2011 the American Physiological Society