Four-dimensional maps of the humansomatosensory systemPietro Avanzinia,b,1, Rouhollah O. Abdollahia,c, Ivana Sartorid, Fausto Caruanae, Veronica Pellicciad, Giuseppe Casacelid,Roberto Maid, Giorgio Lo Russod, Giacomo Rizzolattia,b,e,1, and Guy A. Orbana

aDipartimento di Neuroscienze, University of Parma, 43125 Parma, Italy; bIstituto di Neuroscienze, Consiglio nazionale delle Ricerche, 43125 Parma, Italy;cCognitive Neuroscience Section, Institute of Neuroscience and Medicine, Research Center Jülich, 52428 Julich, Germany; dCentro per la chirurgiadell’Epilessia “Claudio Munari,” Ospedale Ca’Granda-Niguarda, 20162 Milan, Italy; and eBrain Center for Social and Motor Cognition, Istituto Italiano diTecnologia, 43125 Parma, Italy

Contributed by Giacomo Rizzolatti, February 5, 2016 (sent for review November 21, 2015; reviewed by Riitta Hari, Jon H. Kaas, Philippe Kahane,and Karl Zilles)

A fine-grained description of the spatiotemporal dynamics of humanbrain activity is a major goal of neuroscientific research. Limitations inspatial and temporal resolution of available noninvasive recordingand imaging techniques have hindered so far the acquisition ofprecise, comprehensive four-dimensional maps of human neuralactivity. The present study combines anatomical and functional datafrom intracerebral recordings of nearly 100 patients, to generatehighly resolved four-dimensional maps of human cortical processingof nonpainful somatosensory stimuli. These maps indicate that thehuman somatosensory system devoted to the hand encompasses awidespread network covering more than 10% of the cortical surfaceof both hemispheres. This network includes phasic components,centered on primary somatosensory cortex and neighboring motor,premotor, and inferior parietal regions, and tonic components,centered on opercular and insular areas, and involving human parietalrostroventral area and ventral medial-superior-temporal area. Thetechnique described opens new avenues for investigating the neuralbasis of all levels of cortical processing in humans.

intracranial recordings | stereo-EEG | median nerve | cerebral cortex |temporal dynamics

Adetailed description of the spatiotemporal dynamics of hu-man brain activity is a major goal of neuroscientific re-

search. However, it has been impossible so far to attain both highspatial and temporal resolution using the available noninvasiverecording and imaging techniques. Hence, a precise and com-prehensive four-dimensional cartography of human neural activityhas not yet been obtained. High spatial resolution, provided byneuroimaging techniques such as functional magnetic resonanceimaging (fMRI), is crucial for highlighting the topographical orga-nization of specific areas (e.g., somatotopy of sensorimotor areas)as well as identifying the nodes of brain networks endowed withspecific functional properties (1). It is not sufficient, however, toknow which nodes are active; information is also needed aboutthe local dynamics of the nodes, as well as the relative timing oftheir activity, to fully understand human brain functions (2, 3).Even if the temporal resolution of electroencephalography(EEG) and magnetoencephalography (MEG) allowed one toobserve the intra- and interareal dynamics, to date such re-cordings remain too poor in localization power (1–2 cm) (3, 4).Combining EEG and fMRI has been suggested as a solution,using EEG to determine the temporal dynamics within and be-tween the areas identified with fMRI (5). However, the disparatenature of the two signals recorded (hemodynamic for fMRI,electrical for EEG) creates discrepancies in the results that preventprecise matching of these methods (3).Invasive intracranial EEG offers a unique opportunity to observe

human brain activity with an unparalleled combination of spatialand temporal resolution. Depending on the electrodes used, twokinds of recordings can be made: (i) intraparenchymal recordings,also called stereo-EEG (sEEG) (6), obtained using stereotacticallyinserted needle-like electrodes with multiple recording leads; and

(ii) the electrocorticogram (ECoG), obtained from subdural elec-trode grids, covering regions of cortex. The latter technique suffersfrom three disadvantages: (i) the technique only samples fromcortical gyri, missing cortical regions buried in sulci (7); (ii) thetechnique is affected by volume conduction (8); and (iii) thetechnique does not record directly from gray matter, because piaand arachnoid mater lie in between the electrodes and the cortex,reducing the amplitude and making it more difficult to extract high-frequency activity from the recorded signal. Both approaches sufferfrom the so-called sparse-sampling problem (i.e., the limited anduneven coverage of the patients’ brain because the positioning of theelectrodes in any given patient is dictated by clinical criteria), leadingto difficulties in carrying out analyses at the population level.The primary aim of the present study is to show how ana-

tomical and functional data recorded from a large number ofpatients could be combined to generate highly resolved four-dimensional maps of human cortical processing. To demonstratethe feasibility of computing such maps, indicating consistentlyactive cortical nodes and the nodes’ time course at a millisecondscale, we leveraged the advantages of sEEG to investigate so-matosensory processing following electrical stimulation of themedian nerve in nearly 100 patients. fMRI studies have consis-tently reported activation of the primary somatosensory complex(SI), secondary somatosensory complex (SII), and insula in responseto transient nonpainful stimuli (9) but less so of supplementary

Significance

Here, we show how anatomical and functional data recordedfrom patients undergoing stereo-EEG can be combined togenerate highly resolved four-dimensional maps of humancortical processing. We used this technique, which providesspatial maps of the active cortical nodes at a millisecond scale,to depict the somatosensory processing following electricalstimulation of the median nerve in nearly 100 patients. The re-sults showed that human somatosensory system encompassesa widespread cortical network including a phasic component,centered on primary somatosensory cortex and neighboringmotor, premotor, and inferior parietal regions, as well as a toniccomponent, centered on the opercular and insular areas, lastingmore than 200 ms.

Author contributions: P.A., I.S., G.R., and G.A.O. designed research; I.S., V.P., G.C., R.M.,and G.L.R. performed research; P.A., R.O.A., and F.C. analyzed data; and P.A., G.R., andG.A.O. wrote the paper.

Reviewers: R.H., School of Science, Aalto University; J.H.K., Vanderbilt University; P.K.,INSERM U836 and Université Joseph Fourier; and K.Z., Research Center Julich.

The authors declare no conflict of interest.

Freely available online through the PNAS open access option.1To whom correspondence may be addressed. Email: pietro.avanzini@gmail.com orgiacomo.rizzolatti@unipr.it.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1601889113/-/DCSupplemental.

E1936–E1943 | PNAS | Published online March 14, 2016 www.pnas.org/cgi/doi/10.1073/pnas.1601889113

Dow

nloa

ded

by g

uest

on

Janu

ary

24, 2

020

motor area (SMA) (10), anterior cingulate cortex (11), superiorparietal lobule (SPL) and inferior parietal lobule (IPL) (12), andpremotor cortex (13). In addition to providing sensitive detectionof responsive regions, sEEG responses also reveal time courses,which could be relevant to debates concerning cortical dynamics,such as the parallel vs. serial operation of SII relative to SI (14, 15).Here, we show for the first time, to our knowledge, that compre-hensive four-dimensional maps of human cortical somatosensoryprocessing can be computed from high-frequency broadbandgamma activity, thus surpassing previous subdural and EEG/MEGmaps in spatiotemporal precision.

