Symmetric Sensorimotor Somatotopy

Background Functional imaging has recently been used to investigate detailed somatosensory organization in human cortex. Such studies frequently assume that human cortical areas are only identifiable insofar as they resemble those measured invasively in monkeys. This is true despite the electrophysiological basis of the latter recordings, which are typically extracellular recordings of action potentials from a restricted sample of cells. Methodology/Principal Findings Using high-resolution functional magnetic resonance imaging in human subjects, we found a widely distributed cortical response in both primary somatosensory and motor cortex upon pneumatic stimulation of the hairless surface of the thumb, index and ring fingers. Though not organized in a discrete somatotopic fashion, the population activity in response to thumb and index finger stimulation indicated a disproportionate response to fingertip stimulation, and one that was modulated by stimulation direction. Furthermore, the activation was structured with a line of symmetry through the central sulcus reflecting inputs both to primary somatosensory cortex and, precentrally, to primary motor cortex. Conclusions/Significance In considering functional activation that is not somatotopically or anatomically restricted as in monkey electrophysiology studies, our methodology reveals finger-related activation that is not organized in a simple somatotopic manner but is nevertheless as structured as it is widespread. Our findings suggest a striking functional mirroring in cortical areas conventionally ascribed either an input or an output somatotopic function.


INTRODUCTION
Receptors in the skin provide cutaneous information to the primary somatosensory (SI) cortex. Following detailed electrophysiological studies of cortical somatosensory responses in monkeys, investigations of somatosensory organization in human cortex have either relied on superficial recordings from epileptic patients [1], or more recently from fMRI (functional magnetic resonance imaging) [2][3][4][5]. Previous human somatotopy studies using fMRI have been limited not only by scanning resolution [2,4] but by an underlying premise that fMRI reveals discrete activation ''hotspots'' in response to cutaneous stimulation [3,5], like those measured electrophysiologically as the firing output of cells. Here we sought cortical somatosensory responses using highresolution (4T) fMRI and an interpretation of the measured activity as reflecting more the metabolic demands of dendritic input to cells [6].
Our experimental paradigm adapted the phase analysis method used previously to identify retinotopic maps in visual cortex [7]. A reference time course consisting of periods of stimulation interspersed with non-stimulation intervals was created by a ''sliding window'' [8,9] of stimulation that cycled repeatedly over the digit surface (Fig. 1). Any location on the digit surface experienced alternating windows of stochastic stimulation and silence. Neighboring locations were stimulated in like fashion, but with a stimulation time course that was staggered in time. This design was not strictly a block design, except with respect to individual locations on the digit surface. Any correlated cortical activity should have been driven by an on/off stimulation time course, but one whose particular phase lag would indicate its sensitivity to a certain location on the digit surface. The relatively slow progression of the stimulation window over the digit, and the random nature of the jet activations within any stimulation window, mitigated any ''cutaneous rabbit'' illusions of movement [10]. As a within-digit control, we varied the general direction of stimulation, from the base of the finger to the fingertip, or from the tip to the base.
We passively stimulated the hairless surface of the thumb (D1), index (D2) and ring finger (D4) of human subjects using puffs of air delivered through a custom-built apparatus. While stimulating the digits we measured the cortical blood-oxygenation level dependent (BOLD) response using fMRI. We assigned each functional voxel a particular location on the digit surface based on its phase lag of maximal correlation to the stimulation time course.

Subjects
We scanned seven healthy right-handed human subjects (four males; average age 26 years), at the Robarts Research Institute (London, ON) using a 4T magnetic resonance imager (Varian, Palo Alto, CA; Siemens, Erlangen, Germany). Incorrect slice plane positioning excluded the experiments of one subject, and a significant motion artifact excluded the D4 experiments of subject 1.

