One of the basic properties of sensory cortices is their topographical organization. Most imaging studies explored this organization using the positive blood oxygenation level-dependent (BOLD) signal. Here, we studied the topographical organization of both positive and negative BOLD in contralateral and ipsilateral primary somatosensory cortex (S1). Using phase-locking mapping methods, we verified the topographical organization of contralateral S1, and further showed that different body segments elicit pronounced negative BOLD responses in both hemispheres. In the contralateral hemisphere, we found a sharpening mechanism in which stimulation of a given body segment triggered a gradient of activation with a significant deactivation in more remote areas. In the ipsilateral cortex, deactivation was not only located in the homolog area of the stimulated parts but rather was widespread across many parts of S1. Additionally, analysis of resting-state functional magnetic resonance imaging signal showed a gradient of connectivity to the neighboring contralateral body parts as well as to the ipsilateral homologous area for each body part. Taken together, our results indicate a complex pattern of baseline and activity-dependent responses in the contralateral and ipsilateral sides. Both primary sensory areas were characterized by unique negative BOLD responses, suggesting that they are an important component in topographic organization of sensory cortices.
One of the basic properties of sensory cortices is their topographical organization. Recent studies have found that topographic mapping might also be a crucial component of brain organization for higher processing functions such as memory, spatial attention and numerosity (Silver et al. 2005; Hagler and Sereno 2006; Kastner et al. 2007; Harvey et al. 2013). One of the most famous manifestations of topographic organization is the somatotopic organization of the somatosensory-motor system described by Penfield and Boldery (1937) in the primary somatosensory (S1) and primary motor (M1) cortices of humans. This organization has been confirmed using noninvasive imaging techniques (Fox et al. 1987; Kurth et al. 1998; McGlone et al. 2002; Blankenburg et al. 2003; Kell et al. 2005; Chen et al. 2006; Miyamoto et al. 2006; Huang and Sereno 2007; Nelson and Chen 2008; Schweizer et al. 2008; Overduin and Servos 2008a) by showing that tactile stimulation of neighboring parts of the body evokes an increase in the blood oxygenation level-dependent (BOLD) signal in adjacent areas in the contralateral S1. However, the vast majority of these studies have only used up to 3 distinct body parts or focused on distinct areas of the body surface (Kurth et al. 1998; McGlone et al. 2002; Blankenburg et al. 2003; Chen et al. 2006; Huang and Sereno 2007; Nelson and Chen 2008; Schweizer et al. 2008; Overduin and Servos 2008a). Furthermore, phase-locked approaches, which have been proven to be extremely useful for mapping full topographic gradients, are only of limited use in the somatosensory system. These approaches have been used to map several sensory modalities such as the visual system and more recently the auditory and motor system (Engel et al. 1994, 1997; Hertz and Amedi 2010; Striem-Amit et al. 2011; Engel 2012; Besle et al. 2013; Zeharia et al. 2012, 2015). They were also used to study integration of visual and tactile responses in the human homolog of the face area ventral intraparietal (VIP) (Sereno and Huang 2006) using an automated air-puff system. Surprisingly, the whole-body sensory representation has never been tested with the high-resolution neuroimaging techniques available today. Thus, the first goal of the current study was to map full-body gradients in S1 using continuous natural tactile stimulation combined with phase-locking techniques.
The second recently observed and still controversial phenomenon that might characterize topographical mapping is the drop in the signal below baseline observed in the nonstimulated parts of the topographic gradient. This was reported in a study using simultaneous electrophysiological recording and functional magnetic resonance imaging (fMRI) in the visual cortex of macaque monkeys on visual eccentricity maps (Shmuel et al. 2006). The findings indicated that both the multiunit electrode activity and the BOLD signals exhibited a peripheral drop in the signal below the baseline condition when the fovea was stimulated and vice versa. Such negative BOLD responses, or deactivations, have been found in the visual cortex of humans as well (Shmuel et al. 2002; Bressler et al. 2007; Pasley et al. 2007) and were suggested to reflect functional decreased neural activity and/or increased inhibition (Stefanovic et al. 2004; Shmuel et al. 2002, 2006; Bressler et al. 2007; Devor et al. 2007; Kastrup et al. 2008; Boorman et al. 2010; Klingner et al. 2010; Schäfer et al. 2012). Recently, our group extended these findings to the motor domain by showing the same phenomenon in the primary motor area (M1) for synchronous bilateral body movements but not in the supplementary motor area (SMA, Zeharia et al. 2012) or anywhere else in the other 8 homunculi outside M1 (Zeharia et al. 2015). However, it is not clear whether this mechanism is also present in the somatosensory system. Therefore, our second goal was to explore negative BOLD responses in the contralateral hemisphere, and specifically to test whether such negative responses would show a topographically biased organization similar to the results reported for the visual and motor cortices. Since such topographic negative BOLD responses mostly seem to characterize primary sensory and motor areas, we focused here on the primary somatosensory area. Our prediction was that stimulation of the feet would result in negative BOLD responses in the lip area and vice versa and that stimulation of the hands would result in deactivation in both the upper and lower parts of the body gradient.
Interestingly, negative BOLD responses in the somatosensory system have primarily been studied in the context of cross-hemispheric interactions. Specifically, negative BOLD was found in the ipsilateral primary somatosensory and motor cortices after unilateral hand or finger stimulation (Hlushchuk and Hari 2006; Kastrup et al. 2008; Gröschel et al. 2013; Klingner et al. 2010, 2014) and hand movement (Ferbert et al. 1992; Allison et al. 2000; Hamzei et al. 2002; Stefanovic et al. 2004; Newton et al. 2005; Perez and Cohen 2008). These ipsilateral somatosensory evoked deactivations were associated with an elevated sensory threshold (Kastrup et al. 2008) and their amplitude was found to increase with stimulus duration and intensity (Klingner et al. 2010), suggesting an underlying inhibition or dys-facilitation of the nonstimulated hand. However, one limitation of previous studies is that they stimulated one body part (the hands) and tested the positive and negative BOLD responses only in the contralateral and ipsilateral hand representation. Here, by taking a more holistic approach and looking at the pattern of ipsilateral BOLD signals following stimulation of many body parts, we can examine whether this deactivation in the ipsilateral hemisphere is specific to the area of the ipsilateral stimulated body part as suggested by previous studies (Hlushchuk and Hari 2006; Kastrup et al. 2008; Gröschel et al. 2013; Klingner et al. 2010, 2014) or can be detected for all other ipsilateral body parts, thus reflecting a general mechanism for all nonstimulated body parts whether contralateral or ipsilateral.
Materials and Methods
A total of 21 healthy right-handed subjects (8 females) aged 24–37 (mean age 29) with no neurological deficits were scanned in the current study. All participants were scanned in the continuous-periodic somatosensory experiment, with different subsets of subjects participating in the block-design experiment and resting-state scan (specified below). The Tel Aviv Sourasky Medical Center Ethics Committee approved the experimental procedure and written informed consent was obtained from each subject before the scanning procedure.
Experiments and Stimuli
Periodic Experiment (n = 21)
This experiment included 2 scanning sessions of somatosensory stimulation. The natural tactile stimulus was preformed manually using a 4-cm-width paint brush. In each stimulation block, the body surface was stimulated sequentially by brushing the right side of the subjects’ skin surface from the lips, and then continuously from the fingers and palm through the shoulder, waist, knee and down to the toes. This stimulation order was reversed in the second scanning session so each body segment was brushed backwards (i.e., toes to knee, knee to waist). The length of each stimulation cycle was 15 s (from lips to toes or vice versa), which was followed by a 12-s rest baseline (Fig. 1A). Each scanning session included 8 blocks of tactile stimulation. Thirty seconds of silence were added before and after the 8 cycles of tactile stimulation for baseline. The stimulus was passive and subjects were asked to concentrate on the tactile sensation. We chose this continuous stimulus paradigm as it has been shown to be an optimal method for topographical mapping (Sereno et al. 1995; Engel et al. 1997; Engel 2012) that requires relatively short scanning sessions.
In order to further test the pattern of somatotopically organized negative BOLD responses, we conducted a block-design experiment. In each block, the subject's right body was brushed in one of the following locations: lips, fingers and palm, shoulder, waist, knee and foot (in 2 perpendicular directions). A stimulation block lasted 9 s, followed by a 9-s rest baseline. Each body part stimulation was repeated 4 times, in a pseudorandomized order. A subset of 8 subjects participated in this experiment.
In both experiments, the tactile stimulation was delivered by manual brushing of the body. We chose this stimulus method for several reasons. First, since brushing the skin surface imposes a robust tactile stimulus, it is likely to activate both cutaneous and deep tissue receptors in the primary somatosensory area. This stimulation paradigm has been shown to result in a wider activation pattern and to be more effective in recruiting all the subareas of the primary somatosensory cortex, compared with air-puff stimulation (Huang et al. 2012). Second, a broad stimulus has been found to be effective in activating higher somatosensory areas that are characterized by neurons with large receptive fields (Robinson and Burton 1980; Disbrow et al. 2000). Although this stimulation method is somewhat less controlled than an automatic stimulus, and may lead to small variations in stimulus strength and timing, it provides an easy and fairly natural way to stimulate the whole body. Moreover, the relatively slow temporal sampling rate of the magnet entails a temporal resolution threshold in which small deviations from precise timing (on the order of tens to hundreds of milliseconds) will not affect the large-scale somatotopic organization. In order to maintain constant stimulation parameters and minimize potential distortions of the image, we took several steps. The stimulation of all subjects was performed by the same experimenter, who was well trained prior to the scans to maintain a constant pace and pressure during the sessions. Precise timing of stimuli was achieved by auditory cues delivered to the experimenter brushing the subject through fMRI-compatible electrodynamic headphones (MR-Confon). The paramagnetic paint brush was made of wood with a long (60 cm) cane that enabled the experimenter to stimulate all body parts. The experimenter stood next to subject's right body side (the stimulated side), and took extreme measure not to touch the subject or bend into the magnet's bore throughout the whole stimulation cycle. Furthermore, in order to directly test for putative bias in the resulted images caused by the manual stimulation paradigm, we compared the whole-brain distributions of the responses to lips and toe stimulation in both experiments and found no evidence for such a bias (see
Functional and Anatomical MRI Acquisition
The BOLD fMRI measurements were obtained in a whole-body, 3-T Magnetom Trio scanner (Siemens). The scanning session included anatomical and functional imaging. The functional protocols were based on multislice gradient echoplanar imaging (EPI) using Siemens's 12-channels Head Matrix Coil. The functional data were collected under the following timing parameters: time repetition (TR) = 1.5 s, time echo (TE) = 30 ms, flip angle (FA) = 70°, imaging matrix = 80 × 80, field of view (FOV) = 24 × 24 cm (i.e., in-plane resolution of 3 mm). In all, 27–35 slices with slice thickness = 4.1–4.3 mm and 0.4–0.5 mm gap were oriented in the axial position, for complete coverage of the whole cortex. The first 10 images (during the first baseline rest condition) were excluded from the analysis because of nonsteady-state magnetization. High-resolution 3D anatomical volumes were collected using a 3D-turbo field echo (TFE) T1-weighted sequence (equivalent to magnetization prepared rapid gradient echo (MP-RAGE)). Typical parameters were FOV 23 cm (RL) × 23 cm (VD) × 17 cm (AP); Foldover- axis: RL, data matrix: 160 × 160 × 144 zero-filled to 256 in all directions (ca. 1 mm isovoxel native data), TR/TE = 9 ms/6 ms, FA = 8°. Cortical reconstruction included the segmentation of the white matter using a grow-region function embedded in the BrainVoyager QX 2.0.8 (Brain Innovation) software package. The cortical surface was then inflated. Group results were superimposed on a 3D cortical reconstruction of a Talairach normalized brain (Talairach and Tournoux 1988).
Preprocessing of fMRI Data
Data analysis was initially performed using the BrainVoyager QX 2.0.8 software package. Functional MRI data went through several preprocessing steps that included slice scan time correction and head motion correction. No data included in the study showed translational motion exceeding 2 mm in any given axis, or had spike-like motion of >1 mm in any direction. Preprocessing also included high-pass filtering using a general linear model (GLM)-Fourier based approach implemented in the BrainVoyager QX 2.0.8 software package. In this method, the GLM is used to estimate the contribution of low-frequency and linear drifts to the voxels’ time course by a design matrix containing sines and cosines of low frequency (cutoff frequency: 2 cycles/scan), as well as a linear trend predictor. The resulted predicted time course was then subtracted from the original time course to obtain a filtered time course. In the periodic experiment, we also implemented temporal smoothing to improve the signal-to-noise ratio by convolving the fMRI signal in the time domain using a Gaussian kernel with full width at half maximum (FWHM) of 4 s. For group analyses, the functional data also underwent spatial smoothing (spatial Gaussian smoothing, FWHM = 6 mm) to overcome intersubject anatomical variability within and across experiments. Functional and anatomical data sets for each subject were aligned and fit to the standardized Talairach space (Talairach and Tournoux 1988).
In order to map the somatosensory homunculus, we used 2 phase-locking analytic approaches—spectral analysis and cross-correlation. We conducted the 2 pipelines of analysis independently from each other, without using any result from one pipeline for the analysis in the other.
Spectral analyses were conducted with an in-house program using MATLAB (MathWorks). Following standard retinotopy procedures (Engel et al. 1994, 1997), we applied Fourier analysis to the time course of each voxel, locked to the stimulus repetition frequency (Hertz and Amedi 2010; Striem-Amit et al. 2011; Zeharia et al. 2012). The Fourier analysis delivered amplitude and phase values of the time course in each voxel.
Figure 1B shows the Fourier amplitudes in various frequencies in a typical voxel at the postcentral gyrus (PCG). This amplitude was the highest at the repetition frequency (0.037 Hz). For each voxel, both the amplitude and phase parameters were used to construct a pure cosine, which served as a model of the activation:
A Pearson's correlation coefficient was then calculated between the model and the original time course, yielding a correlation coefficient for each voxel.
Figure 1C shows the BOLD signal in the same voxel, along with a pure cosine function made up of the amplitude and the phase values at the repetition frequency. This cosine profile served as a model of the activation for further statistical evaluation of the results. A Pearson's correlation coefficient was calculated between this model cosine and the time course in each voxel and this was used as a direct measure of the voxel's response to the tactile stimulation (regardless of the specific body segment represented in this voxel).
In regions showing high correlation to the stimulus repetition frequency, the phase value was inspected. Phase values were distributed between –π and π, and were linearly transformed to range between 0 and 18, representing time points in each stimulus cycle in TR units. Since the tactile stimulation of a given body segment was always made at a fixed time point in the stimulus cycle, the phase value served as a measure of the relative time in which the body segment was touched. Due to the time delay of the hemodynamic response (Logothetis and Wandell 2004), the phase code does not temporally overlap with the stimulus presentation time. The onset of the first response detected in the PCG was considered to represent the response to the first body part that was stimulated (lips or foot, according to the stimulus session). Similarly, the latest response observed in PCG was assumed to correspond to the last stimulus that was presented. These values formed the phase code corresponding to the specific preferred body area of each voxel, and resulted in individual phase code maps which corresponded to the individual tactile maps.
Phase-locked based paradigms usually use a continuously varies stimulus, which is presented in a constant order in each stimulation cycle (although some group use a related methods in which the stimuli are presented in quasi-random orders; see Hansen et al. 2004 and Vanni et al. 2005). In order to minimize potential order effect and following previous works of retinotopic, cochleotopic and motor mapping (Sereno et al. 1995; Engel et al. 1997; Striem-Amit et al. 2011; Zeharia et al. 2012; Alvarez et al. 2015), we averaged the phase maps obtained from the 2 stimulation direction scans. The data of each stimulation direction were analyzed separately and the 2 scan directions were combined by averaging the resulted maps as follows. First, the phase values of the toe-to-lip maps were flipped such as that the first phase values became the last and the last phase values became the first. Next, the phase values from the individual lip-to-toe maps were averaged with the flipped values from the toe-to-lip maps to yield single-subject averaged maps. The resulting individual maps of 3 individual subjects are shown in
On the group level, the phase maps of the individual subjects were averaged to create a mean phase map. To run a random-effect analysis, we used GLM parameter estimators’ analysis, as follows. First, a GLM analysis was carried out at the single-subject level using the pure cosine model described above as a predictor. The pure cosine model predictor is positively biased, because it is derived from the Fourier analysis. To account for this bias, we applied the same approach to 20 other non stimulus related frequencies for each subject. The average GLM parameter estimator value from these analyses was used as the bias estimator and subtracted from the GLM parameter estimator values in the stimulus representation frequency. The resulting GLM parameter estimator values were then used in a second-level analysis for the group random effect. Finally, the random-effect results were corrected for multiple comparisons using the ClusterThresh plug-in in BrainVoyager QX. This cluster-size threshold correction procedure for multiple comparisons is based on a Monte Carlo simulation approach (Forman et al. 1995), which was extended to 3D data sets (Goebel et al. 2006). This method exploits the fundamental assumption that areas of activity tend to stimulate signal changes over spatially contiguous groups of voxels rather than over sparsely isolated voxels, and it uses randomness to estimate a cluster-extent threshold. The computation of the minimum cluster threshold is accomplished via Monte Carlo simulation of the random process of image generation (1000 iterations in our case), followed by the injection of spatial correlations between neighboring voxels, voxel intensity thresholding and cluster identification. The resulting map becomes “corrected for multiple comparisons” at a desired confidence level (α < 0.05 in our case). The averaged group phase maps were thresholded by both the averaged correlation coefficient and the random effect corrected for multiple comparison maps. The resulting group maps presenting the known somatotopic gradients of S1 and M1 are shown in Figure 1D.
The cross-correlation analysis was carried using the Linear Correlation algorithm implemented in BrainVoyager package. In this procedure, a boxcar function 2 TRs long (the duration of stimulation of one body segment) was convolved with a 2 gamma hemodynamic response function (HRF, Friston et al. 1998); with the following parameters—onset: 0, time to response peak: 5 s, response dispersion: 1, response undershoot ratio: 6, time to undershoot peak: 15, undershoot dispersion: 1—to derive predictors for the analysis. These predictors and the time course of each voxel were cross-correlated, allowing for 12 lags with one TR interval time to account for the stimulation duration of each cycle (10 TRs) and 2 additional TRs to account for potential additional hemodynamic delay. This algorithm calculates for each voxel the correlation coefficients for the 12 predictors, and results with cross-correlation map, that present for each voxel the lag value with the highest correlation coefficient. Similarly to the spectral analysis, the lag value served as a measure of the relative time in which the body segment was touched, and they were color coded from red to blue, to present the somatotopic gradient from lips to toe. In order to allow averaging of the maps from the 2 stimulation directions, the lag values from the toe-to-lips experiment were flipped such as that the first lag values became the last and the last lags values became the first. The cross-correlation maps of the individual subjects were averaged across the 2 scan directions to create a mean cross-correlation map (see
One difference between the cross-correlation and spectral analysis methods is that in spectral analysis the model of activation (a pure cosine in the stimulation frequency) does not include the hemodynamic delay, rather, this delay is reflected in the resulted phase values of each voxel. For example, if a voxel is responsive to lip stimulation its corresponding lag value would be expected to be 1 (after convolving with the hemodynamic function), while its phase value would be expected to be 4 (corresponding to the time of lip stimulation plus a hemodynamic delay, Logothetis and Wandell 2004). Nevertheless, in both methods, the reversal of the lag and phase values and the average of both stimulation directions maintained the topographic gradient and its dynamic range. Furthermore, as the HRF is expected to delay the responses in opposite directions in the 2 scans, another advantage of this average procedure is that it also helps to cancels this phase delays (Sereno et al. 1995). Indeed the results were quite similar in both methods and serve to show that the data are robust regardless of the exact phase-locking approach we used.
In order to assess the negative BOLD responses elicited by specific body segments, we applied a GLM analysis on the periodic and block-design experiments using predictors convoluted with a typical 2 gamma (with similar parameters to the HRF used in the cross-correlation analysis described above). Since we were interested in the somatosensory and motor cortices, the GLM analysis was restricted to these areas by applying a cortical mask to the model. Cross-subject statistical parametric maps were calculated using hierarchical random-effects model analysis (Friston et al. 1999). A statistical threshold criterion of P < 0.05 was set for all results, which were corrected for multiple comparisons using the ClusterThresh plug-in in BrainVoyager as was described in the previous section.
In the periodic experiment, the raw time course of each cycle of responses to the continuous full-body tactile stimuli was divided into 5 brushing segments corresponding to the stimulation of cortically adjacent body parts according to Penfield's classic homunculus (lip, hand, shoulder and upper trunk, waist to knee and knee to toe). This division of the continual stimulus into 5 segments roughly reflected the dermatomal arrangement of the skin surface: lips (trigeminal cranial nerve), arm and trunk (cervical dermatomes) and waist and leg (thoracic and lumbar dermatomes). The responses to brushing the 5 segments were averaged from the 2 sessions with lips-to-toe and toe-to-lips strokes to avoid order and attention effects. Error bars display the standard error across subjects and sessions.
To investigate the phase shifts in the hemodynamic responses as a continuous and to test for negative BOLD responses, we sampled the raw BOLD signal from LH (contralateral to the stimulated body side) S1 and M1 homunculi at sequential points across the group mapping sequence (Engel et al. 1997). The time course of activation from these ROIs was sampled and GLM analysis was conducted in each ROI for the 2 scanning sessions, across all subjects, yielding the GLM parameter estimators of stimulation for each of the body segments versus the rest baseline. In order to apply a similar analysis to the ipsilateral hemisphere, we defined homologous set of S1 ROIs in the right hemisphere (by vertical axis inversion), based on the bilateral symmetries of these homunculi.
Functional Connectivity MRI Acquisition and Data Analysis (n = 13)
A data set of spontaneous BOLD fluctuations for the investigation of intrinsic (rest state; Biswal et al. 1995) functional connectivity was collected while the subjects lay supine in the scanner without any external stimulation or task (data in-plane matrix size, 64 × 64; FOV, 24 × 24 cm; TR = 2000 ms; FA, 90°; TE = 30 ms). To obtain full coverage of the subjects’ brains, 35 slices of 4 mm thickness were used. In all, 160 whole-brain images were collected in one functional scan and the first 2 images of each scan were excluded from the analysis because of nonsteady-state magnetization. Ventricles and white-matter signals were sampled using a grow-region function embedded in the BrainVoyager from a seed in each individual brain. Using MATLAB (MathWorks), ventricle and white-matter time courses were regressed out of the data and the resulting time course was filtered to the frequency band-width of 0.1–0.01 Hz (in which typical spontaneous BOLD fluctuations occur). The resulting data were then imported back into BrainVoyager for group analyses. Single-subject data were spatially smoothed with a 3D 6 mm half-width Gaussian to reduce intersubject anatomical variability. Five ROIs were defined as group peak activations of the five stimulated body segments from the periodic experiments, and served as seeds for whole-brain functional connectivity analysis. Individual time courses from these seed ROIs were sampled from each of the subjects who participated in the rest scans, z-normalized, and used as individual predictors in the whole-brain group connectivity analysis using a GLM with a hierarchical random effect. The minimum significance level of the results was set to P < 0.05 (corrected for multiple comparisons) to inspect the whole-brain connectivity pattern or to P < 0.01 to view the peak connectivity in the ipsilateral hemisphere.
For ROI-based cross-correlation analysis, we also defined homologous set of 5 ROIs in the right hemisphere (by vertical axis inversion). Functional connectivity was obtained by calculating Pearson's linear correlation between the time courses of each pair of ROI, giving a 10 × 10 matrix of correlation coefficients. The Correlation coefficients were Z-transformed, tested for significant using a t-test (for the contrast between the correlation coefficient and zero), and corrected for multiple comparison (Bonferroni correction for P < 0.05). For adjacency analysis, each pair of ROI was assigned with an adjacency index, reflecting the somatotopic proximity between the ROIs (see Fig. 7A for a matrix of all adjacency indexes). Thus, for example, within each hemisphere, adjacency index of one corresponds to consecutive ROIs such as lips and hand, hand and shoulder, shoulder and waist and waist and foot. Adjacency index of 2 corresponds to pairs of ROIs, which are separated by another body segment. These include lips and shoulder, hand and waist and shoulder and knee. Adjacency index of 3 includes the lips and waist and hand and foot pairs, while the lips and foot connection represent index of 4. The same indexing method was applied to the connections between the hemispheres (e.g., the connection between left hemisphere lip ROI and right hemisphere hand ROI was assigned adjacency index of 1). The adjacency index of the connections between homologous ROIs of the 2 hemispheres (e.g., lips and lips, hand and hand) was set to zero. The ROIs’ correlation coefficients were grouped according to their corresponding adjacency indexes to present the connectivity as a function of somatotopic organization.
Using a set of functional MRI experiments, including a unilateral tactile passive stimulation of the whole body, we tested for gradients and topographic biases of both positive and negative BOLD responses in the contralateral and ipsilateral S1. We hypothesized that in addition to the positive responses in the contralateral hemisphere, tactile stimulation of different body parts would result in topographically organized deactivations in areas corresponding to the nonstimulated body parts in both the ipsilateral and contralateral hemispheres. The data were analyzed on several levels and approaches, at both the group and single-subject level: first, to characterize the large-scale somatotopic organization of the tactile responses we implemented 2 phase-locked analysis methods: spectral analysis and cross-correlation. Next, we searched for negative BOLD responses in the contralateral and ipsilateral hemispheres through whole-brain and ROI GLM analysis of the periodic and block-design experiments. Finally, we investigated the connectivity patterns between the somatosensory homunculi, using functional connectivity analysis of resting-state fMRI.
Our first goal was to map the primary somatosensory homunculus following a continuous and periodic tactile stimulation of the right side of the body (Fig. 1A). This periodic design enabled us to apply phase-locked analysis methods, which have been shown to be optimal for detecting areas that contain gradual representations or topographic gradients. Two phase-locked analysis methods were used. First, we applied spectral analysis, in which phase values were inspected in regions showing a high correlation with the stimulus repetition frequency, and served as a measure of the relative time in which the body segment was stimulated (see Materials and Methods and Fig. 1B,C). Second, we applied cross-correlation analysis which yielded a map of lags and correlation values that was also used to inspect the relative representation of the body segments throughout the stimulation cycle. The resulting group maps extracted by these methods are presented in a coronal view (Fig. 1D), and on the inflated cortical surface (Fig. 1E,F). Both methods revealed the detailed somatotopic organization of the primary somatosensory cortex contralateral to the stimulated body side. The whole-body detailed somatotopic maps were also evident at the single-subject level using both phase-locked methods on spatially smooth and nonsmooth data sets (1F, right).
A completely different pattern was revealed in the ipsilateral hemisphere. Although most of the tactile evoked responses were localized in the posterior parietal lobe, posterior to the postcentral sulcus, there was no phase-locked response in the ipsilateral S1. These findings are consistent with previous studies showing that ipsilateral activations are localized posterior to the contralateral ones (Iwamura 2000; Nihashi et al. 2005). The only exception was the response to lip stimulation, which was visible in the ventral ipsilateral cortex. This result was expected, given that this was the only body area in which the stimulus was bilateral.
Next, we tested for negative BOLD responses in the contralateral hemisphere. Since both phase-locked analysis methods use the highest correlation value between the model and the voxel's time course as a measure of the relative time in which the body segment was stimulated, they are not appropriate for detecting simultaneous negative evoked responses. Thus, to do so, we defined 7 consecutive ROIs, from the most ventral lip representation to the most dorsal toe representation in the contralateral S1 (Fig. 2A). The BOLD signal from each ROI was sampled and a group GLM analysis was calculated for the body segments. Figure 2B,C presents the extracted GLM parameter estimator and the average time course for the separate body segments in S1 ROIs. The results also served to test the gradual shift in the representation of different body parts along S1, as depicted by the peak parameter estimator and time course. Importantly, this analysis can also test for significant deactivations in the primary somatosensory homunculus (Zeharia et al. 2012). Specifically, in most of the ROIs along the PCG we found a combined pattern of significant positive and negative BOLD responses to different body parts. For example, whereas in the most ventral ROI within the PCG the lips and hand evoked positive BOLD responses, the waist and foot evoked negative BOLD. This pattern was reversed in the most dorsal ROI of the PCG, where the foot stimulation resulted in positive BOLD and the lips and hand evoked deactivation. A similar trend was also observed in the contralateral M1 (
In the next step, we explored the tactile evoked negative responses in the ipsilateral hemisphere. We showed above that the ipsilateral S1 was not activated using the phase-locked analysis methods (Fig. 1). In order to apply a similar ROI GLM analysis approach, we defined homologous ROIs in the ipsilateral hemisphere using the same 7 ROIs (see Materials and Methods). These ROIs served for the GLM parameter estimation analysis. Figure 3 presents these results for the separate body segments versus the baseline as well as the contrast of the summation of all body segments versus baseline. Both types of analysis revealed several interesting and novel findings. First, when contrasting all body segments with the baseline, a general deactivation was observed in most of the ROIs of the ipsilateral PCG. Second, inspection of the separate responses for the different body segments showed that the evoked negative BOLD responses in the ipsilateral hemisphere were not confined to the homologous area within S1 (as previously suggested), but rather extended to other parts of the homunculi. For example, stimulation of the hand evoked significant negative BOLD in medial ROIs (ROIs 4 and 5 in the PCG), corresponding to hand representation, but also in dorsal ROIs (ROIs 6 and 7 in the PCG), corresponding to waist and foot representations. Furthermore, foot evoked deactivations were evident almost in the entire ipsilateral S1, though they did not reach significance in the most dorsal ROIs. There were 2 exceptions to these results. One is trivial—the lips were the only area in which the stimulation was bilateral and therefore there was a significant activation for the lips in the most ventral ROI, whereas significant deactivations were observed in most of the other ROIs. Second, the only body segment that showed a trend for positive parameter estimator rather than negative (though not significant above baseline) in response to unilateral stimulation was the shoulder and upper trunk. This result is consistent with previous studies showing that midline structures in the axis of the body (such as the trunk) have symmetrical bilateral representations (Fabri et al. 2005; Eickhoff et al. 2008). To summarize, we showed that beside the axial bilateral structures, in all other ROIs not only was the negative BOLD not limited to the ipsilateral stimulated segment, but usually it was not even the peak of the negative BOLD (in term of the beta GLM parameter estimator).
To further characterize the positive and negative BOLD across the 2 hemispheres and to verify the results using another method, we conducted a GLM analysis of the BOLD signal across the entire contralateral and ipsilateral sensory-motor cortices (this was done on a voxel-by-voxel basis unlike the ROI approach above; see Materials and Methods). The general pattern of activation and deactivation in response to lips, hand and foot stimulation is presented in Figure 4. Because negative BOLD was defined as a decrease in the BOLD signal relative to the rest baseline, the negative responses were seen when contrasting stimulation periods with rest. The results revealed a somatotopic organization of the negative BOLD in S1 (Fig. 4A, t(20) = 2.2, P < 0.05, corrected for multiple comparisons). In the contralateral hemisphere, where lip stimulation elicited positive BOLD, hand and foot elicited negative BOLD; where hand stimulation evoked positive BOLD, foot evoked negative BOLD, and where foot stimulation elicited positive BOLD, the lips and hand elicited negative BOLD. In the ipsilateral hemisphere, positive BOLD responses were observed in the most posterior parts of S1 and the posterior parietal cortex whereas the negative BOLD responses were localized mostly to the anterior part S1, along the posterior bank of the central sulcus. The pattern of the negative BOLD in the ipsilateral hemisphere was more extensive compared with the contralateral side, and also included the homologous area of the unilateral stimulated body segment. Thus, the areas in the ipsilateral hemisphere corresponding to the hand and foot were deactivated by other ipsilateral body segments and not only by the hand or foot. The inverse relations of the positive and negative BOLD in S1 are clearly visible when the statistical parametric maps are plotted separately for positive or negative BOLD (Fig. 4B and see4B), reflecting the widespread characteristics of the ipsilateral BOLD.
Finally, to further investigate the connectivity patterns of the bilateral somatosensory homunculi, and to assess the potential connections mediating the bilateral positive and negative responses, we conducted a functional connectivity analysis of resting-state fMRI data. This was done at 2 levels—whole-brain connectivity analysis and ROI-based connectivity between pairs of ROI in S1. For whole-brain analysis, we defined 5 seed ROIs in the contralateral hemisphere representing 5 segments of the body (see Materials and Methods). These ROIs served as seeds for functional connectivity maps of the different body segments (Fig. 5, t(12) = 4.2, P < 0.01, corrected, showing the connectivity of 3 seeds; see
Next, to further characterize these patterns we calculated the ROI-based correlation of the resting-state BOLD signal within and between the homunculi. For that we defined additional 5 homologous ROIs in the right hemisphere. The connectivity matrix of all the connections shows that while all S1 ROIs are positively connected, both within and between hemispheres (Fig. 6A), the connections between the most distant ROIs (lips and foot in both hemispheres and lips with RH-waist) do not reach significance level following correction for multiple comparisons (Fig. 6B). Figure 6C summarizes these connectivity patterns by illustrating a separate diagram of connections for each pair of homologous ROIs. The results show that for most ROIs there is a gradual decrease in connectivity with farther ROIs. For example, within the left hemisphere, the hand is highly connected to the most adjacent ROIs—lips and shoulder (R = 0.53 and 0.62, respectively) less to the waist (R = 0.48) and least for the foot (R = 0.36).
In order to quantify these patterns and test whether the connectivity reflects topographic organization, we looked at the functional connectivity as a function of ROIs’ adjacency along S1. The correlation coefficients of the different ROIs were grouped according to their adjacency index (see Materials and methods for details and Fig. 7A for adjacency matrix) in several levels—pooled over both hemispheres, within and between hemispheres (Fig. 7B–D). The results clearly show that the connections’ strength decreases with increasing adjacency index. The strongest connections are between the somatotopically homologous ROIs (adjacency index zero, lips–lips, hand–hand est.), and then gradually decline until they fail to reach significance level for the far most connections (adjacency index 4, lips–foot). Thus, the resting-state functional connectivity both within and between the homunculi is somatotopically organized.
In this work, we studied the somatotopic organization of the positive and negative BOLD responses in contralateral and ipsilateral primary somatosensory homunculi and found several novel results along with several verifications and extensions of the previous literature. Using continuous and periodic unilateral tactile stimulation of the entire body and applying phase-locked analysis methods, we were able to verify the well-known contralateral topographical organization of S1, and to demonstrate the gradual shift in the representation of different body parts (Figs 1 and 2 and
Mapping the primary somatosensory homunculus
Since its discovery by Penfield and Boldrey (1937), the principal organization of the primary somatosensory cortex has been verified in large number of species, using variety of methods (Kaas et al. 1979; Nelson et al. 1980; Kaas 1983; Shoham and Grinvald 2001; Rothemund et al. 2002; Chen et al. 2005, 2007). The somatotopic organization of S1 in humans has been partially confirmed noninvasively by neuroimaging studies (Kurth et al. 1998; McGlone et al. 2002; Kell et al. 2005; Miyamoto et al. 2006; Huang and Sereno 2007; Nelson and Chen 2008; Schweizer et al. 2008; Overduin and Servos 2008b). However, the vast majority of these studies focused on distinct body part stimulation. Hence, a detailed full-body somatotopic organization in humans has not been fully investigated using imaging techniques. Moreover, using only one or a few body parts makes it impossible to use a phase-locking approach to map whole-body gradients; this was the first goal of this study. Here, using cross-correlation and spectral analysis we mapped Penfield's homunculus in great detail. These maps revealed the phase shifts of activation from the lip activated area in the lateral surface of the PCG to the foot and toe activated area at the most dorsal part of the PCG and hemispheral rim (Fig. 1 for group results; Fig. 2 for phase shift in the time course;
In this study, we used 2 phase-locking analytic approaches for somatotopic mapping: cross-correlation and spectral analysis. These approaches are considered the classical means of defining early visual areas (Sereno et al. 1995; Engel et al. 1994, 1997; Wandell and Winawer 2011) and have also been used to map other topographically organized sensory (Hertz and Amedi 2010; Orlov et al. 2010; Striem-Amit et al. 2011) and motor (Zeharia et al. 2012, 2015) gradients. They yield a highly robust topographic mapping and can pinpoint gradients undetected by standard GLM analysis (Engel 2012; Zeharia et al. 2015). In the somatosensory system, these approaches have been used to map detailed local gradients of the arm (Servos et al. 1998), fingers (Overduin and Servos 2008b; Sanchez-Panchuelo et al. 2010; Mancini et al. 2012; Besle et al. 2013), or in higher somatosensory face selective area in the VIP cortex (Sereno and Huang 2006). Moreover, in contrast to most of the previous studies which used air-puff stimulation, we applied a more natural tactile stimulation of passive brushing of the skin. This stimulation paradigm has been shown to result in a wider activation pattern and to be more effective in the recruitment of the primary somatosensory cortex than air-puff stimulation (as also observed by Huang et al. 2012). However, it is important to note that despite the power of these phase-locked methods in noninvasive mapping of topographic gradients, they are not the optimal method to infer the exact locus of representation of each body part along the topographic gradient (although this kind of localization is possible, see Besle et al. 2013). Such a validation of the results could be achieved using other methods such as GLM analysis, which enables a direct contrast between the responses to different body parts.
Our results demonstrate that a natural tactile stimulus combined with phase-locked analysis provide a powerful and sensitive tool for the mapping of whole-body somatosensory gradients at the group and single-subject levels. These experimental paradigm and analytical approaches should in our view be the current golden standard to map whole-body gradients in both basic research and future clinical setting. For example, as is done in retinotopic mapping, somatotopic gradients can be used to define brain areas and their boundaries, the connectivity between different somatosensory areas as well as in plasticity studies. Furthermore, beyond improved resolution and clearer gradients, this periodic experimental design requires a shorter scanning time and fewer repetitions compared with other experimental designs, and thus has some clear practical benefits especially in potential clinical setting (e.g., mapping somatosensory areas in pre-operative procedures, guiding deep brain stimulation, Yamgoue et al. 2016). A highly interesting future application of this paradigm might be in the rapidly growing and exciting field of somatosensory restoration. Our paradigm for example can be a useful tool in studying the neural correlates and plasticity processes evoked by artificial somatosensory feedback such as somatosensory prosthetics or intracortical microstimulation of the somatosensory cortex (Tabot et al. 2015; Tyler 2015; Delhaye et al. 2016) and might help in advancing these as a procedure to augment and restore somatosensation.
Negative BOLD responses in the contralateral homunculus
Negative BOLD responses have been studied most extensively in the visual system of both primates and humans (Tootell et al. 1998; Shmuel et al. 2002; Smith et al. 2004; Bressler et al. 2007; Pasley et al. 2007; Wade and Rowland 2010). These studies have shown that the negative BOLD responses in the early visual cortex can be elicited in areas surrounding the positively activated areas. Specifically, stimulation of central visual field evoked negative BOLD responses, which were coupled to decreased local field potentials and multiunit activity in the surrounding peripheral areas (Shmuel et al. 2006). Such a precise pattern of deactivation has been shown to contain stimulus-specific distributed information and was suggested to play an important role in the resulting visual percept (Bressler et al. 2007). Our results document a similar phenomenon following tactile stimulation of the body. For the first time, we show topographically organized negative BOLD responses in the areas of the nonstimulated body parts along the primary somatosensory and motor homunculi. The current literature on negative BOLD responses in the somatosensory system has focused on unilateral stimulation of the hand, which was shown to result in deactivation of the homologous area in the ipsilateral hemisphere (Hlushchuk and Hari 2006; Kastrup et al. 2008; Schäfer et al. 2012; Gröschel et al. 2013; Klingner et al. 2014). We extend these findings and show that negative BOLD responses are a defining characteristic of the contralateral Penfield homunculus (Figs 2 and 4,Zeharia et al. 2012), but not in any of the other recently discovered motor homunculi (Zeharia et al. 2015). A similar trend in M1 was also evident in our results, though both positive and negative BOLD signals were lower and less significant (
Taken together, the results from the visual and motor domains as well as the current somatosensory study indicate that negative BOLD responses are manifested in primary sensory-motor areas, and suggest that they are an important component of the BOLD signal that reflects a neuronal sharpening mechanism of the signal-to-noise ratio. These topographic negative responses seem to specifically characterize primary cortical areas, as higher order sensory and motor fields contain neurons with larger and overlapping receptive fields. It would be interesting to test whether additional topographic sensory areas such as the primary auditory cortex contain the same sharpening mechanism. In fact, though this was not tested directly, hints were found in a previous work by our group which identified multiple cochleotopic maps in the human auditory cortex (Striem-Amit et al. 2011). Specifically, closer inspection of the time courses from different points along the auditory gradients (see Figure 6 in Striem-Amit et al. 2011) suggests that such topographic negative BOLD responses are evident also in the auditory cortex.
Negative BOLD responses in the ipsilateral hemisphere
Here, by using unilateral body stimulation and inspecting the pattern of ipsilateral BOLD signal following stimulation of different body parts, we also extend the results in the literature on deactivation in the ipsilateral hemisphere. This allowed us to further compare the patterns of both hemispheres, since one limitation of previous works by our group in the motor domain (Zeharia et al. 2012) was that the movements were bilateral and could not address this issue. We found that the negative BOLD responses were not limited and did not peak in the ipsilateral homologous area of the stimulated body segment. Rather, the ipsilateral deactivations extended to areas of the nonstimulated body parts, resulting in an overlap of the negative BOLD in the ipsilateral hemisphere across large portions of the homunculus (Figs 3 and 4,Hlushchuk and Hari 2006; Kastrup et al. 2008; Schäfer et al. 2012; Gröschel et al. 2013; Klingner et al. 2010, 2014). However, a closer inspection of the published data shows that in many cases the pattern of ipsilateral deactivations appears to have been more extensive than the representation of the hand (see for example Figures 2 and 3 in Kastrup et al. 2008 and Gröschel et al. 2013), which is consistent with the present results. Since these studies did not include stimulation of other body parts, the precise pattern of the deactivation could not be tested directly. Another potential explanation for the some level of discrepancy of the results could be related to differences in the experimental setup. For example, previous works used median nerve stimulation with various intensities and time intervals (Kastrup et al. 2008; Klingner et al. 2010; Schäfer et al. 2012), which were associated with alteration in the strength of the negative BOLD response (Klingner et al. 2010). In contrast, we applied a more ecological stimulus of brushing the body surface, which was found to evoke robust positive BOLD responses in all the subareas of the somatosensory cortex (Huang et al. 2012) and therefore, could also result in more pronounced negative responses as described here. The fact that these widespread ipsilateral patterns were also observed in the block-design paradigm of single subjects (
The physiological source and functional role of the negative BOLD
Even though the precise physiological mechanism underlying negative BOLD responses is still not fully understood, increasing evidence suggests that these responses reflect a functional measure of neuronal deactivation (Ferbert et al. 1992; Hamzei et al. 2002; Bressler et al. 2007; Kastrup et al. 2008; Klingner et al. 2010; Schäfer et al. 2012; Mullinger et al. 2014). Some of this evidence comes from studies of the somatosensory system, which provides a good model to study the interaction between the BOLD signal and physiological parameters of blood flow, blood volume, metabolic rate and neuronal activity. First, the negative BOLD responses found in the ipsilateral hemisphere (Figs 3 and 4,Harel et al. 2002; Shmuel et al. 2002; Devor et al. 2007) as both cerebral hemispheres belong to different vascular territories (Nair 2005). Our results evidencing somatotopically distributed negative BOLD (together with the recent findings from the motor domain; Zeharia et al. 2012) show that this hemodynamic blood-stealing phenomenon also cannot account for the negative responses within each hemisphere (the contralateral negative homunculus), since 2 different blood vessels supply the ventral-lateral and dorsal-medial parts of the homunculus (Kiernan and Barr 2009). Second, direct intracortical recordings have demonstrated an inhibition of neuronal activity in S1 in response to ipsilateral stimulation (Lipton et al. 2006; Boorman et al. 2010). In humans, simultaneous BOLD, electroencephalography and cerebral blood flow recording further support a neuronal mechanism underling the negative BOLD responses (Mullinger et al. 2014).
The functional role and the precise pathway mediating the negative BOLD response to unilateral somatosensory stimulation are still largely unclear. While all subareas of the primary somatosensory cortex are bidirectionally connected to the thalamus, which could exert an excitatory or inhibitory influence on specific brain regions (Sherman 2007), intracortical and intercortical connections seem the most likely pathways. Although interhemispheric connections were considered to take place mainly between bilateral S2 (Picard et al. 1990; Sutherland 2006), recent studies indicate that such inhibitory connections could arise as well between the primary somatosensory homunculi (Iwamura 2000; Ragert et al. 2011). Specifically, there is evidence of interhemispheric transfer between homologous areas of S1, mainly BA2s (Iwamura et al. 2001; Hlushchuk and Hari 2006; Klingner et al. 2011), which in turn, are densely connected to BA1/BA3b. Unilateral tactile stimulation of the hand activates the central and posterior regions of the corpus callosum (Fabri et al. 2014), and resection of the posterior corpus callosum in humans eliminates the ipsilateral responses (Fabri et al. 1999).
Resting-state functional connectivity patterns characterized by strong bilateral consistency and somatotopic organization
A growing body of evidence suggests that resting-state functional connectivity patterns can be a powerful marker of the baseline state of a system (see Harmelech and Malach 2013 for review). These patterns have been shown to correlate between areas that are parts of the same functional network (Smith et al. 2009; see Fox and Raichle 2007 for review), closely mimicking anatomical (although not necessarily monosynaptic) connectivity. In both motor and visual systems, resting-state functional connectivity patterns characteristically show strong bilateral consistency, resulting in highly symmetrical and topographic networks (Biswal et al. 1995; van den Heuvel and Hulshoff Pol 2010; Yeo et al. 2011; Zeharia et al. 2015; Dawson et al. 2016). Here, we show the same phenomenon in the somatosensory system. Using whole-brain analysis, we found that the peak functional connectivity from different body segments in the contralateral hemisphere was somatotopically localized to the homological areas in ipsilateral S1 (Fig. 5). This finding is also supported by the ROI-based analysis (Figs 6 and 7B), which shows that the strongest connections between the hemispheres are between homologous ROIs, with decreasing connectivity to adjacent body segments. Thus, the topographical gradient of ipsilateral S1 could be reconstructed by the connectivity patterns from seeds in the contralateral hemisphere. Although these findings do not necessarily point to direct interhemispheric connections, they support a common functional connectivity, which might imply interhemispheric integration.
Furthermore, comparing the group evoked negative BOLD responses in the ipsilateral hemisphere (Fig. 4) with the resting-state connectivity maps (Fig. 5) reveals that both tend to peak around the central sulcus, probably corresponding to areas BA3/BA1 in S1. Though the spatial resolution of the data presented here were not sufficient to specifically draw the borders between S1 subareas, this putative overlap between the negative responses and the connectivity pattern stress the dynamic change in tactile processing across the 2 hemispheres, which ranges from symmetrical streams of inputs (i.e., during bimanual tactile exploration) and lateralized activation patterns evoked by unilateral sensation.
The functional connectivity within hemispheres was also found to be topographically organized, as the correlation coefficient decreased with increased separation between ROIs (Figs 6 and 7C). The connections between ROIs representing the extremities of the homunculus (lips and foot—adjacency index of 4) did not reach significance level in both hemispheres. These findings are in agreement with previous works studying functional connectivity in the sensorimotor system (Yeo et al. 2011; Long et al. 2014). Long and colleagues (2014) reported that resting-state fMRI-based parcellation of the sensorimotor cortex identified somatotopically arranged clusters, corresponding to the representation of upper–middle–lower limb segmentation. Similarly to our finding, they report that the foot–tongue connection presented no significant correlation. The lack of significant connections between these body segments could be explained by the spatial distance between these regions, and might also be explained by the fact that these body segments are rarely behaviorally coupled. This explanation is supported by recent studies that stressed the contribution of prior coactivation to the functional connectivity of cortical networks at different scales (Harmelech et al. 2013; Vidyasagar et al. 2014).
Our results of decreased connectivity with increased separation between ROIs parallel findings of retinotopically organized functional connectivity in the visual system. In a recent study, Dawson and colleagues (2016) investigated the organization of resting-state functional connectivity within and between lower visual areas. Within a visual area, they found a consistent decrease in correlation with increasing eccentricity separation in V1 and V2 (and to a lesser extent also in V3). They also showed that the connectivity between visual areas is retinotopically organized as pairs of ROIs with similar eccentricity showed higher correlation than pairs of ROIs distant in eccentricity. Since in this study we focused on the primary somatosensory cortex, such investigation of the functional connectivity between different somatosensory areas was beyond our scope. However, future studies aiming to explore the somatotopic organization and functional connectivity pattern of additional somatosensory areas might contribute to our understanding of the basic organization principles in these areas.
In this work, we studied the somatotopic organization of the positive and negative BOLD responses in contralateral and ipsilateral S1. We verified the Penfield homunculus in detail and found that the homunculi are characterized by a combination of positive and negative BOLD patterns, which differ in the contralateral and ipsilateral hemispheres. Our results indicate a complex pattern of baseline and activity-dependent responses in the contralateral and ipsilateral sides, in which negative BOLD responses characterize both primary sensory-motor areas. This suggest that the negative BOLD is an important component of the hemodynamic response, reflecting a basic mechanism that underlies the sharpening of tuning curves of large populations of neurons, as was previously demonstrated for topographic gradients in both the visual and motor cortices.
An European Research Council grant (grant number 310809 to A.A.); The James S. McDonnell Foundation scholar award (grant number 220020284 to A.A.); The Edmond and Lily Safra Center for Brain Sciences (ELSC) Vision center grant; the Maratier family; The Charitable Gatsby Foundation; the Levzion Program for PhD students (to Z.T.) and by the Hebrew University Hoffman Leadership and Responsibility Fellowship Program (to R.G).
Conflict of Interest: None declared.