ResultsRecordings were obtained from 17,009 leads in 99 patients, ofwhich 11,983 (5,678 in the left hemisphere and 6,305 in the right)were localized in the cortical gray matter according to the ana-tomical reconstruction procedure (Fig. 1 and Figs. S1–S3). Ofnote, 48.8% of cortical leads (2,733 in the left hemisphere and3,112 in the right) were localized in sulcal areas [i.e., corticalregions characterized by a negative folding index (Materials andMethods)], many of which are not accessible to ECoG recordings.The sampling density maps computed for the two hemispheres(Fig. 1D and Figs. S2 and S3) demonstrate the extensive coverageof cortical sheet, with peaks bilaterally in the middle temporalgyrus and in the mesial temporal region (32 leads/cm2 in the left,48 leads/cm2 in the right hemisphere). Only the frontal and oc-cipital tips of the hemispheres as well as the cortical crowns werepoorly represented because of the obligatory orthogonal insertionof electrodes and to the anatomical and vascular constraints (pres-ence of frontal bone sinus and superior sagittal sinus, respectively).

Statistical analysis revealed that 1,146 of the leads exploring graymatter presented a significant broadband gamma power (50- to150-Hz) increase (see Fig. 1C for an example of single lead analysison gamma-band time course) in response to the median nervestimulation (489 in the left hemisphere, 657 in the right; Fig. S4).By localizing each lead, we computed the overall responsiveness(in % of responsive leads) maps for both hemispheres (Fig. 2), inprinciple comparable with neuroimaging data in that timing is notconsidered. High responsiveness values were found in SI. In par-ticular, the highest proportion of responsive leads was present bi-laterally in areas 3a (left 81%, right 94%) and 3b (left 95%, right87%). Lower values characterized areas 1 (left 60%, right 60%)and 2 (left 55%, right 64%). These results are in agreement withprevious MEG results (16). Opercular (OP) and insular regionswere active in both hemispheres, peaking at 50% in OP1 (17).Responsive regions included mainly OP1, OP2, and the long gyri ofinsular cortex (LgI), with both LgI active in the right hemisphereand the posterior LgI active in the left hemisphere. OP3, OP4, aswell as short gyri of the insula (SgI) were poorly responsive (Fig. 2and Fig. S4). OP1 and OP4 [corresponding to areas S2 (secondarysomatosensory area) and PV (parietal ventral area), as defined byDisbrow et al. (18)] are the main parts of classically defined SII.High proportions of responsive leads were observed in the

motor system, including the primary motor area (maximumvalues around 85%), large sectors of dorsal and ventral premotorcortex, and SMA. Anterior intraparietal sulcus was reliably re-sponsive in both hemispheres, approaching 60% in its mostrostral extent. Further posterior parietal cortex activations werebiased in favor of the left hemisphere in supramarginal gyrus andin a region located between OP4 and area 44 resembling parietalrostroventral (PR) area, described by Disbrow et al. (19) in the

Fig. 1. Anatomical and functional data analysis. (A) Axial MR slice of a patient, after the coregistration with the postimplantation CT, on which the entire Y′electrode implanted in the left hemisphere is visible. The reconstruction of the Y′ electrode is shown on the right side, with the leads numbered 1–18 from thetip, for a total distance of 50 mm. Leads are represented by sets of voxels, colored red if located in cortical gray matter and green otherwise. (B) Midthicknesssurface of the fs_LR brain template with all leads located in gray matter of the left hemisphere (from 49 patients) indicated as black dots. For four nodes,circles with 0.5 cm (thick line) and 1 cm (thin line) of radius represent the masks. (C, Upper) Average time frequency plot (100 trials) in response to mediannerve stimulation in one lead of one patient. (C, Lower) Time course of the average gamma power (50–150 Hz) for the same channel, reported as z scoresbased on the prestimulus interval. Black asterisks indicate time bins with gamma power significantly exceeding baseline (P < 0.001). (D) Sampling density ofthe left hemisphere computed from data in B. The color scale is expressed in the number of recording leads per cm2. Cytoarchitectonic areas (1–4, 6, OP1–4,and 44) and anatomically defined areas (long and short gyri of insula) are indicated by white lines, and functionally defined regions are indicated in blue[confidence ellipses of phAIP, DIPSA, and DIPSM from Jastorff et al. (78)] or green [MT cluster from Abdollahi et al. (72)].

Avanzini et al. PNAS | Published online March 14, 2016 | E1937

NEU

ROSC

IENCE

PNASPL

US

Dow

nloa

ded

by g

uest

on

Janu

ary

24, 2

020

monkey. This area may correspond to OP6 of Amunts and co-workers (20, 21). Interestingly, a weak but reliable responsiveness(around 20%) was found in right middle temporal (MT) cluster(22), extending dorsally into middle temporal gyrus (MTG).This activation partially overlaps with hOc5 (23).To rule out the strong responsiveness of motor and premotor

areas being due to the intensity of stimulation, just above motorthreshold, we computed overall responsiveness for a second setof recordings, obtained in 70 patients, in which the intensity was20% below the motor threshold. No difference in the propor-tions of active leads (overall responsiveness >10%) of sensory(BA2) and motor (BA4) areas were found between the supra-and subthreshold data in the left hemisphere (χ2 = 3.97; notsignificant). A significant difference (after correction for twocomparisons) was found in the right hemisphere (χ2 = 14.75; P <0.01), but the proportion of active leads decreased more in BA2than in BA4 (Table S1) with subthreshold stimulation.Clustering (k means) the time course of gamma band activity

(Table S2) over all responsive leads indicated five to be the

optimal number of clusters (average silhouette value: 0.448). Thecentroids obtained for each cluster (Fig. 3B) depicted three syn-chronous patterns (strong-, middle-, and weak-phasic), charac-terized by a steep rise after stimulus delivery, a peak in the thirdtime bin (20–30 ms) and a return to baseline level within 50 msfrom stimulation. A prolonged cluster (green line) exhibited asimilar time course in the rising phase, but persisted two time bins(20 ms) longer than phasic clusters. A completely different patternwas shown by a tonic cluster (blue line) lasting more than 200 ms,whose peak was consistently lower and later than all other clusters.Quality of the clustering was indicated by the small percentage ofnegative silhouette values (less than 5% of leads, most of thembelonging to middle-phasic and prolonged clusters, whose peakamplitudes were very similar). These leads were excluded fromsubsequent cluster-related analyses.The relative responsiveness (percentage of responsive leads

belonging to a cluster) maps showed a remarkable topographicaldistribution across the two hemispheres. As shown in Fig. 3 A–C,the strong phasic cluster (see red palette) covered the middle

Fig. 2. Overall responsiveness maps. Overall responsiveness (responsive leads as a percentage of total explored leads per disk) maps for the left (A) and right(B) hemispheres. Only surface nodes with values exceeding 10% were shown. The same conventions as in Fig. 1 are used.