Scanning parameters
We used a custom 14-cm quadrature surface coil centered on the left frontoparietal region to acquire functional and anatomical images, and a custom transmit-receive cylindrical birdcage coil to acquire a full-brain anatomical volume in a separate scanning session. A T2*-weighted gradient echo-echo planar imaging (EPI) pulse sequence was used in functional image acquisition (TR = 750 ms, 4 shots, TE = 15 ms, FA = 40u, FOV = 19.2619.2 cm). Nine contiguous pseudocoronal 5-mm thick functional slice planes (1.561.5 mm resolution in-plane) were sampled, and 72 such volumes were collected over each 216-s experiment. A T1-weighted volume was acquired in the same session, consisting of 64 contiguous pseudocoronal 1-mm thick anatomical slice planes (0.7560.75 mm resolution in-plane; TR = 12 ms, TI = 500 ms, TE = 6 ms, FA = 11u, FOV = 19.2619.2 cm) approximating the orientation of each subject's central sulcus. The full-brain T1-weighted volume included 256 contiguous axial 0.94-mm thick slice planes (TR = 12 ms, TI = 500 ms, TE = 6 ms, FA = 11u, FOV = 24.0624.0 cm).

Stimulation equipment
We used a custom Plexiglas frame to deliver pulses of air via an array of shallow 0.5-cm wide depressions designed to conduct reflected air away from the hand surface. Connected to 16 of these depressions were flexible tubes leading to a series of manifoldmounted air valves. The pattern of these 16 connections was chosen for each subject to optimize the distribution of the active jets under the digits. In particular, the jet locations were spaced regularly about 1.5 cm apart in a straight line under D1 (4 jets), D2 (6) and D4 (6). The tubing length and aperture was consistent across channels. The air valves were computer-controlled and received air from a compressor at a steady-state pressure of 240-275 kPa. Subjects' eyes were closed; a bite bar stabilized the head.

Experimental paradigm
In six experiments on each subject we separately stimulated three digits (D1, D2, and D4) and two directions of stimulation (tip-tobase and base-to-tip). Each subject participated in all experiments in a pseudorandom sequence such that no two successive experiments would involve the same digit. Each experiment was defined as six successive 36-s stimulation cycles. Within each stimulation cycle a ''window'' of stimulation (equal to half of the available jets) rotated through all the jets spanning a digit. During each 18-s half-cycle the potentiated jets were activated in a pseudorandom sequence, in bursts of between five and ten 40-ms air puffs each separated by 60 ms. Each such burst was confined to a single jet, but could be followed by another burst at the same location or at another jet potentiated within the window. These bursts were separated by random inter-burst intervals of 1-2 s. The resulting puff and burst stimulation frequencies were within the range of superficial cutaneous receptors [11].

Anatomical analysis
Our anatomical regions of interest included areas 4, 3a, 3b, and 1. These were defined for each subject based on visual inspection of full-brain cortical anatomy. Areas 4 (along the anterior wall of the central sulcus and onto the precentral gyrus) and 3b (on the posterior wall of the central sulcus) lay facing each other across the fundus of the central sulcus (area 3a), and converged medially at the paracentral lobule. Area 1 lay on the crown of the postcentral gyrus and extended posteriorly to the rostral lip of the postcentral sulcus.

Functional analysis
Prior to our correlation analysis, we preprocessed the functional images in BrainVoyager (Brain Innovation, 2000) within the frequency domain by removing linear trends in vascular activity and highpass filtering at 0.014 Hz (cf. [12]). Functional data were not spatially averaged. The 18-s on/off half-cycles defined a reference time course that we blurred and shifted by a 5-s hemodynamic response function. The reference time course was then shifted by iterative 3-s increments, and correlated with the functional data at each delay. Functional voxels were color-coded according to the phase lag giving maximal correlation (Fig. 1B). Note that we used twelve phase values, each corresponding to a particular 3-s shift of the 36-s stimulation cycle. This number of phase values exceeded the number of stimulation locations on the digits (four or six). We chose to use twelve phase values in order to retain the 3-s temporal resolution of the BOLD signal, and because it was the lowest common denominator of four and six and thus allowed us to use the same phase analysis parameters and color scale for all digits. After coregistering a subject's functional and anatomical slice planes, phase-coded functional data were assigned to one of the regions of interest (ROIs). Voxels above a correlation threshold (r$0.3) were tested for uniformity of phase using Rayleigh's statistic [13].

