Previous reports conflict as to the role of primary somatosensory neocortex (SI) in tactile detection. We addressed this question in normal human subjects using whole-head magnetoencephalography (MEG) recording. We found that the evoked signal (0–175 ms) showed a prominent equivalent current dipole that localized to the anterior bank of the postcentral gyrus, area 3b of SI. The magnitude and timing of peaks in the SI waveform were stimulus amplitude dependent and predicted perception beginning at ∼70 ms after stimulus. To make a direct and principled connection between the SI waveform and underlying neural dynamics, we developed a biophysically realistic computational SI model that contained excitatory and inhibitory neurons in supragranular and infragranular layers. The SI evoked response was successfully reproduced from the intracellular currents in pyramidal neurons driven by a sequence of lamina-specific excitatory input, consisting of output from the granular layer (∼25 ms), exogenous input to the supragranular layers (∼70 ms), and a second wave of granular output (∼135 ms). The model also predicted that SI correlates of perception reflect stronger and shorter-latency supragranular and late granular drive during perceived trials. These findings strongly support the view that signatures of tactile detection are present in human SI and are mediated by local neural dynamics induced by lamina-specific synaptic drive. Furthermore, our model provides a biophysically realistic solution to the MEG signal and can predict the electrophysiological correlates of human perception.
- computational model
- dendritic processes
- conscious perception
- network dynamics
- somatosensory cortex
The role of primary sensory cortex in conscious perception is debated (Danckert and Goodale, 2000; Ress et al., 2000; Crick and Koch, 2003; Ergenoglu et al., 2004; Stoerig, 2006), and previous reports conflict as to the presence of neural correlates of tactile detection in primate primary somatosensory cortex (SI). In recent single-unit recording studies in monkeys, de Lafuente and Romo (2005, 2006) found that action potential firing rate in areas 3b and 1 did not predict “hit” and “miss” trials at threshold. Furthermore, monkeys can relearn a tactile detection task after comprehensive SI lesions (LaMotte and Mountcastle, 1975, 1979). In contrast, when somatosensory evoked potentials are recorded from the postcentral gyrus (PoCG) in macaques, the N1 and P2 peak magnitudes, at 50–65 and 105–130 ms, respectively, predict detection (Kulics, 1982; Kulics and Cauller, 1986; Cauller and Kulics, 1991). In human studies, response magnitude and coherence across electrodes above rolandic cortex have been reported to predict detection (Meador et al., 2002; Palva et al., 2005). Tactile spatial attention in humans also recruits what appear to be SI-specific effects (Drevets et al., 1995; Bauer et al., 2006). Furthermore, studies of the impact of ongoing activity at the time of stimulus presentation in humans show that functional magnetic resonance imaging (fMRI) blood–oxygen level-dependent signals localized to the anterior PoCG are greater on hit than miss trials (Moore et al., 2007), and intermediate amplitudes of prestimulus power (10, 20, and 40 Hz) in magnetoencephalography (MEG) signals measured above the sensorimotor cortex predict detection (Linkenkaer-Hansen et al., 2004).
These discrepancies may be attributable to differences in recording approach. Single-unit recordings are advantageous for reporting individual neuron activity but are difficult to obtain in humans, limiting observations to highly trained monkeys. Furthermore, action potential recordings do not provide information about subthreshold cellular and field potential activity or the activity of smaller neurons that may be coincident with sensory information processing. Surface recordings using MEG or EEG provide access to neurophysiological signals in human subjects, but the cortical population dynamics inducing the recorded signal have not been conclusively defined. Furthermore, these studies often do not localize activity to specific cortical areas.
We examined cortical correlates of somatosensory perception in human SI by combining MEG and computational neural modeling. By calculating the equivalent current dipole (ECD) at SI, we were able to localize an early and robust source to the anterior PoCG, putative area 3b. The timing and magnitude of this SI waveform predicted perception of the tactile stimulus beginning at 70 ms after stimulus.
To provide a framework for understanding the MEG signal and SI neural correlates of perception, we developed a realistic laminar network model. Longitudinal intracellular currents were sampled from detailed compartmentalized pyramidal neurons (PNs) in the supragranular and infragranular layers to calculate the net ECD produced by the population. We found that the timing, magnitude, and direction of the peaks in the observed SI evoked response could be accurately reproduced in the model when a sequence of lamina-specific exogenous drive to the SI network was simulated. The model also predicted that the observed differences in SI between detected and nondetected stimuli resulted from subtle changes in the timing and amplitude of an ∼70 ms input to the supragranular layers and a later ∼135 ms output from the granular layers. These findings provide evidence that SI activity predicts tactile detection and that these predictive signals in SI depend on the pattern of exogenous, lamina-specific input.
Materials and Methods
Seven neurologically healthy, right-handed, 18- to 45-year-old adults were studied.
Brief taps were delivered to the subject's right hand in the form of a single cycle of 100 Hz sine wave (10 ms duration) via a custom piezoelectric device. Subjects rested their hand on a Delrin frame that held a piezoelectric parallel to the finger (Noliac ceramic multilayer bender plate 32 × 7.8 × 1.88 mm). A deflection stroke drove a Delrin contactor into the fingertip (7 mm diameter contractor presented within a 1 cm circular rigid surround).
Individual subjects' thresholds were obtained before imaging using a parameter estimation by sequential testing (PEST) convergence procedure (Dai, 1995; Leek, 2001). This algorithm obtained the threshold by first using a strong stimulus (100% amplitude, 350 μm deflection), followed by a weak stimulus (50% amplitude) and a blank stimulus. If the subject reported detection of the strong and the weak stimulus, then the next set of stimuli was given the 50% value set for the strong stimulus and a new value at one-half the distance between the weak and blank stimulus for the middle stimulus. If the subject did not detect the middle stimulus, then the maximum stimulus remained the same, the middle stimulus became the new minimum stimulus, and the new value for the middle stimulus was obtained by averaging the maximum and new minimum. This procedure was repeated until the change in the amplitude of movement for the piezoelectric between trials was <5 μm. Before MEG recording, the stimulus strength was set to the threshold value obtained during the PEST.
During MEG imaging, for 70% of presented trials, stimulus strength was maintained at a perceptual threshold (T) level (50% detection) using a dynamic algorithm. If two correct responses were made, the threshold-level voltage sent to the piezoelectric was decreased by 0.005 V (a change of ∼4.5 μm in piezoelectric movement), and the correct response counts were reset to zero. In contrast, if three incorrect responses were made, the voltage was increased by 0.005 V. The deflection amplitudes for all threshold stimuli across all subjects were between 200 and 230 μm. Suprathreshold (ST) stimuli (10% of all trials; 350 μm deflection; 100% detection) and null trials (20%) were randomly interleaved with the threshold stimuli and were excluded from the correct response counts used in the dynamic algorithm. Trial duration was 3 s. Each subject underwent eight runs with 120 trials. Trial onset was indicated by a 60 dB, 2 kHz auditory cue delivered to both ears for a duration of 2 s. During the auditory cue, the 10 ms finger tap stimulus was delivered between 500 and 1500 ms, in 100 ms intervals, from trial onset. The number of trials of a given latency to tap was randomly distributed during each run. After the cessation of the auditory cue, subjects reported detection or nondetection of the stimulus with button presses using the second and third digits of the left hand, respectively. The auditory cue ended ≥500 ms after tactile stimulus and 1000 ms before the next trial began.
MEG data acquisition.
Using a 306-channel magnetoencephalograph (VectorView; Elekta NeuroMag, Helsinki, Finland), neuromagnetic responses were recorded with 306 sensors arranged in triplets of two planar gradiometers and a magnetometer at 102 sites. In addition to MEG, the vertical and horizontal electrooculogram (EOG) was recorded with electrodes placed close to the left eye. Four head-position indicator coils were placed on the subject's head to coregister the subject's anatomical MRI and the MEG sensors. The data were sampled at 600 Hz with the bandpass set to 0.01–200 Hz. The responses were averaged on-line for quality control. In the off-line analysis, the data were reaveraged using a bandpass of 0.03–200 Hz. Epochs with EOG peak-to-peak amplitude exceeding 150 μV were excluded from the analysis.
MEG source analysis.
The aim of our source analysis was to locate a primary current-dipole source at the contralateral SI area and to find the time course of this source, taking into account the presence of other active areas. Previous studies have shown that somatosensory evoked fields (SEFs) can be well approximated with a small number (typically three, for median nerve stimulation) of serially and simultaneously active ECD sources (Brenner et al., 1978; Hari and Forss, 1999). For tactile stimulation, a primary current source is typically initially observed in SI with later activation having additional contributions from a second source, likely representing SII (Forss et al., 1994b; Hari and Forss, 1999; Hoechstetter et al., 2000, 2001; Kakigi et al., 2000). To isolate the contribution from SI to the SEF, we used the following approach.
Initial inspection of the early field patterns of the suprathreshold trials demonstrated activity consistent with SI activation contralateral to the stimulus, and later field patterns indicated a source in the contralateral parietal operculum, presumably SII. The activation of this second source was weak and inconsistent across subjects. No consistent activity over the ipsilateral SII or in other brain areas was observed. Therefore, we modeled the data with two dipoles and optimized this fit with help of the signal-space projection (SSP) method (Tesche et al., 1995; Uusitalo and Ilmoniemi, 1997) as follows. First, using a least-squares fit with the dipole forward solution calculated using the spherically symmetric conductor model (Sarvas, 1987; Hämäläinen and Sarvas, 1989), we found an initial ECD at the peak activity in the suprathreshold stimulus signals from one data run (average of 12 trials; mean, 68 ms; SD, 8 ms). With help of the anatomical MRI, coregistered with the MEG, we were able to confirm that this source localized to SI in four of seven subjects (Fig. 1 A, B). In the second step, the contribution of the SI ECD was removed from the data using the SSP method, and a second ECD was fitted to the residual data (Tesche et al., 1995; Uusitalo and Ilmoniemi, 1997; Nishitani and Hari, 2000). Importantly, the same projection operator that was applied to the data to remove the effect of the SI source was also taken into account in the forward solution when the second dipole was fitted (Uusitalo and Ilmoniemi, 1997). The goodness of fit of the two-dipole model was larger than 70% in all fit data during peak responses. In the final step, the effect of the second ECD, which was not the focus of the present study, was removed from the data using SSP, the SI ECD was refitted to the residual, and its waveform was recalculated. The ECD localizations fitted to the suprathreshold data were used to model all responses. To account for adaptation and learning effects with training, only the last 100 trials for a given response were considered for analysis.
For two of seven subjects, large baseline rhythmic activity interfered with the localization of peak responses. For these subjects, the initial source was placed in the anterior bank of the postcentral gyrus in an area consistent with the localization of the finger representation of area 3b (Fig. 1 C) (Penfield and Rasmusson, 1950; Uematsu et al., 1992; White et al., 1997; Yousry et al., 1997; Sastre-Janer et al., 1998; Moore et al., 2000), and the second dipole source was placed in the parietal operculum. For one subject, we were unable to obtain anatomical MRI data. Dipole localization was determined by field contours on the spherical head model and showed an initial source above the predicted position of contralateral SI and a second source over the predicted position of the parietal operculum. The source localizations used for each of the three subjects discussed above produced evoked dipole waveforms consistent with those of the other subjects (see Results).
Computational neural model
Several lines of evidence predict that postsynaptic intracellular longitudinal currents within the long apical dendrites of synchronized cortical pyramidal cells are the main contributors to MEG primary current sources (Hämäläinen et al., 1993; Okada et al., 1997; Murakami et al., 2002, 2003; Ikeda et al., 2005; Murakami and Okada, 2006). As such, we incorporated realistic pyramidal neuron morphology and physiology into our SI model and calculated the net ECD from the longitudinal intracellular currents within these pyramidal neurons. The complete SI network contained pyramidal neurons and inhibitory interneurons in the supragranular and infragranular layers, as depicted in Figure 2 A. This laminar construction enabled us to study the influence of a specific physiologically based sequence of exogenous synaptic inputs, defined by the laminar location of their postsynaptic effects (Figs. 2 C,D), on the ECD produced by the SI population.
The individual neurons and network architecture were implemented as follows.
Pyramidal neuron morphology and physiology.
The model contained PNs with somata in layers 2/3 (L2/3) and layer 5 (L5) (Fig. 2). Simulations of L2/3 and L5 PN morphology and physiology were adapted from Bush and Sejnowski (1993), whose code is available in NEURON software format at http://senselab.med.yale.edu/senselab/modeldb/. Individual L2/3 and L5 PNs were constructed with a small number of compartments (eight and nine, respectively) and maintain the morphology of the real PNs on which they are based [digitized HRP-filled L2 and L5 PNs from the cat visual cortex (Koch et al., 1990)]. The accurate morphological construction of the dendrites allowed for a spatially accurate network model. The compartmentalization of each of the PNs is shown in Fig. 2 A and described in Table 1. In our model, a scaling factor of 1.3 was applied to the length and diameters of the dendritic compartments used by Bush and Sejnowski (1993) to account for increases in dendritic length and volume in human somatosensory neurons, as predicted by larger cortical thickness (Geyer et al., 1997; Fischl and Dale, 2000) and an increase in the number of dendritic spines and arborization (Elston et al., 2001). The membrane resistance was increased and membrane capacitance was decreased by the same scaling factor (R m = 23,474 Ω · cm2 for L5 and L2/3; C m = 0.85 and C m = 0.6195 μF/cm2 for L5 and L2/3, respectively) to maintain the input resistances in the cells of 45 MΩ for the L5 and 110 MΩ for L2/3 (Douglas et al., 1991). The axial resistance for each cell was R a = 200 Ω · cm (Segev et al., 1992).
Active currents in L2/3 PNs included a fast sodium current (I Na), a delayed rectifier potassium (I Kdr) current, an adapting potassium current (I M), and a leak current (I L). The L5 PNs contained the same currents with the addition of a calcium current (I Ca) and a potassium-activated calcium current (I KCa). The kinetic equations and NEURON code used for each of these currents were as used by Mainen and Sejnowski (1996) and downloaded from http://senselab.med.yale.edu/senselab/modeldb/. The maximal conductances of each current were constant throughout the soma and dendrite (Stuart and Sakmann, 1994; Bekkers, 2000; Korngreen and Sakmann, 2000; Migliore and Shepherd, 2002) and were chosen to produce adapting spikes in the L2/3 PNs and bursting in the L5 PNs to current injected in the soma (1 nA for 100 ms) (Fig. 2 B) representative of neurons classified as regular spiking and intrinsically bursting, respectively (Silva et al., 1991; Moore and Nelson, 1998; Zhu and Connors, 1999). The maximal conductance parameters (in S/cm2) were as follows: L2/3 parameters, g Na = 0.15; g K = 0.01; g M = 0.00025; g L = 0.0000426; L5 parameters, g Na = 0.14; g K = 0.01; g M = 0.0002; g Ca = 0.00006; calcium decay time constant = 20 ms; g KCa = 0.2 × 10−9. The ionic reversal potentials (in mV) were as follows: E Na = −50; E K = E M = 77; E L = −65.
Inhibitory interneurons (INs) were included to simulate their postsynaptic effect on the network. They were modeled with single compartments and contained only fast sodium (I Na) and potassium currents (I Kdr) to create spiking activity, as in other network models (Jones et al., 2000; Garabedian et al., 2003; Pinto et al., 2003). Parameters regulating the IN dynamics were length = 39 μm; diameter = 20 μm; R a = 200 Ω · cm; C m = 85 μF/cm2; g Na = 0.12 S/cm2; g K = 0.036 S/cm2; g L = 0.003 S/cm2; E Na = −50 mV; E K = 77 mV; E L = −54.3 mV.
Local synaptic architecture.
Each modeled cortical layer contained 10 PNs and three INs (Thomson et al., 2002). Local excitatory and inhibitory synapses within the SI cortical column model were constructed as shown in Figure 2 A. Connection lines are schematic representations of axonal-to-dendritic input. Axons were not explicitly modeled. The local synaptic architecture was based on an abundance of animal studies and, in particular, studies of the mouse/rat somatosensory cortex (Bernardo et al., 1990a,b) [for review, see Thomson et al. (2002), Thomson and Bannister (2003), and Bannister (2005)]. Inhibitory synaptic connections onto PNs were located on the soma (Somogyi et al., 1983; Kisvarday et al., 1985; Freund et al., 1986), and excitatory synapses contacted the basal and apical oblique dendrites (Deuchars et al., 1994; Lubke et al., 1996; Thomson and Bannister, 1998; Feldmeyer et al., 2002). Fast and slow excitatory (AMPA/NMDA) and inhibitory (GABAA/GABAB) synapses were simulated using an α function that was “turned on” by the soma of the presynaptic cell crossing a voltage threshold (0 mV). The synaptic dynamics were defined by the following rise/decay time constants and reversal potentials, respectively: AMPA, 0.5/5 ms, 0 mV; NMDA, 1/20 ms, 0 mV; GABAA, 0.5/5 ms, −80 mV; GABAB, 1/20 ms, −80 mV. The strengths of the synaptic connections within the local network were defined with a Gaussian spatial profile, with a delay incorporated into the synaptic connection between two cells defined by an inverse Gaussian (Jones, 1986; Kaas and Garraghty, 1991). The maximum synaptic conductances and Gaussian weight space constants (number of cells from center) are listed in Table 2 along with the minimum synaptic delay and corresponding Gaussian delay space constant. The three INs in each layer were regularly distributed in space among the 10 PNs.
Exogenous drive to the local network was excitatory only (Cauller and Connors, 1994; Cauller et al., 1998; Guillery and Sherman, 2002) and was defined by the laminar location in SI of its synaptic effects based on general principles of cortical circuitry (Rockland and Pandya, 1979; Friedman and Jones, 1980; Felleman and Van Essen, 1991; Jones, 2001; Douglas and Martin, 2004). One source of drive emerged from the granular layer, layer 4 (L4), and contacted the L2/3 neurons, with a delayed and weaker connection to the infragranular L5 neurons (for specific poststimulus dendritic compartments, see Fig. 2 C). Activity in L4 is modeled to reflect drive from the thalamus and was based on several studies of intracranial laminar electrophysiological recordings of evoked responses in SI, including responses to vibrissa and thalamic stimuli in rodents (Di et al., 1990; Barth and Di, 1991; Castro-Alamancos and Connors, 1996; Kandel and Buzsáki, 1997; Douglas and Martin, 2004), trigeminal stimulation in piglets (Ikeda et al., 2005), and tactile (Kulics and Cauller, 1986; Cauller and Kulics, 1991) and median nerve stimuli in awake monkeys (Peterson et al., 1995; Lipton et al., 2006). A second source of drive to the SI network contacted the distal apical dendrites in the supragranular layers of each neuronal population (Fig. 2 D). This connection could be representative of input from higher-order cortical areas or nonspecific thalamic sources (Rockland and Pandya, 1979; Friedman et al., 1980; Felleman and Van Essen, 1991; Jackson and Cauller, 1998; Jones, 2001; Douglas and Martin, 2004).
The sources of SI drive were modeled as spike generators with a predefined temporal profile, such that on a trial-by-trial (n = 100 trials) basis the presynaptic spike times were chosen from a Gaussian distribution to reproduce a physiologically realistic sequence of input to SI as follows: initial granular layer (GR) output, mean, 25 ms; SD, 2.5 ms; subsequent exogenous input to the supragranular layer (SGR), mean, 70 ms; SD, 6 ms; second later wave of GR output, mean, 135 ms; SD, 7 ms (Table 3) (all GR to L5 inputs had an additional fixed 5 ms conduction delay). The mean input timings in this sequence were chosen to be consistent with laminar recordings of electrical activity in SI during tactile (Kulics and Cauller, 1986; Cauller and Kulics, 1991; M. L. Lipton and C. E. Schroeder, personal communication), thalamic (Kandel and Buzsáki, 1997), and trigeminal (Ikeda et al., 2005) stimulation. The spatial distribution of the synaptic weights were uniform, as in Table 3. A baseline noise level was incorporated into the model by injecting a random amount of current into each neuronal compartment at each time step taken from a uniform distribution between −0.3 and 0.3 nA. This random noise was included to create heterogeneity in spike timing across the population and trials.
In each simulation presented, the biophysical properties of the cortical neurons, the intrinsic and exogenous synaptic connectivity, and the sequence in which exogenous drive arrived to the model (GR-SRG-GR) were fixed, based on previous physiological and anatomical findings. Our initial simulations with the model gave close to accurate matching with the MEG data, and parameter fitting was done empirically on only the relative strength and/or relative timing of the exogenous drives (given the fixed overall sequence) to best fit the MEG data. Specifically, we manually adjusted these parameters until our quantitative estimate of the current dipole moment (given in nA · m) matched those observed experimentally when multiplied by a constant network scaling factor of 3000. This scaling factor was fixed for each simulation.
Calculation of net current dipole.
The SI ECD was calculated as the net sum across the population of the intracellular currents flowing within the PN dendrites in a direction perpendicular to the longitudinal axis of the apical dendrite multiplied by the corresponding length of the dendrite.
More simplified models, often referred to as “neural mass models,” have been used to investigate signals underlying event-related MEG/EEG responses in the human brain (David et al., 2005, 2006a,b; Lee et al., 2006; Riera et al., 2006, 2007). These models generally assume that the MEG signal is proportional to some superposition of the model neurons' membrane potentials, which contribute to the mean state, and have succeeded in capturing basic features of MEG evoked responses. Consistencies between such simplified models and those that rely on details of the PN dendritic morphology and physiology remain to be investigated. However, our results and those from more reduced preparations (Murakami et al., 2002, 2003; Murakami and Okada, 2006) suggest that these nonspecific models will obscure predictions about primary ECD sources that depend crucially on intracellular dendritic currents.
All simulations were performed using the shareware software program NEURON available at http://www.neuron.yale.edu/neuron/. A fixed-time-step implicit Euler integration method was used with a time increment dt = 0.025 ms. Results shown were smoothed with a 333 ms Hamming window. After publication, the code that produced all simulated data in this paper will be available on the ModelDB website, http://senselab.med.yale.edu/senselab/modeldb/, and we refer the reader here for equations regulating active current kinetics.
Features of SI evoked responses
A brief suprathreshold stimulus applied to the D3 digit tip evoked a robust and consistent response above the contralateral somatosensory cortex. The peak activity (∼70 ms) of the estimated SI ECD reliably localized to the postcentral gyrus in the hand area of SI, consistent with other studies of tactile stimulation to the fingers (Forss et al., 1994b; Hoechstetter et al., 2000, 2001; Kakigi et al., 2000; Braun et al., 2002; Druschky et al., 2003; Iguchi et al., 2005). Figure 1 shows representative examples of the location and orientation of the ECD on T1-weighted MRI images. Figure 1, A and B, is derived from the inverse solution technique described in Materials and Methods, and Figure 1 C shows an example of one of two subjects for which the ECD was placed in area 3b. In each subject, the dipole was localized in the lateral edge of the hand area defined by the “Ω”-shaped passage of the central sulcus. Previous studies have shown that this position is consistent with the location of the finger representation in humans (Penfield and Rasmusson, 1950; Uematsu et al., 1992; White et al., 1997; Yousry et al., 1997; Sastre-Janer et al., 1998; Moore et al., 2000).
The SI ECD responses evoked by suprathreshold- and threshold-level stimuli showed reliable peaks with consistent polarity in all subjects. Figure 3 displays the SI response in each individual subject for suprathreshold- (Fig. 3 A, left) and threshold (Fig. 3 A, right)-level stimuli, averages across all subjects (Fig. 3 B), and individual and mean differences between suprathreshold- and threshold-level responses (Fig. 3 C). The average response to the suprathreshold stimulus (Fig. 3 B, red curve) had an initial peak with positive polarity at ∼25 ms (labeled M25), followed by a peak with negative polarity at ∼35 ms (M35) and another smaller but consistent positive peak at ∼50 ms (M50). The largest peak occurred at ∼70 ms (M70) with a negative polarity, followed by two peaks with positive polarity at ∼100 ms (M100) and ∼135 ms (M135) (Forss et al., 1994b; Hoechstetter et al., 2001; Druschky et al., 2003).
When comparing suprathreshold- and threshold-level responses, several features of the evoked response were stimulus amplitude dependent. The amplitude of the dynamic threshold-level stimulus was always at least 120 μm smaller than the fixed stimulus amplitude of the suprathreshold stimulus (see Materials and Methods). The earliest peak evident in the threshold-level response (Fig. 3 B, blue curve) was the M70. The M25-M35-M50 peaks in the threshold-level response are likely a result of the lack of sufficient signal-to-noise to detect the subtle early signal from the weak tactile stimulus. Although the observed M70 peak in the threshold-level response had the same negative polarity of the M70 peak of the suprathreshold-level response, the magnitude was smaller (M70 = minimum in 50–100 ms, paired t test, p < 0.001; ST, −130 ± 46 nA · m; T, −58 ± 27 nA · m) and the latency was longer (p = 0.0661; ST, 72 ± 9 ms; T, 80 ± 6 ms), resulting in a smaller onset slope to peak (slope from response at 50 ms to M70, p < 0.0002; ST, −6 ± 3 nA · m/ms; T, −2 ± 2 nA · m/ms). Two subsequent peaks with positive polarity, M100 and M135, were also evident in the threshold response. The magnitudes of these peaks and the area under the curve between them were smaller for the threshold-level response (M100 = maximum in 75–125 ms, p < 0.02; ST, 110 ± 80 nA · m; T, 26 ± 23 nA · m; M135 = maximum 125–175 ms, p < 0.02; ST, 127 ± 71 nA · m; T, 47 ± 30 nA · m; average magnitude, 100–150 ms, p < 0.04; ST, 82 ± 68 nA · m; T, 18 ± 19 nA · m). Furthermore, the latency of the M135 was longer (p < 0.02; ST, 140 ± 9 ms; T, 155 ± 8 ms) for the threshold-level response. The latency of the M100 peak was not significantly different for the threshold-level response (p = 0.44; ST, 113 ± 9 ms; T, 116 ± 7 ms). However, the onset slope to the M100 peak was smaller (slope from M70 to response at 100 ms, p < 0.02; ST, 7 ± 4 nA · m/ms; T, 3 ± 2 nA · m/ms).
Figure 3 C shows that the stimulus amplitude-dependent differences were not driven by the response of an “outlier” subject. Displayed are the differences between the suprathreshold- and threshold-level responses for each subject (thin colored curves), along with the mean difference across all subjects (thick black curve), the mean excluding the subject with the smallest difference (thick dark gray curve; excludes red curve), and the mean excluding the subject with the largest mean differences (thick light gray curve; excludes cyan curve). The excluded subjects were based on the mean absolute difference from 0 to 175 ms (smallest = 23.71 nA · m; median = 49.65 nA · m; largest = 96.2 nA · m).
The SI evoked response predicts perception
We used an adaptive stimulus algorithm, such that the threshold-level stimulus amplitude was maintained dynamically at ∼50% detection rate (Dai, 1995; Leek, 2001). Because misses [i.e., nonperceived (NP) trials] in this design will, on average, occur for lower-amplitude stimuli, differences observed in stimulus response for perceived (P) and nonperceived trials could simply reflect peripheral input amplitude. Therefore, we analyzed the correlation between the SI evoked response and perception in two ways, each of which showed statistically significant differences between perceived and nonperceived trials (Fig. 4 A,B).
In the first analysis, we averaged the last 100 trials of perceived and nonperceived stimuli to minimize within-session training effects (Fig. 4 A) (n = 7). This dataset provided a strong statistical sample and showed that the onset slope from the M70 to the M100 peak was larger for perceived (dark blue) than for nonperceived (light blue) trials (slope from 70 to 100 ms, p < 0.002; P, 2 ± 2 nA · m/ms; NP, 0.5 ± 2 nA · m/ms), and the magnitudes of the M100 and M135 peaks and area under the curve between them were larger (M100 = maximum in 95–110 ms, p < 0.04; P, 31 ± 27 nA · m; NP, 1 ± 38 nA · m; M135 = maximum 120–175 ms, p < 0.02; P, 68 ± 26 nA · m; NP, 43 ± 32 nA · m; average magnitude, 100–150 ms, p < 0.02; P, 32 ± 22 nA · m; NP, 1 ± 20 nA · m). Other trends that did not reach significance were present in SI. These included shorter latencies to the M70 and M135 peaks for the perceived stimuli (M70, p = 0.0918; P, 74 ± 11 ms; NP, 84 ± 10 ms; M135, p = 0.1253; P, 143 ± 16 ms; NP, 154 ± 12 ms).
In the second analysis, for each subject, we subsampled the dataset and analyzed an equal number of hit and miss trials from an equal stimulus amplitude (Fig. 4 B). The amplitude selected for analysis was the level that was most commonly encountered by a given subject (mean number of trials, 19; SD, 9; n = 6 subjects; for one subject, the data file containing stimulus amplitudes was corrupted). The signals generated during perceived and nonperceived trials in this “equal” condition also showed significant differences between hit and miss trials. The mean response for perceived trials showed a well defined M70 peak (minimum, 70–145 ms; magnitude, −118 ± 89 nA · m; latency, 90 ± 5 ms), whereas nonperceived trials did not. When comparing the M70 peak for perceived trials with the temporally aligned 90 ms (mean latency for perceived trials) response for nonperceived trials, the magnitude was larger for perceived trials (M70 magnitude, p < 0.05; P, −118 ± 89 nA · m; NP, −6 ± 33 nA · m), as was the onset slope to the M100 (slope M70 to M100 = magnitude at 125 ms, p < 0.05; P, 5 ± 4 nA · m/ms; NP, 0.003 ± 2.101 nA · m/ms). The significant differences in the M100 and M135 magnitude, and area under the curve between them, that emerged in the larger dataset showed a trend, but did not reach statistical significance, in the equal condition analysis (M100 = maximum at 125 ms, p = 0.1151; P, 51 ± 80 nA · m; NP, −9 ± 42 nA · m; M135 = maximum 125–175 ms; magnitude, p = 0.0787; P, 123 ± 72 nA · m; NP, 70 ± 33 nA · m; average magnitude, 100–150 ms, p = 0.146; P, 27 ± 52 nA · m; NP, −9 ± 15 nA · m).
In a third analysis (Fig. 4 C), we sorted the threshold-level responses by stimulus amplitude to test whether threshold-level stimulus amplitude was a better predictor of response magnitude. We found that the average responses to the largest (L) and smallest (S) threshold-level stimulus amplitudes (Fig. 4 C, dark and light green curves, respectively) (100 trials each; n = 6 subjects) showed no statistically significant difference at any time in the 0–175 ms response (M70 magnitude = minimum 50–100 ms, p = 0.08; L, −81 ± 51 nA · m; S, −65 ± 39 nA · m; M70 latency, p = 0.5; L, 81 ± 13 ms; S, 76 ± 14 ms; slope M70 to M100 = maximum 95–110 ms, p = 0.16; L, 4 ± 3 nA · m/ms; S, 3 ± 3 nA · m/ms; M100 magnitude, p = 0.9; L, 21 ± 21 nA · m; S, 20 ± 46 nA · m; M135 magnitude = maximum 120–175 ms, p = 0.5; L, 62 ± 36 nA · m; S, 53 ± 31 nA · m; M135 latency, p = 0.4; L, 152 ± 9 ms; S, 147 ± 16 ms; average magnitude, 100–150 ms, p = 0.11; L, 25 ± 16 nA · m; S, 14 ± 21 nA · m). This analysis shows that detection probability is a better predictor of variance in the SI response than peripheral stimulus amplitude.
Computational neural modeling
The ECD inverse modeling technique that we used allowed us to infer the magnitude and direction of the net current flow from a dipole source within a local cortical region. In the observed SI ECD responses, positive polarity corresponds to net current flow anterior, which in the case of area 3b corresponds to flow toward the cortical surface, and negative polarity corresponds to net current flow toward the white matter. Confidence that the early ECD signal arises from somatosensory area 3b, without contribution from other neighboring somatosensory areas, comes from the fact that, unlike the surrounding somatosensory areas 1 and 3a, the orientation of area 3b in the postcentral gyrus is tangential to the nearby inner skull surface, which makes it particularly well situated for MEG imaging (Hämäläinen et al., 1993). The magnitude of the ECD is commonly believed to result from the net current flow in the dendrites of a synchronous population of PNs, with the polarity determined by the direction of this current flow (Hari et al., 1980; Hämäläinen et al., 1993; Okada et al., 1997; Ikeda et al., 2005). Thus, to understand the local cortical dynamics that induced the magnitude and polarity of the observed ECD responses, we developed a biophysically realistic laminar cortical model of a local SI network, representing area 3b. Simulations of net current flow within PNs in the model were used to make predictions on the underlying network dynamics that define the evoked responses and their relation to perception.
A schematic of the local SI network is shown in Figure 2 A. The model consisted of excitatory PNs and INs in L2/3 and L5, shown in light green and red, respectively. To maintain accurate morphology, the PNs were modeled with multiple dendritic compartments, and they contained active conductances in their soma and dendrites to reproduce realistic spiking patterns to somatic injected current (1 nA) (Fig. 2 B). The INs were represented as single compartment spiking neurons in each layer and were included to simulate their postsynaptic effects in the network. Synaptic connections in the model included local synapses (Fig. 2 A), and exogenous drive simulated as (1) output from activity in the GR, representative of activity induced by the lemniscal thalamus (Fig. 2 C), and (2) input to the SGR, representative of input from higher-order cortex or nonspecific thalamic sources (Fig. 2 D). The GR output and input to the SGR were each modeled as a spike train generator that consisted of a single spike on each trial with a mean latency across trials (n = 100) taken from a Gaussian distribution, as described in Materials and Methods and Table 3.
All parameters, including intrinsic currents, connectivity, and cell morphology, were derived from published experimental and modeling studies of mammalian neocortex, primarily from the somatosensory system (for relevant literature, see Materials and Methods and Discussion). The model construction enabled us to investigate the influence of a physiological sequence of exogenous drive on the ECD produced by the local SI network. The relative timing and strength of the exogenous drives were the only “free” parameters in the model that were adjusted to accurately simulate the MEG evoked signal.
Polarized current flow driven by a sequence of exogenous drive to SI predicts the temporal sequence of the evoked response
Having set the intrinsic network parameters and synaptic architecture in the model, we reproduced the evoked response by simulating a specific sequence of excitatory drive to the SI model (Fig. 5). The sequence consisted of output from activity in the GR emerging at ∼25 ms, followed by input to the SGR at ∼70 ms and a subsequent late wave of GR output at ∼135 ms. Each input time was chosen from a Gaussian distribution across trials (Table 3), as shown schematically in green in Figure 5 A with an arrow indicating the mean input times. This pattern of synaptic drive could be interpreted as initial feedforward input arriving from the lemniscal thalamus to the GR, followed by feedback input from a higher-order cortical area or nonspecific thalamic neurons to the SGR, followed by a late lemniscal thalamic input to the GR induced as part of a thalamocortical loop of activity (see Discussion). Figure 5 A (red curve) shows the net current dipole from the PN population resulting from this pattern of synaptic drive averaged over 100 simulations, with suprathreshold-level synaptic conductances as in Table 3. The feedforward output from the GR at ∼25 ms (Gaussian mean, 25 ms; SD, 2.5 ms) induced the entire initial dipole volley containing the M25-M35-M50 peaks. We describe the origin of each induced peak separately.
The separate contribution to the population average suprathreshold response from L2/3 (gray) and L5 (black) PN populations are shown in Figure 5 B (top) along with representative contributions from a single L2/3 and L5 PN during a single trial (bottom). As shown in Figure 5 B (top), the M25 peak was primarily induced by the activity of the L2/3 PNs because the GR output to L5 PNs was weaker and delayed (Bannister, 2005). The 25 ms input to the L2/3 PNs was strong enough to induce somatic spikes that back-propagated into the apical and basal dendrites, as shown in the example neuron in Figure 5 B (bottom left) (Stuart and Sakmann, 1994; Murakami et al., 2002). The contributions from the apical (red) and basal (green) dendrites to the total dipole current (blue) produced by a single L2/3 PN are shown separately, as well as the somatic membrane potential (black). The back-propagation of current up the apical dendrites in the L2/3 PNs created a net current flow in the cell, and in the nearly synchronous population, that was directed upward, and hence the net dipole had positive polarity at the M25 peak (Fig. 5 A, red curve). The timing and duration of the M25 peak were set by the kinetics of active currents in the L2/3 PN dendrites, the slight heterogeneity in spike timing across the population, and the Gaussian distribution of GR output times across trials.
A subsequent repolarization of the L2/3 dendrites created a net downward current flow and induced the negative M35 peak (Fig. 5 B, top, gray curve; bottom left, red curve). The absolute magnitude of this peak was diminished by the initial current flow in the L5 PNs that was delayed and opposite in direction (Fig. 5 B, top, black curve). Although the drive to the L5 neurons was not strong enough to create spiking, it did induce currents that back-propagated up the L5 dendrites creating a net positive current dipole in the population. The duration of this peak was again influenced by the kinetics of the L2/3 PN active dendritic currents, as well as the time course of inhibitory synaptic activity near the soma.
The subsequent positive M50 peak was induced by a combination of continued synaptic activity within the L2/3 network from the initial drive, via the local network connections, and the residual excitatory synaptic activity in the L5 PNs.
The large negative M70 peak (Fig. 5 A, red curve) was reproduced by simulating excitatory drive to the SGR at ∼70 ms (Gaussian mean, 70 ms; SD, 6 ms). This SGR input induced downward currents in the apical dendrites of each PN population (Fig. 5 B, top). Although the strength of input was the same to each population, the activity induced in the L5 PNs dominated the resulting net dipole because of the larger length and diameter of their dendrites (Fig. 5 B, top, black trace). The strength of the SGR input was sufficient to induce bursting activity in the somata of the L5 PNs. An initial downward current in the apical dendrites that lasted ∼10 ms preceded each spike in the burst (Fig. 5 B, bottom right). These downward currents created a net negative response across the population and hence the negative M70 peak. The duration of this peak was again set by a combination of active dendritic current dynamics, slight heterogeneity in burst timing across the population, and the Gaussian distribution of feedback input timing across trials.
After each spike in the burst, there was an active back-propagation of current up the apical dendrites of the L5 PNs (Fig. 5 B, bottom right), and this upward current flow created the positive M100 peak in the average dipole. The duration of the bursting activity in the L5 PNs was an essential feature in defining the duration of the M100 peak activity (Fig. 5 A, red curve).
After the M100 peak activity, the simulated dipole current returned to baseline, creating an apparent discrepancy between the simulated evoked dipole (Fig. 5 A) and the observed average evoked dipole (Fig. 3 B), which fell toward, but did not completely return to, baseline after the M100 peak. This inconsistency arose from the fact that in our model the bursting activity of the L5 PNs ended nearly synchronously across the population. We predict that simulations of a larger network model that incorporates more heterogeneity and perhaps subclusters of synchronous networks would rectify this discrepancy.
To reproduce the positive M135 peak, a late GR output at ∼135 ms (Gaussian mean, 135 ms; SD, 7 ms) was simulated. This drive acted in conjunction with the low level of ongoing spiking network activity to create a more synchronous activity profile. Again this spiking created a dominant active back-propagation of current up the apical dendrites of the PNs, and hence a net current dipole with positive polarity (Fig. 5 A, red curve).
The arrows located next to the PNs in the network schematics in Figure 5 B (top; similarly in Fig. 5 C) emphasize the fact that positive polarity in the SI response was created by a net intracellular current flow up the apical dendrites of the PNs induced by back-propagation of spiking activity (M25, M100, and M135) and excitatory synaptic inputs on the basal and oblique dendrites (M50), whereas negative polarity was created by a net intracellular current flow down the apical dendrites of the PNs induced by postspike repolarization of the apical dendrites (M35) and excitatory synaptic input to the most distal apical dendrites (M70). To accurately reproduce the magnitude of the observed MEG dipole peaks, which were on the order of 50–100 nA · m for the suprathreshold-level stimulus, the model results were multiplied by a scaling factor of 3000 (Fig. 5 A,B). Thus, because there were 10 PNs in L2/3 and L5, the model predicts that ∼60,000 PNs contribute to the net dipole created from the brief tactile stimulus used. This prediction is in agreement with previous predictions by Murakami and Okada (2006), who used simulations of single neocortical PNs to predict that ∼50,000 synchronously firing neurons may be simultaneously active to produce an observable MEG response. Furthermore, Murakami et al. (2002, 2003) have shown in models of hippocampal neurons that back-propagation of PN action potentials can create measurable dipole currents directed away from the soma and that extracellular stimulation of apical dendrites can induce measurable currents in the opposite direction.
Weak exogenous drive reproduced the threshold-level response
By decreasing the strength of each of the synapses activated by the exogenous drive to the network, while keeping the same input sequence and timing (GR ∼25 ms/SGR ∼70 ms/late GR ∼135 ms), and all other network parameters fixed as in the simulation of the suprathreshold-level stimuli, the waveform of the observed dipole evoked from threshold-level stimulation emerged (Fig. 5 A, blue curve; for synaptic conductances, see Table 3, threshold level, nonperceived). With this simple manipulation, the differences between the simulated threshold and suprathreshold-level responses were almost entirely consistent with those seen experimentally. First, the weaker (50% smaller than the suprathreshold simulation) initial GR output at ∼25 ms was not strong enough to induce initial M25-M35-M50 peaks. This output did create a small positive peak before 50 ms because of back-propagation of current in the dendrites of both the L2/3 and L5 PN populations (Fig. 5 C). The subsequent negative deflection from repolarization of the dendrites could be seen in the L2/3 population (Fig. 5 C, gray curve); however, this effect was negligible in the total population response (Fig. 5 C, blue curve). This early peak was not distinguishable from the prestimulus level of activity in the experimental data (Figs. 3, 4), likely because of a lack of signal-to-noise in the MEG measurement. We anticipated that the initial peak could be eliminated in the model by increasing the level of baseline noise and/or decreasing the initial GR output strength. However, based on the known physiology of input to the cortex from the periphery, we assumed that the initial thalamic input/GR output was present at threshold, albeit weakly. Second, the magnitude of the M70 peak was smaller and, although there was no manipulation in the timing of input, the longer latency and smaller onset slope of the induced M70 peak emerged naturally from the underlying network dynamics, which exhibited a slower and smaller recruitment of net activity from the weaker (25% smaller than the suprathreshold simulation) SGR input (M70 magnitude, p < 0.001; ST, −543 ± 7 nA · m; T, 70 ± 29 nA · m; latency, p < 0.001; ST, 76 ± 5 ms; T, 79 ± 7 ms; SD, 7; onset slope, p < 0.001; ST, −22 ± 6 nA · m/ms; T, −2 ± 1 nA · m/ms). As in the suprathreshold simulation, the activity in the L5 PN population dominated the M70 peak because of the larger length and diameter of the L5 PN dendrites (Fig. 5 C, compare gray and black curves). The decrease in network activity also recreated the smaller magnitude and onset slope to the M100 peak (M100 magnitude, p < 0.003; ST, 453 ± 2 nA · m; T, 66 ± 67 nA · m; onset slope, p < 0.001; ST, 26 ± 9 nA · m/ms; T, 3 ± 2 nA · m/ms), which was again dominated by the L5 network (Fig. 5 C, black curve). Last, the weaker (7% smaller than the suprathreshold simulation) late GR output at ∼135 ms also induced a smaller and slower network response, recreating the smaller magnitude and increased latency of M135 peak (M135 magnitude, p < 0.001; ST, 860 ± 56 nA · m; T, 230 ± 68 nA · m; latency, p < 0.001; ST, 137 ± 7 ms; T, 142 ± 6 ms), and the mean difference in the area under the curve between the M100 and M135 peaks (p < 0.001; ST, 26 ± 8 nA · m; T, 8 ± 4 nA · m).
Timing and magnitude of SGR input and late GR output predict perception
Significant differences between the experimentally observed SI evoked responses from the perceived and nonperceived threshold-level stimuli emerged beginning with the M70 peak (Fig. 4). Our model results predict that the M70 peak and subsequent activity come from exogenous drive to the SGR of SI at ∼70 ms, followed by a later output from the GR at ∼135 ms (Fig. 5). As such, an additional prediction from the model was that observed differences in the SI evoked responses for perceived and nonperceived trials arose from differences in the timing and magnitude of these drives. Specifically, we predicted that on perceived trials the SGR input and late GR output arrive earlier and are stronger. By simulating these effects in our SI model, we reproduced the observed differences in the evoked response with perception.
Figure 6 shows the SI evoked responses for simulated perceived (dark blue) and nonperceived (light blue) trials averaged over 100 trials. The results in this figure are comparable with those in Figure 4 A, where several statistically significant differences between perceived and nonperceived trials emerged. The parameters reproducing the response from nonperceived trials were the same as those producing the threshold-level response in Figure 5 A. To reproduce the response from perceived trials, the mean latencies of the SGR input and late GR output were each decreased by 5 ms (to 65 and 130 ms, respectively), and the synaptic strengths were increased by 5 and 30% from the default values, respectively (see Table 3). These parameter changes reproduced the observed correlates of perception in SI, including an increase in M70 magnitude with perception (significant in the equal stimulus amplitude data in Fig. 4 B) (model statistics, p < 0.001; P, −70 ± 29 nA · m; NP, −93 ± 44 nA · m), the larger onset slope to the M100 peak (slope from M70 to M100, p < 0.003; P, 27 ± 34 nA · m/ms; NP, 16 ± 17 nA · m/ms), and the larger magnitude of the M100 and M135 peaks and mean area under the curve between them (M100 magnitude, p < 0.001; P, 112 ± 84 nA · m; NP, 65 ± 67 nA · m; M135 magnitude, p < 0.001; P, 253 ± 64 nA · m; NP, 230 ± 23 nA · m; area between 100 and 150 ms, p < 0.001; P, 13 ± 4 nA · m; NP: 8 ± 4 nA · m). Surprisingly, although the mean latency of the SGR input and late GR output were each decreased by 5 ms to simulate perceived trials, consistent with the experimental data, on the average the latency difference of the M70 remained significant, whereas the latency difference of the M135 peak emerged only as a trend (M70 latency, p < 0.001; P, 76 ± 7 ms; NP, 79 ± 7 nA · m; M135 latency, p = 0.1075; P, 136 ± 7 ms; NP, 137 ± 6 ms).
We combined human MEG and computational modeling to investigate the cortical dynamics underlying tactile detection. This combination of in vivo experimental data and biophysically constrained modeling led to four key results. First, brief tactile stimuli to a subject's fingertip evoked a consistent SI MEG signal that could be well localized to the position of area 3b. Second, the SI response predicted detection beginning ≥70 ms after stimulus. Third, simulation of an SI evoked response with a cortical model reproduced all major peaks recorded, providing a solution to the MEG signal based in cell-type- and lamina-specific activity patterns. This modeling led to the novel prediction that the polarity and magnitude of peaks in the SI response were induced by a sequence of exogenous excitatory drive consisting of GR output at ∼25 ms, followed by SGR input at ∼70 ms, followed by a late GR output at ∼135 ms. This sequence may be interpreted as initial input from the periphery through the lemniscal thalamus to GR, followed by feedback input from higher-order cortex or nonspecific thalamus to SGR, followed by a second wave of lemniscal thalamic input to GR. Using this network, we were also able to simulate differences in the magnitude and onset timing of peaks in the SI evoked response for suprathreshold- and threshold-level stimuli by decreasing the strength of exogenous drive. Fourth, specific manipulations of the model reproduced differences in observed MEG response that correlated with perception. Specifically, perceived threshold-level stimuli were characterized by SGR inputs and late GR outputs with earlier latencies and stronger amplitudes.
Because the synaptic architecture in the model was based on general principles of anatomy and physiology in a laminated cortex, our results are likely to be applicable to the interpretation of evoked MEG signals from primary sensory areas during a range of cognitive tasks, specifically in the interpretation of studies that report changes in peak latency and/or magnitude as a marker for changes in perception, attention, and brain disease.
Origin of evoked response and modulation with detection
Our model provides a realistic network that predicts the electrophysiological origin of the evoked SI response and its correlates with human perception. These predictions are consistent with previous studies. Studies focused specifically on investigating the neural origin of primary somatosensory evoked surface potentials (SEPs) (Towe, 1966; Arezzo et al., 1979, 1981; Allison et al., 1980; Gardner et al., 1984; Desmedt, 1988; Peterson et al., 1995) and fields (SEFs) (Hari and Forss, 1999; Ikeda et al., 2005) have primarily analyzed peaks emerging <50 ms after stimulus. A consensus from studies in which electrophysiological activity was recorded from multiple cortical layers is that the net laminar current flow that produces these early evoked responses derives from the postsynaptic current flow in PNs in the supragranular and infragranular layers (Towe, 1966; Peterson et al., 1995; Ikeda et al., 2005). Our data are consistent with a recent study by Ikeda et al. (2005) that showed with simultaneous MEG and intracranial laminar recordings in piglet SI that the first recorded SEF peak (N20, analogous to M25) from trigeminal nerve stimulation is generated by intracellular currents in two populations of excitatory neurons that appear to fire initially in the soma and produce back-propagating spikes toward distal apical dendrites. The predicted origins of the M35 and M50 peaks have yet to be verified.
A key hypothesis from our study is that tactile detection correlates with the SI evoked response beginning at the M70 peak. The M70 and M100 were reproduced in the model by simulating excitatory synaptic input at ∼70 ms to the SGR (layers I/II). The importance of this input is supported by the work of Kulics and Cauller, who recorded laminar profiles of local field potential and multiunit activity from the PoCG of awake monkeys during detection of cutaneous electrical stimulation (Kulics and Cauller, 1986; Cauller and Kulics, 1991; Jackson and Cauller, 1998). Their data showed an SEP between 50 and 65 ms, N1, corresponding to our M70 that was created by excitatory input to layers I/II that induced spiking in the deep layers. This result is consistent with our findings that excitatory input to the dendrites of the L2/3 and L5 neurons in the SGR induced spiking and bursting activity. Laminar profiles of SI activity that indicate excitatory input to the SGR at consistent latencies are also observed during tactile stimulation in anesthetized rats (Di et al., 1990; Barth and Di, 1991) and during brief tactile finger stimuli in awake monkeys (Lipton and Schroeder, personal communication).
There are two primary candidate sources of exogenous input to supragranular SI. One is “higher-order” cortical inputs, including possibly those from SII (Friedman et al., 1980; Cauller et al., 1998; Jackson and Cauller, 1998) or the frontal cortex (FC) and posterior parietal cortex (PPC) (Mauguiere et al., 1997a,b; Staines et al., 2002; Golmayo et al., 2003; Bauer et al., 2006). Human studies that show a peak in evoked activity in the parietal operculum at 70 ms during median nerve and tactile stimuli support SII as a source of M70 input (Karhu and Tesche, 1999; Hoechstetter et al., 2000, 2001). We investigated SII activation in the current dataset, but the signal-to-noise ratio from this second source modeled was too low to allow quantitative interpretation (data not shown), and no other sources could be estimated. Activations of the FC and PPC with somatosensory stimuli may require median nerve stimulation (Forss et al., 1994a,b, 1996; Mauguiere et al., 1997a,b). A second possible source of SGR input is a class of widespread calbindin-staining thalamic neurons that project diffusely to these layers (Jones, 2001).
To recreate the M135 activity in the SI model, a second output from the GR was necessary. This second wave of lemniscal thalamic input is consistent with the induction of a reverberatory thalamic-cortical-thalamic cycle (Kandel and Buzsáki, 1997; Jones, 2001; Guillery and Sherman, 2002; Nicolelis and Fanselow, 2002; Sherman and Guillery, 2002). Laminar profiles of activity that show patterns consistent with input to the GR at 100–140 ms latencies have been observed in SI of awake rats during thalamic stimulation (Kandel and Buzsáki, 1997) and in awake monkeys after tactile stimulation (Lipton and Schroeder, personal communication).
The correlates of detection in SI were reproduced in the model by decreasing the latency and increasing the strength of the later (∼70 and ∼135 ms) exogenous drives. The decrease in latency was modeled as a 5 ms decrease in the mean of the Gaussian distribution of presynaptic spike timing over trials. Thus, a prediction of our findings is a decrease in firing latency in the presynaptic areas that drive the M70/M100 and M135 responses. The increase in strength was modeled as an increase in the synaptic conductance of the targeted PN dendrites. Several mechanisms could underlie this increase, including an increase in the active conductance of the synapse, possibly via neuromodulators (McCormick, 1992; Hasselmo, 1995; Sarter et al., 2005; Yu and Dayan, 2005), or changes in the presynaptic source such as an increase in synchrony and/or firing rate. Support for such changes has been observed in SII during somatosensory detection and attention tasks (Garcia-Larrea et al., 1991; Hsiao et al., 1993; Hoechstetter et al., 2000; Steinmetz et al., 2000; Eimer and Forster, 2003; de Lafuente and Romo, 2006; Moore et al., 2007) and in FC during tactile detection (de Lafuente and Romo, 2005, 2006; Palva et al., 2005), and implicated in a class of thalamic neurons that project directly to SGR (Jones, 1998, 2001).
SI correlates of detection: comparison with previous studies
Other studies have investigated tactile detection in awake primates and have found correlations between SI activity and perception. A MEG study by Palva et al. (2005) used a fixed-amplitude electrical current with an immediate motor report. Although the signal was not localized with ECD methods, the rectified response (30–150 ms) from sensors over sensorimotor cortex was larger in magnitude for perceived stimuli. Presuming that the activity they recorded was generated in part in SI, these results predict a subset of our findings. In the present study, the ECD SSP method allowed localization of activation to area 3b, and computational neural modeling demonstrated a specific sequence of cellular events through which the SI circuit can produce these changes.
Kulics and Cauller (Kulics, 1982; Kulics and Cauller, 1986; Cauller and Kulics, 1991) observed a correlation between the magnitude and latency of the N1 and detection in awake monkeys. Given the similarity in the origin and character of the N1 and M70, our findings are in good agreement. Our finding that the magnitude of the M135 correlates with tactile detection is also in agreement with the correlation they observed between the magnitude of the later P2 peak (105–130 ms) and behavioral response latency.
Our findings are in apparent disagreement, however, with recent studies from macaque monkeys trained to tactile detection (de Lafuente and Romo, 2005, 2006). Our MEG signal and the activity in the model that provided a robust fit to the data predict robust increases in action potential firing rate on detected trials in SI, whereas these previous studies reported no difference in SI for hits and misses. This discrepancy could have many explanations, including differences in stimuli, subject training, and/or species differences. However, the strong agreement between the present findings and other data from well trained monkeys (Kulics, 1982; Kulics and Cauller, 1986; Cauller and Kulics, 1991) and rats (Krupa et al., 2004) argues against the latter interpretations. Another possible difference is the mode of electrophysiological recording conducted, the measurement of single neuron action potential activity versus synchronous activity across populations of neurons. That said, our modeling strongly suggests that differences in action potential firing rate should be present in large L5 PNs, putatively a population that is frequently recorded using thresholded extracellular recording techniques. Additional studies targeted to reconciling these data are required.
State-dependent regulation of detection
Recent tactile detection studies reported an impact of ongoing rhythmic activity in rolandic cortex on perception (Linkenkaer-Hansen et al., 2004; Palva et al., 2005), such that detection was dependent on power and phase locking in multiple frequency bands, with a dominant dependence in the α range (8–14 Hz) (see also Worden et al., 2000). Gamma frequency activity (30–80 Hz) has also typically been associated with attention and perceptual success (Meador et al., 2002; Gonzalez Andino et al., 2005; Bauer et al., 2006), as have nonrhythmic baseline magnitudes (Martin et al., 2006). Analysis of state properties is an essential and extensive topic. As such, it was beyond the scope of the present report, which focused on poststimulus response dynamics, accurate biophysical modeling, and their relation to conscious perception.
This work was supported by National Institutes of Health Grants P41RR14075, K25MH072941, 1R01-NS045130-01, and T32 GM007484; National Science Foundation Grant 0316933; the Athinoula A. Martinos Center for Biomedical Imaging; the McGovern Institute for Brain Research; and the MIND (Mental Illness and Neurodiagnostic Discovery) Institute. We thank Michael Sikora and Michael Hines for excellent technical support in implementation of network structure in NEURON software code. We also thank Charles Schroeder and Michael Lipton for valuable discussions of preliminary data from laminar recordings in awake monkeys, Seppo Ahlfors for insightful comments on MEG equivalent current dipole sources, and Dahlia Sharon and Talia Konkle for thorough remarks on this manuscript.
- Correspondence should be addressed to Stephanie R. Jones, Massachusetts General Hospital, Athinoula A. Martinos Center for Biomedical Imaging, 149 13th Street, Suite 2301, Charlestown, MA 02129.