Fig. 3. Clustering of time courses. (Center) Time courses (centroid ± SD) of the five clusters, as indicated (Inset). (Left and Right) Relative responsiveness(leads belonging to one cluster as a percentage of total number of responsive leads per disk) maps of left (L) and right (R) hemispheres (middle portion) forstrong-phasic (yellow red), delayed-phasic (green), and tonic (blue) clusters. Only nodes with values exceeding 20% are shown. The same conventions as inFig. 1 are used.

E1938 | www.pnas.org/cgi/doi/10.1073/pnas.1601889113 Avanzini et al.

Dow

nloa

ded

by g

uest

on

Janu

ary

24, 2

020

strip of primary somatosensory areas (peaks of 81% in left 3b and70% in right 3a), primary motor area and anterior part of inferiorparietal cortex, corresponding to the caudal sectors of PF, PFt,and PFm (24, 25). The prolonged cluster (green palette) was re-stricted to a small sector of dorsal premotor area bilaterally (bothreaching 33%), area 4 in the right hemisphere and SMA in the lefthemisphere. The tonic cluster (see blue palette) was highly specificfor secondary somatosensory areas (OP1 and OP2) and long gyriof insular cortex (peaks around 80% and 90%). Also, activity inhuman PR, ventral premotor cortex (bilaterally), and MT cluster(right hemisphere) exhibited a tonic time course. Middle- andweak-phasic clusters (Fig. S5) showed less focal distributions, in-stead spreading centrifugally from the primary sensory areas withthe middle-phasic cluster confined to closer areas and the weak-phasic cluster to more distant ones. Of note, these two clusterslargely avoided the regions belonging to the three focal clusters.The region of interest (ROI) analysis confirmed that all four

primary somatosensory areas are characterized by a strong-phasictime course, with areas 3a, 3b, and 1 similar in terms of peakamplitude (Fig. 4A). A repeated-measures ANOVA returnedsignificant main effects for both time and ROIs, and a significantROI-by-time interaction [F(3,177) = 2.787; P < 0.0001]. However,post hoc comparisons did not reveal any consistent differenceacross ROIs. The same analysis conducted on areas of the motorsystem (Fig. 4B) showed significant main effects for both time andROIs, and a significant ROI-by-time interaction [F(4,236) = 6.545;P < 0.0001]. Interestingly, post hoc analysis identified three timeintervals differentiating the motor ROIs: (i) around the peak re-sponse (synchronous for all ROIs), the medial premotor areashows significantly lower power values than area 4, putative hu-man anterior intraparietal area (phAIP), and dorsal premotorarea (PMd); (ii) in the 50- to 70-ms interval, phAIP shows nosustained power, in contrast to all other ROIs, with PMd pre-senting the highest values; (iii) between 100 and 200 ms, power inthe ventral premotor area (PMv) and the medial premotor areasignificantly exceeded that of PMd, probably reflecting the pres-ence of tonic cluster activity (Fig. 3). A repeated-measures ANOVArun on secondary somatosensory areas returned a significant ROI-by-time interaction [F(2,118) = 4.367; P < 0.0001]. Post hocanalysis showed a significant difference between OP1 and long gyriof insula, and between OP1 and other OPs combined in the 20- to30-ms period, with OP1 presenting a stronger phasic responserelative to other ROIs. In a late interval (100–200 ms), however,gamma power recorded from LgI exceeded that in opercular areas.The motor temporal pattern was tested also for the sub-

threshold dataset, with the specific aim of validating that thesustained behavior of PMd between 50 and 70 ms was unrelatedto thumb twitch induced by the median nerve stimulation. TheROI-by-time interaction again reached significance [F(4,236) =1.294; P < 0.005] with subthreshold stimulation. Post hoc anal-ysis, comparing mean gamma power across ROIs in the 50- to70-ms time window confirmed the prolonged activation in motorand premotor areas relative to phAIP (Fig. S6).Clustering and ROI analyses have provided valuable insights

into the dynamics of hand somatosensory processing. Bothanalyses have limitations because the first lumps large numbersof leads together, and the second blends leads selected on apriori anatomical criteria. To provide data-driven evidence, weapplied map computation directly to gamma power amplitude ateach time bin. In this way, we constructed a temporal sequenceof maps, one every 10 ms, depicting the distribution of averageamplitude (considering only responsive leads, see Fig. S7 A–C)over cortex. To complement this information, we time-resolved theoverall responsiveness maps (Fig. 2) to achieve time dependentresponsiveness maps (Fig. S7 D–F). Both datasets contain distinctand independent information, related to variations in time ofstrength and reliability of response. Hence, we combined them bymeans of a pointwise multiplication, thus obtaining weighted am-

plitude maps over time. Because this operation is performed foreach time-bin, a movie that links together the entire sequence vi-sualizes the 4-dimensional dynamics of the cortical activity in leftand right hemispheres (Movies S1–S4). By inspecting four framesof such movie (Fig. 5) one can easily see: (i) the prevalence of thestrong phasic response in area 3b, 3a and 1 at 20–30 ms (Fig. 5A);(ii) the marked involvement of posterior part of dorsal premotorcortex visible around 40–50 ms after stimulation (Fig. 5B); (iii) thedominance of premotor cortex and OP1 at 60–70 ms after stimu-lation, when primary sensory areas are already on the decrease(Fig. 5C); and, finally, (iv) the restriction of late activity (more than100 ms following the onset of the electrical stimulation) to OP1,OP2 and LgI.

Fig. 4. ROI analysis. Average (±SE) time course of leads in somatosensoryareas (A), motor regions (B), and opercular/insular areas (C). Areas are listed(Insets). Black marks under the curves indicate significant post hoc compar-isons (P < 0.002). Number of responsive leads for each graph: 34 in area 1, 54in area 2, 21 in area 3a, 50 in area 3b (A); 68 in area 4, 157 in PMd, 45 inmedial premotor areas, 42 in PMv, and 28 in phAIP (B); and 133 in OP1, 54 inother OPs, and 38 in LgI (C).

Avanzini et al. PNAS | Published online March 14, 2016 | E1939

NEU

ROSC

IENCE

PNASPL

US

Dow

nloa

ded

by g

uest

on

Janu

ary

24, 2

020

DiscussionIn the present study, by localizing intracerebral electrodes inindividual hemispheres, assessing the gamma broadband activityof all recording leads, and mapping the results to a commontemplate, we obtained a comprehensive four-dimensional car-tography of the cortical processing of median nerve stimulation.Beyond its neuroscientific relevance, this knowledge may allow amore careful interpretation of the sEEG activity and a betterdesign of the sEEG implantations based on the ictal clinicalsemiology.This result was achieved thanks to the group analysis of in-