Cortical flattening
In order to quantify our qualitative impression that functional activity in the precentral gyrus tended to mirror that of the postcentral gyrus, we flattened the pericentral cortical area into a 2D manifold using the Isomap algorithm [14,15; see also 16], a nonlinear dimensionality reduction algorithm. Isomap is a variant of classical multidimensional scaling (MDS [17]), and like MDS, it computes a non-sparse distance matrix for all points in the workspace to find a lower-dimensional representation of the data. Unlike MDS, the distances are not Euclidean but are the sums of the geodesic distances, which accumulate in transition from the center of one local neighborhood to the next. Isomap thus allows a highly-folded surface to be flattened even when voxels on either side of a sulcus may be nearly adjacent in Euclidean space. Isomap also does not assume that the extracted surface is an intrinsically linear 2D sheet, but instead generalizes to a larger class of nonlinear manifolds [14], which likely include the cortical surface or patches thereof [16]. For our purposes we defined the neighborhood of a voxel as its 10 nearest neighboring voxels. This neighborhood parameter e was the minimal size that explained 99% of the variance in geodesic distance estimates, averaged across subjects [15]. Although distances along the 2D manifold are in arbitrary units and cannot be directly equated with Euclidean distances, for plotting purposes we have labeled the approximate rostrocaudal and lateromedial orientations after rotating the manifolds so that the fundus (here, a line joining the medial and lateral termini of area 3a) was vertical, primary motor cortex (MI) was shown on its left and SI on its right. By visual observation of the extracted manifolds we confirmed that the relative positions of the cortical areas were preserved.

Symmetry analysis
Our analysis involved the following steps. 1) We plotted the functional data (above the r$0.3 threshold) over the flattened map of the pericentral cortex, and binned the data into grid cells (arbitrarily set at 20 units 2 ) tiling this surface. Within each grid cell, we took the average phase delay value of all functional voxels within the bin; if none, the grid cell was considered not ''active''. 2) We computed the degree of phase similarity in each pair of active MI (area 4) and SI (3b/1) grid cells that lay an equal distance along the fundus and an equal distance either rostral or caudal to the fundus. Given our uncertainty in assigning area 3a voxels to pre-or postcentral cortex, and in the relationship of this area to MI and SI (see Discussion), we have excluded it from this calculation. The similarity values were normalized such that if two grid cells of a pair had identical phase delays, their phase similarity was 1; if the activity was in anti-phase (e.g. 0u and 180u), the similarity was 21; if phase delays were 90u and 180u, the similarity was 0.5; etc. We then computed an average symmetry score for each experiment, over all pairs of active grid cells in that experiment. 3) To determine a threshold for ''significant'' symmetry, we ran 1000 Monte Carlo simulations for each experiment. In each simulation, we randomly scrambled the grid cells within the pericentral map before performing step 2. We took the 95 th percentile of the resulting distribution of 1000 simulated symmetry scores as the threshold for declaring a true symmetry score obtained from that experiment significant.