tracranial recordings from a large number of patients. Similarattempts to solve the sparse-sampling problem have been re-cently made using ECoG (26–31). However, although useful forinvestigating activity from areas near the crowns of gyri, suchrecordings are blind to neuronal activity from deep regions (7).It is worth noting that half of the leads recorded in the presentstudy by sEEG were localized in sulci, including major sulci forsomatosensory processing like the sylvian and the rolandic fis-sures. Previous sEEG studies (32–35) reported maps derivedfrom population analyses, but the limited number of patientsprevented them from obtaining the full picture of corticalprocessing.High-frequency broadband gamma activity (50–150 Hz) was

extracted from sEEG as an index of cortical activity, beingspatially and functionally more specific than other frequencybands and presumably reflecting population spiking activity(36, 37), and thus allowing networks to be investigated at mil-lisecond time-scales (33). A comparison between gamma timecourse and ERP responses from the same leads is presented inFigs. S8 and S9. To obtain a precise spatial localization forcontinuous maps, we minimized the circular kernel size by usingthe geodesic distance (38) and logistic weighting function. Ofnote, masks centered at 1.5-cm geodesic distance were com-pletely independent, having no shared nodes.

Responsiveness Maps Document the Wide Extent of the CorticalSomatosensory Network. The overall responsiveness (Fig. 2) andweighted amplitude maps (Fig. 5) reveal the cortical networkinvolved in processing somatosensory stimuli and document itstime course. This network encompasses several areas the in-volvement of which remained controversial so far. Indeed, besideprimary somatosensory cortices and SII/insula region (39), werevealed the activation of motor and premotor cortex (13), SMA(10, 40), superior and inferior parietal lobules (12, 41, 42). Thus,despite the very stringent response criteria we used, computedmaps indicate a widespread network (over 10% of the surface inboth hemispheres) devoted to processing median nerve input. Thisresult is even more notable given that the present report is limitedto stimulation of only one of the three nerves innervating the hand.In addition, two additional regions were identified as responding

to median nerve stimulation. The first is located in both hemi-spheres in the same sector of parieto-ventral region that has beencalled human PR (19). This area appears to be the counterpart ofmonkey PR, which is connected to IPL, area 4, and premotorcortex (19, 43). Disbrow et al. (18) reported that human PRresponded to passive tactile stimulation, although in only 60% oftested participants. Investigating the temporal relationship be-tween areas S2 and PR using MEG, Hinkley et al. (44) concludedthat both S2 and PR contribute to the somatomotor integrationneeded for manual exploration and object discrimination, with S2activity preceding that of PR. The tonic time course observed inhuman PR in the present study supports the involvement of thisregion in these discriminative functions.The second site was revealed by leads exploring the ventral

medial–superior–temporal area (pMSTv), part of the MT cluster(22), and neighboring posterior MTG in the right hemisphere(Fig. 2B). This region is classically considered involved in high-order visual processing, whereas its contribution to somatosensoryprocessing is controversial. van Kemenade et al. (45) reportedthat different neuronal populations in human middle temporal

Fig. 5. Weighted amplitude maps over time. Four weighted amplitude maps of the left hemisphere at time bins indicated. The color range is dynamicallyadjusted from frame to frame. Note how the activity peaks in primary sensory areas in first frame, spreads to dorsal premotor area (second frame), remainsactive longer (third frame), whereas residual activity after 100 ms is recorded mostly from parietal operculum (OP1) and long gyri of insula (last frame). Thesame conventions as in Fig. 1 are used.

E1940 | www.pnas.org/cgi/doi/10.1073/pnas.1601889113 Avanzini et al.

Dow

nloa

ded

by g

uest

on

Janu

ary

24, 2

020

(hMT)+/V5 process tactile and visual motion information. Incontrast, Jiang et al. (46) reported only a weak (2% of voxels) re-sponse to tactile stimulation. Given the high sensitivity of in-tracranial recordings, we were able to show that indeed somatosensoryinformation reaches pMSTv, suggesting a possible integration be-tween tactile and visual information (47). If pMSTv functions inthe same capacity as its monkey counterpart, this integration mayserve to control tracking of moving objects not only with eyes, asdescribed by Newsome et al. (48), but also with arms (49).

Four-Dimensional Maps Provide Previously Unidentified InformationAbout Interareal Dynamics. There is an ongoing debate of whetherSI and SII operate serially or in parallel (14, 15). Both views arecompatible with anatomical findings demonstrating that SI isconnected to SII via cortico–cortical connections (50, 51) andthat various thalamic nuclei project in parallel to SI and SII(52). Neuroimaging studies addressed this issue using causalmodeling approach, but with mixed results: Kalberlah et al.(53) and Khoshnejad et al. (54) supported serial processing,whereas Chung et al. (55) and Liang et al. (56) the parallelmodel. Even intracranial studies yielded conflicting results.Barba et al. (57) reported a negative ERP peaking at 30 ms inparietal operculum following the first negative component oc-curring at 20 ms waveform, thus favoring the serial interpretation.On the contrary, Karhu and Tesche (58) observed synchronousneuronal population activity in contralateral SII area 20–30 ms afterstimulation, coinciding in time with the first responses generatedin SI, thus supporting the parallel hypothesis.In the present study, the phasic component in OP1 coincided in

time with the phasic activity in SI, thus favoring the view of aparallel propagation of sensory information to the two areas.However, because of the 10-ms time bin of our analysis, this con-clusion should be taken with caution. Interestingly, our findingsdistinguished different patterns of activation in peri-sylvian areas,with OP1 presenting the strongest phasic component followed by atonic pattern, and OP2, OP3, and OP4 exhibiting mainly tonicactivity (Fig. 4C). The relative timing of phasic and tonic activitiesin operculoinsular areas is consistent with previous MEG data(59, 60). Hence, OP1 appears to play the role of a hub dispatchinginformation to neighboring opercular cortices (61). To settle thispoint, further studies are needed using for example functionalconnectivity techniques like cortico–cortical evoked potentials (62).Finally, our data showed that insular cortex was consistently

activated by median nerve stimulation and presented a markedlytonic time course (Fig. 3). This activation was limited to theposterior part of insula, anatomically corresponding to the in-sular long gyri. This finding is in agreement with the metaanalysisby Kurth et al. (63), which demonstrated that this part of theinsula is a specific sector related to somatosensory processing.Our finding of a tonic activity of this sector may be related to theinvolvement of posterior insula in the integration of somato-sensory inputs with other sensory modalities (64).

Materials and MethodsParticipants. Stereo-EEG data were collected from 99 patients (52 male, 47female) suffering from drug-resistant focal epilepsy. Only patients present-ing with no anatomical alterations (n = 78) or with small abnormalitiesoutside of the sensorimotor areas (n = 21), as evident on MR, were included.Fifteen of the patients with positive MR showed alterations of the temporallobe [4 hippocampal sclerosis, 4 minimal periventricular nodular hetero-topia, 2 cavernomas of temporal pole, 2 focal cortical dysplasia (FCDII), and 3post ischemic injuries]. Four patients presented alteration of the posteriorparietal lobe (no overlap with active ROIs reported in the study), of whom,two patients had FCDII, one had cavernoma, and one had post-ischemicinjury. Finally, one patient presented a FCDII of the occipital lobe, and onepatient had an anterior periventricular nodular heterotopia.

These patients were stereotactically implanted with intracerebral elec-trodes as part of their presurgical evaluation, at the Centro per la chirurgiadell’Epilessia “Claudio Munari.” Implantation sites were selected solely on

clinical grounds, using seizure semiology, scalp-EEG, and neuroimaging asguide. Patients were fully informed regarding the electrode implantationand sEEG recordings. The present study received the approval of the EthicsCommittee of Ospedale Ca’Granda-Niguarda (ID 939-2.12.2013) and in-formed consent was obtained. Intracerebral recordings were performedaccording to sEEG methodology to define the cerebral structures involved inthe onset and propagation of seizure activity (65, 66). No seizure occurred,no alterations in the sleep/wake cycle were observed, and no additionalpharmacological treatments were applied during the 24 h before the ex-perimental recording. Neurological examination was unremarkable in allcases; in particular, no motor or sensory deficit was found in any patient.

Electrode Implantation. Most implantations were unilateral, because clinicalevidence generally indicates the hemisphere generating the seizures. Only8 of the 99 patients were implanted bilaterally, resulting in a total of107 implanted hemispheres. A number of depth electrodes (range: 9–19;average: 13) were implanted in different regions of the hemisphere usingstereotactic coordinates. Each cylindrical electrode had a diameter of 0.8 mmand consisted of eight to eighteen 2-mm-long contacts (leads), spaced1.5 mm apart (DIXI Medical, Besancon, France).

Immediately after the implantation, cone-beam computed tomography(CBCT) was obtained with the O-arm scanner (Medtronic) and registered topreimplantation 3D T1-weighted MR images. Subsequently, multimodalscenes were built with the 3D Slicer software package (67), and the exactposition of each lead was determined, at the single patient level, looking atmultiplanar reconstructions (68). Following clinical conventions, all leads areidentified by a letter corresponding to the electrode shaft, followed by anumber starting from the tip of the electrode.

Median Nerve Stimulation. The day after the implantation, patients wereadmitted to the neurology ward, to undergo clinical and neuropsychologicaltests to functionally characterize the recording leads. To map leads involvedin processing hand somatosensory information, a median nerve stimulationtest was administered to the patients lying in bed with the eyes closed. Themedian nerve opposite to the recorded hemisphere was stimulated at thewrist, using 100 constant-current pulses (0.2-ms duration) at 1 Hz. The in-tensity and exact site of stimulation were varied until an observable thumbtwitch was obtained. The motor threshold in our sample ranged from 3.2 to5.8 mA. The stimulation intensity was set at 10% above the motor threshold.As a control, most of the patients (n = 70) were also tested with a stimulationintensity 20% below the motor threshold.

Anatomical Reconstruction of Electrodes. The aim of this computational stagewas to precisely locate the recording leads in individual cortical surfaces usingthe multimodal explorations performed in each patient, and import theselocations into a common template. Each patient underwent a structural MRI(Achieva; Philips Medical Systems) before electrodes implantation. The T1images were segmented using FreeSurfer software (69) to identify the pialand the white matter surfaces for the native mesh of each patient. Thequality of segmentation was verified by visual inspection of the resultingsurfaces. The midthickness surface, i.e., the average surface lying betweenthe pial and white matter surfaces, was extracted. Moreover, starting fromthe pial surface, the sulci pattern was extracted through a procedure thatevaluates the normalized geometric depth of each point [folding index 70)].The more negative this parameter, the deeper the point is buried in a sulcus.After the implantation, each patient underwent a volumetric brain CT, tolocate precisely the recording leads using the artifacts generated on CTimages (Fig. 1B). The MRI and CT datasets were coregistered [FMRIB’s linearimage registration tool (FLIRT), 6 df, mutual information (71)] to get theanatomical brain data and the implanted electrodes into the same co-ordinate space. Thresholding of the CT signal intensity was used to segre-gate the recording leads (Fig. 1B) and to reconstruct their position in the CTvolume. During this process, the knowledge of the number of implantedelectrodes, the number of leads on each electrode, and the electrodes’ entryand target points were used as constraints to ensure a reliable reconstruction.

Subtracting white-matter from pial surfaces yielded a ribbon surfacecorresponding to the gray matter thickness. Because the spatial resolution ofCT images exceeded that of theMRI imaging and the CT images were reducedto binary information (1: presence of a recording lead; 0: absence of it), theribbon was oversampled at a resolution of 0.4 mm in all directions. The in-tersection of CT images with the ribbon surface detected the contact pointslocated in the cortical gray matter (Fig. 1B, red dots). Thus, at this stage, alead located in the gray matter is represented by the number of 0.064-mm3

voxels it shares with the gray matter. For those leads in the gray matter, thecentroid of the artifact was computed and projected onto the nearest nodes

Avanzini et al. PNAS | Published online March 14, 2016 | E1941

NEU

ROSC

IENCE

PNASPL

US

Dow

nloa

ded

by g

uest

on

Janu

ary

24, 2

020

of the midthickness native surface, thus obtaining a single surface noderepresenting the projection of a lead located in gray matter.

Finally, the individual midthickness surface was resampled to match thenumber of nodes (163.842) of the template (Fs-LR-average) and coregisteredwith this template using Freesurfer_to_fs_LR pipeline (brainvis.wustl.edu/wiki/index.php/Caret:Download). In this step, single-lead representations were en-hanced by arbitrarily adding six nodes (on the hexagonal mesh) surroundingthe original single nodes on the native surfaces, to ensure that all leads in theindividual surfaces were maintained on the template surface. To quantify thegoodness of fit between the native meshes and the template, we computedtwo local indices [i.e., deformation and distortion, evaluated for each node ofthe mesh (72)], indexing the shrinkage and translation of each node during thetransformation. These indices were used to calculate the surface area of nodesin a group taking into account individual cortical surface area.

In conclusion, the recording leads located in the gray matter are repre-sented in three different formats in the reconstruction procedure: as thenumber of 0.064-mm3 native voxels and as a fixed number (7) of nodes inthe native and in the template meshes. Whereas the first format providesthe precise location in the cortical depth, the latter two provide precise in-formation about location in the cortical surface (at midthickness).

SEEG Data Recording and Processing. For each implanted patient, the initialrecording procedure included the selection of an intracranial reference,which was chosen by clinicians using both anatomical and functional criteria.The reference was computed as the average of two adjacent leads bothexploring white matter. These leads were selected time-by-time because theydid not present any response to standard clinical stimulations, includingsomatosensory (median, tibial, and trigeminal nerves), visual (flash), andacoustical (click) stimulations. Nor did the leads’ electrical stimulation evokeany sensory and/or motor behavior. The sEEG trace was recorded with aNeurofax EEG-1100 (Nihon Kohden System) at 1-kHz sampling rate.

Clinicians visually inspected recordings to verify for ictal epileptic dis-charges (IEDs) during the stimulation protocol. In 87 patients, the ROIs weredevoid of any IED. In the remaining 12 patients, sparse IEDs were recordedduring the stimulation protocol. In these cases, however, false-positive re-sponses are very unlikely because none of the patients included showed reflexIEDs, and thus IEDs were not synchronized to the stimulus. This absence wasconfirmed by visual inspection of the quality of the data averaged.