RESULTS
Our fMR imaging and subsequent analysis were directed at primary sensorimotor cortex contralateral to the right hand. Using the full-brain anatomical volume from each subject we highlighted the lateromedial extent of several cortical regions, including Brodmann's areas 3a, 3b and 1 (SI), and area 4 (MI). We found in most experiments that the entire ROI was dominated by voxels having only a few of the possible phase delay values. In particular, the distributions of voxels that were correlated to the reference time course in each digit6direction stimulation condition were significantly nonuniform for both D1 and D2 experiments (p,0.01, using Rayleigh's statistic [13] for the circular distribution of r$0.3 voxels), but only marginally for the D4 experiments (p = 0.03). Remarkably, this trend of nonuniform phase response upon D1 and D2 stimulation was true not only in each of SI areas 3b and 1 but also in MI and area 3a in between. In cortex outside of these areas, the population of correlated voxels in these experiments was not significantly tuned (p.0.05).
The D1/D2 phase value distributions appeared to peak near the delay corresponding to stimulation of the digit tips (Fig. 2). Moreover, while the peaks of the tip-to-base and base-to-tip distributions were both within the phase delay range corresponding to fingertip stimulation, they were also significantly different from each other (p,0.01). Recall that each fingertip was stimulated at two locations (Fig. 1C). When the cyclical stimulation window proceeded in a tip-to-base fashion-i.e. when it contacted  the most distal fingertip location before more proximal digit locations-the cortical sensorimotor activation appeared to correlate predominantly with stimulation at the more distal of the two fingertip jets (Fig. 2, top). Conversely, when the stimulation window approached the fingertip in a base-to-tip direction, the greater part of the sensorimotor activation appeared to correlate with stimulation at the base of the fingertip (Fig. 2,  bottom). It thus appeared that there was an enhanced population BOLD response limited to the part of the fingertip initially contacted in each cycle of stimulation.
The digit-related activation in both motor and somatosensory cortices was not only extensive and biased towards fingertip representation, but was also markedly symmetric with respect to the central fundus. Such symmetry was occasionally evident in individual pseudocoronal slice planes that happened to include both pre-and postcentral gyri (Fig. 3). This symmetry was not only the result of activation clumps that spanned the banks of the central sulcus; such symmetric activation was also evident for activation on the pre-and postcentral gyral surface some distance from the sulcus (e.g. Fig. 3, top). However, the convoluted path of the central sulcus-which can cross into such 2D slices multiple times-makes it difficult to perceive the topographical distribution of the functional activity around the sulcus. Nevertheless, these images do demonstrate the tuning of the functional activity to fingertip stimulation, as can be seen in comparing the tip-to-base (Fig. 3A) and base-to-tip experiments (Fig. 3B).
To visualize and quantify the symmetrical layout of the sensorimotor activation, we mapped the functional data to transformed, 2D views of the pre-and postcentral gyri (Fig. 4,  top, showing the same two experiments as in Fig. 3). These flattened representations of the pericentral cortex were found using the Isomap algorithm (see Methods). We binned the functional data of MI (area 4) and SI (areas 3b and 1) into square grid cells rostral or caudal to the fundus (Fig. 4, middle). Then, we computed the degree of similarity between the averaged phase values within each pair of active grid cells lying across the fundus from each other (Fig. 4, bottom). In the experiments shown, the mean symmetry score across all such pairs of grid cells was 0.38 in (A) and 0.37 in (B). Each of these values exceeded the 95 th percentile of the distribution of symmetry scores obtained from 1000 Monte Carlo simulations of each experiment (see Methods), which was 0.16 in both cases. r Figure 4. Activation tuning and symmetry in flattened representations of pericentral cortex. Shown are the same two experiments represented in Figure 3. In each plot the cortex is shown as a 2D surface, oriented in approximate rostrocaudal and lateromedial directions. The top plot shows functional data, representing phase lag according to the reference bar at top, superimposed over MI (rostral) and SI (caudal) voxels in gray, as per the legend used in Figures 2 and 3. (Area 3a, between these regions at the fundus of the central sulcus, is not shown.) The middle plot shows these same data binned into grid cells tiling the pericentral cortex. ''Active'' grid cells are indicated by the average phase lag across voxels within the bin. Grid cells representing pericentral cortex but lacking any active functional voxels are shown as gray (again, as per the Fig. 2 and 3  Indeed, the average phase symmetry score across all pairs of active grid cells was in many cases significant, in comparison to the simulated data sets. In particular, in the tip-to-base and base-to-tip digit 1 experiments, we found MI/SI symmetry scores in excess of threshold in 5/6 and 5/6 subjects, respectively. In digit 2 experiments, 5/6 (tip-to-base) and 4/6 (base-to-tip) subjects demonstrated significant MI/SI symmetry. In digit 4 experiments (one of them invalidated by motion artifact; see Methods), we found significant symmetry scores in only 2/6 and 3/5 subjects. Thus while the BOLD response following stimulation at least of digits 1 and 2 may have been globally tuned to fingertip phase lags, at a finer level most of these response maps could also be characterized as relatively symmetric mosaics of activation on the two sides of the central sulcus.

DISCUSSION
Our experimental methodology has previously allowed us to locate discrete somatotopic maps of the digits in areas 3b and 1 [9]. These maps were defined as regions of connected voxels displaying a strong correlation to the pattern of stimulation across the digit surface, and a reversed pattern of peak correlations when the stimulation direction was reversed. However, this analysis did not reflect the larger pattern of functional activation in the region of the central sulcus. While our delineation here of areas 1, 3b, 3a, and 4 was crude for lack of precise cytoarchitectonic divisions [18], the populationlevel activity suggests in any case that these areas had a relatively common response to digit stimulation. Notably, the activation in each area appeared to be weighted towards fingertip stimulation of the thumb and index finger (Fig. 2), in accordance with the high sensitivity of the distal finger pad and the relative importance of D1 and D2 in precision grips and other behaviors.
Unlike areas 3b and 1, areas 3a and 4 are conventionally ascribed roles in motor output [19]. Area 3a is thought to transmit kinesthetic rather than cutaneous afferents to both motor cortex [20] and area 1 [21]. But in addition to muscle spindle input to area 3a, in monkeys cutaneous responses here can emerge with training [22] and may physically parallel those in SI [20]. The presence of somatosensory inputs to area 4 has also been neglected (as has been the motor function of parietal areas, from which a significant proportion of motor corticospinal fibers are derived, in humans [23] as in monkeys [24]). Nevertheless, single-unit recordings from awake monkeys reveal cutaneous inputs to MI that appear to segregate in modality-specific maps, with neighboring representations being differentially responsive to cutaneous and joint receptor input [25]. In humans, despite fMRI evidence that a vibrotactile stimulus can elicit precentral gyrus activity [26,12], few researchers have tried to localize cutaneous functions to the frontal lobe.
We observed not only robust and fingertip-weighted activation on both sides of the central sulcus following stimulation of digits 1 and 2, but a symmetrical pattern to this activation. This is the first observation known to us of activation symmetry across functiondefined modalities. Within-modality mirror symmetry has recently been reported within tonotopic maps of primary auditory cortex [27] and object representations within occipito-temporal cortex [28].
This symmetry was evident not only despite the different functions of SI and MI, but also despite a less orderly somatotopic organization in precentral relative to postcentral areas. For instance, in monkeys the area 4 and 3a representations appear to be more fractured than the SI maps [20,22], perhaps consistent with the involvement of these areas in coordinated muscle recruitment. Within humans, a somatosensory homunculus has not been defined precentrally, although the motor output of MI is known to be somatotopically organized as shown by electrical stimulation in epileptic patients [1,29].
We suggest that both non-somatotopic organization in precentral cortex and the nature of the BOLD signal may underlie previous investigators' inability to resolve discrete foci of hemodynamic activity here upon digit stimulation [9,12]. In contrast to the spikingdefined somatotopy found with electrophysiological mapping, fMRI appears to reflect a more distributed, dendritic-level processing of neural inputs [6]. The symmetry of SI and MI responses we observed may reflect common sensory input to these areas. Although the symmetry may also reflect a diffusion of the hyperoxic response [30] in the vasculature on either side of the Rolandic artery, Young et al. [18] found that in the human left hemisphere, the resting regional blood flow within the area 3b classical ''hand'' area was correlated to blood flow within both anterior MI and area 3a.
In conclusion, while there may have been hotspots of BOLD response with some degree of somatotopic organization, the broader pattern of activity we observed upon cutaneous stimulation of the digits, in particular of the thumb and index fingers, was widely and symmetrically distributed on both sides of the central sulcus. In the case of thumb and index finger stimulation, the responses were also tuned to fingertip stimulation and were globally modulated by the direction of stimulation along the finger surface.