The recordings from all leads in the gray matter were filtered (band-pass:0.015–500 Hz; notch: 50 Hz) to avoid aliasing effects and decomposed intotime–frequency plots using complex Morlet’s wavelet decomposition. Powerin the gamma (50- to 150-Hz) frequency band was extracted in a windowextending from 100 ms before to 500 ms after the electrical stimulation, andsubdivided into 60 nonoverlapping 10-ms bins. Following previous in-tracranial studies (35, 73, 74), gamma power was estimated for 10 adjacentnonoverlapping 10-Hz frequency bands.

To identify the leads responsive to median nerve stimulation, the gammaband power in each poststimulus bin was compared with baseline using a ttest. Significance was Bonferroni corrected for 50 comparisons (P = 0.001),and to decrease the false-positive ratio, only leads with significant gammaincreases in at least three time bins were designated as responsive. To nor-malize data across patients and leads, power in poststimulus bins wastransformed into z scores relative to the prestimulus interval.

A k mean clustering was applied to the gamma band time course of theresponsive leads to group them according to the time course of the response,independently of their location. Nineteen clustering procedures were performedimposing increasing numbers of clusters (3–20) and computing silhouette values(75) to evaluate clustering validity. The optimal clustering was defined by themaximal average silhouette value, and leads with negative silhouette valueswere not considered in the evaluation of the optimal clustering using a re-peated-measurement ANOVA with time and cluster as factors. t tests (P < 0.001to account for 50 time bins) were used to evaluate post hoc comparisons.

Finally, to compare cortical areas, leads were grouped by their location in aROI-based analysis to evaluate differences among primary sensory areas (3a,3b, 1, and 2), motor areas (4, medial premotor area, PMd, PMv, and phAIP),

and the secondary somatosensory areas (OP1, OP2, OP3, OP4, and LgI). TheseROIs were defined either by cytoarchitectonical or functional criteria (seeContinuous Maps). For each of these groups, we performed a two-way re-peated-measure ANOVA with time and area as factors. Significant interac-tions were explored in post hoc analysis with planned-comparisons t tests(P < 0.001 to account for the 50 time bins).

Continuous Maps. To provide a continuous view of the topographic pattern ofactive leads, we built a circular mask based on the geodesic distance betweentwo cortical points (76) (i.e., the minimum pathway within the gray matterconnecting the source and the target nodes). For each cortical node, we definedthe nodes within a 1-cm geodesic distance from the original node andweighted the contribution of each node by a sigmoid function. Node weightwas defined as a logistic function with unitary amplitude, a steepness of 2 and amidpoint at 7.5 mm. As a result, each node of the cortical mesh was associatedwith a collection of surrounding nodes (averaging 806 and 813 nodes for theright and left hemispheres), with all nodes within 5 mm of the origin maximallyweighted, whereas those between 5 and 10 mm were gradually reduced inweight, to avoid edge effects.

By this approach, we computed five different functional variables:

i) Cortical sampling density [i.e., number of explored leads per cm2, using thefixed number (7) of nodes per lead and the average surface of a disk].

ii) Overall responsiveness [i.e., the number of responsive nodes (Nr) as apercent of the number of explored nodes (Ne) within a disk]. Given thefixed number of nodes representing a lead (7), this variable is equivalentto the number of responsive leads as a percent of the number of ex-plored leads within a disk. This time-independent variable provides anoverall picture of the cortical responsiveness, ranging from 0% to 100%,directly comparable to neuroimaging studies. Data were thresholded at10% so as to exclude the contribution of sparsely responsive regions.

iii) Relative responsiveness [i.e., the number of nodes (leads) exhibiting aspecific temporal response pattern in percent of the number of respon-sive nodes (leads) within a disk]. This variable indexes the degree towhich an area responds with a specific temporal pattern, informationto which neuroimaging studies are completely blind. Results werethresholded at a 1/N value, where N represents the number of clusters,to show only areas in which proportions of single clusters exceed chance.

iv) The response amplitude for each time bin (i.e., the average z scoreacross all responsive nodes in the disk). This variable indexes thestrength of a response, regardless of how reliably it occurs within a disk.No minimum threshold was set for the z score, but results were maskedfor overall responsiveness exceeding 10%.

v) The weighted amplitude for each time bin, i.e., a combination of re-sponse amplitude and overall responsiveness, depicting the variationover the cortical surface of the response amplitude weighted bythe responsiveness:

NrNe

�overall

responsiveness

�*ΣANr

�responseamplitude

�=ΣANe

�weightedamplitude

�,

where Ne is the number of explored nodes (leads) and Nr is the number ofresponsive nodes (leads) in a disk. Computation was restricted to regionswhose overall responsiveness was higher than 10%.

These maps were plotted using CARET software (77) (www.nitrc.org/projects/caret) and directly compared with retinotopic regions in the occipital cortex(72) to confidence ellipses of rostral intraparietal sulcus (78), the most rostral ofwhich corresponds to phAIP, and to cytoarchitectonic sensorimotor regions(79), opercular areas (17), and area 44 (80). The subdivision between dorsal andventral premotor areas was made according to Tomassini et al. (81).

ACKNOWLEDGMENTS. The authors thank Dr. S. Raiguel for revising theEnglish version of the manuscript. This study was supported by EuropeanResearch Council (ERC) “Cogsystem” project no. 250013 (to G.R.) and ERC“Parietalaction” project no. 323606 (to G.A.O.).

1. Van Essen DC (2013) Cartography and connectomes. Neuron 80(3):775–790.2. Kopell NJ, Gritton HJ, Whittington MA, Kramer MA (2014) Beyond the connectome:

The dynome. Neuron 83(6):1319–1328.3. Hari R, Parkkonen L (2015) The brain timewise: How timing shapes and supports brain

function. Philos Trans R Soc Lond B Biol Sci 370(1668):pii: 20140170.4. Burle B, et al. (2015) Spatial and temporal resolutions of EEG: Is it really black and

white? A scalp current density view. Int J Psychophysiol 97(3):210–220.5. Debener S, Ullsperger M, Siegel M, Engel AK (2006) Single-trial EEG-fMRI reveals the

dynamics of cognitive function. Trends Cogn Sci 10(12):558–563.

6. Cardinale F, Lo Russo G (2013) Stereo-electroencephalography safety and effective-

ness: Some more reasons in favor of epilepsy surgery. Epilepsia 54(8):1505–1506.7. Noy N, et al. (2015) Intracranial recordings reveal transient response dynamics during

information maintenance in human cerebral cortex. Hum Brain Mapp 36(10):3988–4003.8. Whitmer D, Worrell G, Stead M, Lee IK, Makeig S (2010) Utility of independent

component analysis for interpretation of intracranial EEG. Front Hum Neurosci 4:184.9. Ferretti A, et al. (2007) Cortical brain responses during passive nonpainful median

nerve stimulation at low frequencies (0.5-4 Hz): An fMRI study. Hum Brain Mapp

28(7):645–653.

E1942 | www.pnas.org/cgi/doi/10.1073/pnas.1601889113 Avanzini et al.

Dow

nloa

ded

by g

uest

on

Janu

ary

24, 2

020

10. Manganotti P, et al. (2009) Steady-state activation in somatosensory cortex after changesin stimulus rate during median nerve stimulation. Magn Reson Imaging 27(9):1175–1186.

11. Arienzo D, et al. (2006) Somatotopy of anterior cingulate cortex (ACC) and supple-mentary motor area (SMA) for electric stimulation of the median and tibial nerves: AnfMRI study. Neuroimage 33(2):700–705.

12. Ruben J, et al. (2001) Somatotopic organization of human secondary somatosensorycortex. Cereb Cortex 11(5):463–473.

13. Boakye M, Huckins SC, Szeverenyi NM, Taskey BI, Hodge CJ, Jr (2000) Functionalmagnetic resonance imaging of somatosensory cortex activity produced by electricalstimulation of the median nerve or tactile stimulation of the index finger.J Neurosurg 93(5):774–783.

14. Rowe MJ, Turman AB, Murray GM, Zhang HQ (1996) Parallel organization of so-matosensory cortical areas I and II for tactile processing. Clin Exp Pharmacol Physiol23(10-11):931–938.

15. Forss N, Hietanen M, Salonen O, Hari R (1999) Modified activation of somatosensorycortical network in patients with right-hemisphere stroke. Brain 122(Pt 10):1889–1899.

16. Kaukoranta E, Hämäläinen M, Sarvas J, Hari R (1986) Mixed and sensory nerve stim-ulations activate different cytoarchitectonic areas in the human primary somato-sensory cortex SI. Neuromagnetic recordings and statistical considerations. Exp BrainRes 63(1):60–66.

17. Eickhoff SB, Grefkes C, Zilles K, Fink GR (2007) The somatotopic organization of cy-toarchitectonic areas on the human parietal operculum. Cereb Cortex 17(8):1800–1811.

18. Disbrow E, Roberts T, Krubitzer L (2000) Somatotopic organization of cortical fields inthe lateral sulcus of Homo sapiens: Evidence for SII and PV. J Comp Neurol 418(1):1–21.

19. Disbrow E, Litinas E, Recanzone GH, Padberg J, Krubitzer L (2003) Cortical connectionsof the second somatosensory area and the parietal ventral area in macaque monkeys.J Comp Neurol 462(4):382–399.

20. Amunts K, et al. (2010) Broca’s region: Novel organizational principles and multiplereceptor mapping. PLoS Biol 8(9):e1000489.

21. Amunts K, Zilles K (2012) Architecture and organizational principles of Broca’s region.Trends Cogn Sci 16(8):418–426.

22. Kolster H, Peeters R, Orban GA (2010) The retinotopic organization of the humanmiddle temporal area MT/V5 and its cortical neighbors. J Neurosci 30(29):9801–9820.

23. Malikovic A, et al. (2007) Cytoarchitectonic analysis of the human extrastriate cortexin the region of V5/MT+: A probabilistic, stereotaxic map of area hOc5. Cereb Cortex17(3):562–574.

24. Caspers S, et al. (2006) The human inferior parietal cortex: Cytoarchitectonic parcel-lation and interindividual variability. Neuroimage 33(2):430–448.

25. Caspers S, et al. (2008) The human inferior parietal lobule in stereotaxic space. BrainStruct Funct 212(6):481–495.

26. Miller KJ, et al. (2007) Real-time functional brain mapping using electrocorticography.Neuroimage 37(2):504–507.

27. Burke JF, et al. (2013) Synchronous and asynchronous theta and gamma activityduring episodic memory formation. J Neurosci 33(1):292–304.

28. Conner CR, Chen G, Pieters TA, Tandon N (2014) Category specific spatial dissociationsof parallel processes underlying visual naming. Cereb Cortex 24(10):2741–2750.

29. Davidesco I, et al. (2013) Spatial and object-based attention modulates broadband high-frequency responses across the human visual cortical hierarchy. J Neurosci 33(3):1228–1240.

30. Dykstra AR, et al. (2012) Individualized localization and cortical surface-based regis-tration of intracranial electrodes. Neuroimage 59(4):3563–3570.

31. Kadipasaoglu CM, et al. (2015) Development of grouped icEEG for the study ofcognitive processing. Front Psychol 6:1008.

32. Juphard A, et al. (2011) Direct evidence for two different neural mechanisms for readingfamiliar and unfamiliar words: An intra-cerebral EEG study. Front Hum Neurosci 5:101.

33. Lachaux J-P, Axmacher N, Mormann F, Halgren E, Crone NE (2012) High-frequencyneural activity and human cognition: Past, present and possible future of intracranialEEG research. Prog Neurobiol 98(3):279–301.

34. Ossandón T, et al. (2012) Efficient “pop-out” visual search elicits sustained broadbandγ activity in the dorsal attention network. J Neurosci 32(10):3414–3421.

35. Vidal JR, et al. (2010) Category-specific visual responses: An intracranial study com-paring gamma, beta, alpha, and ERP response selectivity. Front Hum Neurosci 4:195.

36. Popivanov ID, Jastorff J, Vanduffel W, Vogels R (2014) Heterogeneous single-unitselectivity in an fMRI-defined body-selective patch. J Neurosci 34(1):95–111.

37. Ray S, Maunsell JHR (2011) Different origins of gamma rhythm and high-gammaactivity in macaque visual cortex. PLoS Biol 9(4):e1000610.

38. Kadipasaoglu CM, et al. (2014) Surface-based mixed effects multilevel analysis ofgrouped human electrocorticography. Neuroimage 101:215–224.

39. Kaas JH (2004) Somatosensory System. The Human Nervous System, eds Paxinos G,Mai JK (Elsevier Academic, New York), 2nd Ed, pp 1059–1092.

40. Korvenoja A, et al. (1999) Activation of multiple cortical areas in response to so-matosensory stimulation: Combined magnetoencephalographic and functionalmagnetic resonance imaging. Hum Brain Mapp 8(1):13–27.

41. Forss N, Salmelin R, Hari R (1994) Comparison of somatosensory evoked fields toairpuff and electric stimuli. Electroencephalogr Clin Neurophysiol 92(6):510–517.

42. Gelnar PA, Krauss BR, Szeverenyi NM, Apkarian AV (1998) Fingertip representation inthe human somatosensory cortex: An fMRI study. Neuroimage 7(4 Pt 1):261–283.

43. Padberg J, Disbrow E, Krubitzer L (2005) The organization and connections of ante-rior and posterior parietal cortex in titi monkeys: Do New World monkeys have anarea 2? Cereb Cortex 15(12):1938–1963.

44. Hinkley LB, Krubitzer LA, Nagarajan SS, Disbrow EA (2007) Sensorimotor integrationin S2, PV, and parietal rostroventral areas of the human sylvian fissure. J Neurophysiol97(2):1288–1297.

45. van Kemenade BM, et al. (2014) Tactile and visual motion direction processing inhMT+/V5. Neuroimage 84:420–427.

46. Jiang F, BeauchampMS, Fine I (2015) Re-examining overlap between tactile and visualmotion responses within hMT+ and STS. Neuroimage 119:187–196.

47. Ricciardi E, Bonino D, Pellegrini S, Pietrini P (2014) Mind the blind brain to understandthe sighted one! Is there a supramodal cortical functional architecture? NeurosciBiobehav Rev 41:64–77.

48. Newsome WT, Wurtz RH, Komatsu H (1988) Relation of cortical areas MT and MST topursuit eye movements. II. Differentiation of retinal from extraretinal inputs.J Neurophysiol 60(2):604–620.

49. Ilg UJ, Schumann S (2007) Primate area MST-l is involved in the generation of goal-directed eye and hand movements. J Neurophysiol 97(1):761–771.

50. Jones EG, Powell TP (1969) Connexions of the somatic sensory cortex of the rhesusmonkey. I. Ipsilateral cortical connexions. Brain 92(3):477–502.

51. Jones EG, Powell TP (1969) Connexions of the somatic sensory cortex of the rhesusmonkey. II. Contralateral cortical connexions. Brain 92(4):717–730.

52. Almeida TF, Roizenblatt S, Tufik S (2004) Afferent pain pathways: A neuroanatomicalreview. Brain Res 1000(1-2):40–56.

53. Kalberlah C, Villringer A, Pleger B (2013) Dynamic causal modeling suggests serialprocessing of tactile vibratory stimuli in the human somatosensory cortex–an fMRIstudy. Neuroimage 74:164–171.

54. Khoshnejad M, Piché M, Saleh S, Duncan G, Rainville P (2014) Serial processing inprimary and secondary somatosensory cortex: A DCM analysis of human fMRI data inresponse to innocuous and noxious electrical stimulation. Neurosci Lett 577:83–88.

55. Chung YG, et al. (2014) Intra- and inter-hemispheric effective connectivity in thehuman somatosensory cortex during pressure stimulation. BMC Neurosci 15:43.

56. Liang M, Mouraux A, Iannetti GD (2011) Parallel processing of nociceptive and non-nociceptive somatosensory information in the human primary and secondary so-matosensory cortices: Evidence from dynamic causal modeling of functional magneticresonance imaging data. J Neurosci 31(24):8976–8985.

57. Barba C, Frot M, Mauguière F (2002) Early secondary somatosensory area (SII) SEPs.Data from intracerebral recordings in humans. Clin Neurophysiol 113(11):1778–1786.

58. Karhu J, Tesche CD (1999) Simultaneous early processing of sensory input in humanprimary (SI) and secondary (SII) somatosensory cortices. J Neurophysiol 81(5):2017–2025.

59. Hari R, et al. (1993) Functional organization of the human first and second somato-sensory cortices: A neuromagnetic study. Eur J Neurosci 5(6):724–734.

60. Simões C, Jensen O, Parkkonen L, Hari R (2003) Phase locking between human primaryand secondary somatosensory cortices. Proc Natl Acad Sci USA 100(5):2691–2694.

61. Mazzola L, Faillenot I, Barral F-G, Mauguière F, Peyron R (2012) Spatial segregationof somato-sensory and pain activations in the human operculo-insular cortex.Neuroimage 60(1):409–418.

62. David O, et al. (2013) Probabilistic functional tractography of the human cortex.Neuroimage 80:307–317.

63. Kurth F, Zilles K, Fox PT, Laird AR, Eickhoff SB (2010) A link between the systems:Functional differentiation and integration within the human insula revealed by meta-analysis. Brain Struct Funct 214(5-6):519–534.

64. Mesulam MM (1998) From sensation to cognition. Brain 121(Pt 6):1013–1052.65. Munari C, et al. (1994) Stereo-electroencephalography methodology: Advantages and

limits. Acta Neurol Scand Suppl 152(Suppl.c):56–67.66. Cossu M, et al. (2005) Stereoelectroencephalography in the presurgical evaluation of

focal epilepsy: A retrospective analysis of 215 procedures. Neurosurgery 57(4):706–718.67. Fedorov A, et al. (2012) 3D Slicer as an image computing platform for the Quanti-

tative Imaging Network. Magn Reson Imaging 30(9):1323–1341.68. Dale AM, Fischl B, Sereno MI (1999) Cortical surface-based analysis. I. Segmentation

and surface reconstruction. Neuroimage 9(2):179–194.69. Reuter M, Schmansky NJ, Rosas HD, Fischl B (2012) Within-subject template estimation

for unbiased longitudinal image analysis. Neuroimage 61(4):1402–1418.70. Van Essen DC, Drury HA (1997) Structural and functional analyses of human cerebral

cortex using a surface-based atlas. J Neurosci 17(18):7079–7102.71. Jenkinson M, Smith S (2001) A global optimisation method for robust affine regis-

tration of brain images. Med Image Anal 5(2):143–156.72. Abdollahi RO, et al. (2014) Correspondences between retinotopic areas and myelin

maps in human visual cortex. Neuroimage 99:509–524.73. Caruana F, et al. (2014) Human cortical activity evoked by gaze shift observation: An

intracranial EEG study. Hum Brain Mapp 35(4):1515–1528.74. Caruana F, Sartori I, Lo Russo G, Avanzini P (2014) Sequencing biological and physical

events affects specific frequency bands within the human premotor cortex: An in-tracerebral EEG study. PLoS One 9(1):e86384.

75. Rousseeuw PJ (1987) Silhouettes: A graphical aid to the interpretation and validationof cluster analysis. J Comput Appl Math 20:53–65.

76. Knutsen AK, et al. (2010) A newmethod to measure cortical growth in the developingbrain. J Biomech Eng 132(10):101004.

77. Van Essen DC (2012) Cortical cartography and Caret software. Neuroimage 62(2):757–764.

78. Jastorff J, Begliomini C, Fabbri-Destro M, Rizzolatti G, Orban GA (2010) Coding ob-served motor acts: Different organizational principles in the parietal and premotorcortex of humans. J Neurophysiol 104(1):128–140.

79. Geyer S, Schormann T, Mohlberg H, Zilles K (2000) Areas 3a, 3b, and 1 of humanprimary somatosensory cortex. Part 2. Spatial normalization to standard anatomicalspace. Neuroimage 11(6 Pt 1):684–696.

80. Amunts K, et al. (1999) Broca’s region revisited: Cytoarchitecture and intersubjectvariability. J Comp Neurol 412(2):319–341.

81. Tomassini V, et al. (2007) Diffusion-weighted imaging tractography-based parcella-tion of the human lateral premotor cortex identifies dorsal and ventral subregionswith anatomical and functional specializations. J Neurosci 27(38):10259–10269.

Avanzini et al. PNAS | Published online March 14, 2016 | E1943

NEU

ROSC

IENCE

PNASPL

US

Dow

nloa

ded

by g

uest

on

Janu

ary

24, 2

020