QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

HOME  |   SEARCH  |   ARCHIVE  |   SUBSCRIBE  |   CONTACT  |   HELP

The Journal of Neuroscience, October 3, 2007, 27(40):10751-10764; doi:10.1523/JNEUROSCI.0482-07.2007

This Article Abstract Full Text (PDF) Submit an eLetter Alert me when this article is cited Alert me when eLetters are posted Alert me if a correction is posted Citation Map Services Email this article to a friend Related articles in J. Neurosci. Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Jones, S. R. Articles by Moore, C. I. Search for Related Content PubMed PubMed Citation Articles by Jones, S. R. Articles by Moore, C. I.

 Previous Article  |  Next Article 

Behavioral/Systems/Cognitive
Neural Correlates of Tactile Detection: A Combined Magnetoencephalography and Biophysically Based Computational Modeling Study

Stephanie R. Jones,¹ Dominique L. Pritchett,² Steven M. Stufflebeam,¹ Matti Hämäläinen,¹ and Christopher I. Moore^1,2

¹Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Charlestown, Massachusetts 02129, and ²McGovern Institute for Brain Research, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139

	Abstract

Top Abstract Introduction Materials and Methods Results Discussion References

 
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.

Key words: computational model; magnetoencephalography; dendritic processes; conscious perception; network dynamics; somatosensory cortex

	Introduction

Top Abstract Introduction Materials and Methods Results Discussion References

 
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

Top Abstract Introduction Materials and Methods Results Discussion References

 
MEG experiment
Subjects. Seven neurologically healthy, right-handed, 18- to 45-year-old adults were studied.

Stimulus paradigm. 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 x 7.8 x 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. 1A, 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.

View larger version (88K): [in this window] [in a new window]   Figure 1. Localization of SI activity. Estimated SI ECD localizations (blue dots) and orientations (blue lines) overlaid on the subjects' structural MRI brain images. The response evoked by stimulus to the third digit of the right hand was localized in the anterior bank of the contralateral postcentral gyrus, the position of area 3b. A, B, Example localizations and orientations from two subjects using the inverse solution technique described in Materials and Methods. C, Example localization for one of two subjects for which the SI ECD was placed (see Materials and Methods).

 
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. 1C) (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 2A. 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. 2C,D), on the ECD produced by the SI population.

View larger version (41K): [in this window] [in a new window]   Figure 2. Model SI network architecture. Ten PNs and 3 INs were included per layer. Excitatory (dark green) and inhibitory (red) synaptic connections were set as depicted. Bold outlined dendrites were contacted. A, Local synapses. Within-layer PN-to-PN synapses (not shown) were also present on dark green outlined dendrites. Each set of synaptic weights had a Gaussian spatial profile (Table 2). B, Spiking patterns evoked by somatic injected current (1 nA, 100 ms; no synaptic input). C, Connection pattern of output from the GR. The black arrow is only schematic, because lemniscal thalamic input was not explicitly modeled. D, Connection pattern of exogenous input to the SGR, presumably from a higher-order cortical and/or nonspecific thalamic neurons. The output from GR and input to SGR were modeled as spike train generators with a predetermined temporal profile and synaptic strength (Table 3).

 
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. 2A 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 · cm² for L5 and L2/3; C_m = 0.85 and C_m = 0.6195 µF/cm² 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).

View this table: [in this window] [in a new window]   Table 1. Dimensions of pyramidal neuron compartments in micrometers

 
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. 2B) 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/cm²) 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 x 10^–9. The ionic reversal potentials (in mV) were as follows: E_Na = –50; E_K = E_M = 77; E_L = –65.

Inhibitory interneurons. 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/cm²; g_Na = 0.12 S/cm²; g_K = 0.036 S/cm²; g_L = 0.003 S/cm²; 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 2A. 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 (GABA_A/GABA_B) 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; GABA_A, 0.5/5 ms, –80 mV; GABA_B, 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.

View this table: [in this window] [in a new window]   Table 2. Local network synaptic connection parameters

 
Exogenous drive. 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. 2C). 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. 2D). 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.

View this table: [in this window] [in a new window]   Table 3. Exogenous synaptic input parameters

 
Parameter fitting. 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.

Simulations. 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.

	Results

Top Abstract Introduction Materials and Methods Results Discussion References

 
MEG experiments
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 1C 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. 3A, left) and threshold (Fig. 3A, right)-level stimuli, averages across all subjects (Fig. 3B), and individual and mean differences between suprathreshold- and threshold-level responses (Fig. 3C). The average response to the suprathreshold stimulus (Fig. 3B, 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).

View larger version (30K): [in this window] [in a new window]   Figure 3. SI evoked responses. A, SI ECD responses from individual subjects (n = 7; baseline, mean 0–20 ms subtracted for display purposes) for suprathreshold-level (left) and threshold-level (right) stimuli. B, Average of suprathreshold-level (red curve) and threshold-level (blue curve) stimuli over all subjects. Consistent peaks emerged in the grand averages from 0 to 175 ms as labeled. Early peaks were not observed in the threshold response. C, Colored curves, Difference between suprathreshold- and threshold-level responses from individual subjects. Black curve, Mean difference over all subjects. Dark gray curve, Mean difference over subjects excluding subject with smallest difference (red curve). Light gray curve, Mean difference over subjects excluding subject with largest difference (cyan curve). The mean differences show that the stimulus amplitude differences were not driven by the response of an outlier subject.

 
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. 3B, 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 3C 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. 4A,B).

View larger version (35K): [in this window] [in a new window]   Figure 4. The magnitude and timing of the SI evoked response predict perception. A, Threshold-level stimulus SI evoked responses averaged over perceived (dark blue) and nonperceived (light blue) trials for varying stimulus amplitudes that were dynamically maintained at 50% perceptual threshold (trials = 100; n = 7). On perceived trials, the onset slope from the M70 to the M100 peak was larger, and the magnitudes of the M100 and M135 peaks and area under the curve between them were larger. B, Average threshold-level responses comparing perceived and nonperceived trials that had equal stimulus amplitudes (trials = 19 ± 9; n = 6; for details, see Results). On perceived trials, the magnitude of the M70 peak and onset slope from M70 to M100 were larger in this comparison. C, Average threshold-level responses sorted by stimulus amplitude. Dark green, Larger stimulus amplitudes; light green, smaller stimulus amplitudes (trials = 100 trials; n = 6). There were no statistically significant differences in this comparison. Error bars represent SEM.

 
In the first analysis, we averaged the last 100 trials of perceived and nonperceived stimuli to minimize within-session training effects (Fig. 4A) (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. 4B). 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. 4C), 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. 4C, 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 2A. 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. 2B). 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. 2A), and exogenous drive simulated as (1) output from activity in the GR, representative of activity induced by the lemniscal thalamus (Fig. 2C), and (2) input to the SGR, representative of input from higher-order cortex or nonspecific thalamic sources (Fig. 2D). 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 5A 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 5A (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.

View larger version (26K): [in this window] [in a new window]   Figure 5. Simulated SI evoked responses for suprathreshold- and threshold-level stimuli. A, The timing and location of synaptic inputs sets SI ECD polarity in the model. Red curve, Suprathreshold response. Output from activity in the GR at 25 ms reproduces the initial M25-M35-M50 peaks, exogenous SGR input at 70 ms created the subsequent M70 and M100 peaks, and second later GR output at 135 ms induced the M135 peak. Input times were selected over trials from Gaussian distributions displayed schematically in green with arrows marking the earliest input times (Table 3) (average of 100 trials). Blue curve, Threshold response. Decreasing the synaptic strengths by 50% (GR), 25% (SGR), and 7% (late GR) reproduced the waveform for the threshold-level stimulus (compare with Fig. 3B). B, Top, Contributions to the net current dipole during the suprathreshold response separated by layer. Bottom, Example of the activity from an individual L2/3 (left) and L5 (right) PN on a single trial; top traces, separate contributions from the basal (green) and apical (blue) dendritic compartments to the total current dipole (red) produced by the neuron; bottom traces, somatic membrane potential showing action potentials (black). C, Contributions to the net current dipole during the threshold response separated by layer. B, C, Schematics of network architecture drawn at each peak and arrows describe the direction of the net intracellular current flow within the pyramidal neurons that determines the polarity of the peak.

 
M25 origin
The separate contribution to the population average suprathreshold response from L2/3 (gray) and L5 (black) PN populations are shown in Figure 5B (top) along with representative contributions from a single L2/3 and L5 PN during a single trial (bottom). As shown in Figure 5B (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 5B (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. 5A, 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.

M35 origin
A subsequent repolarization of the L2/3 dendrites created a net downward current flow and induced the negative M35 peak (Fig. 5B, 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. 5B, 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.

M50 origin
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.

M70 origin
The large negative M70 peak (Fig. 5A, 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. 5B, 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. 5B, 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. 5B, 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.

M100 origin
After each spike in the burst, there was an active back-propagation of current up the apical dendrites of the L5 PNs (Fig. 5B, 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. 5A, red curve).

After the M100 peak activity, the simulated dipole current returned to baseline, creating an apparent discrepancy between the simulated evoked dipole (Fig. 5A) and the observed average evoked dipole (Fig. 3B), 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.

M135 origin
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. 5A, red curve).

The arrows located next to the PNs in the network schematics in Figure 5B (top; similarly in Fig. 5C) 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. 5A,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. 5A, 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. 5C). The subsequent negative deflection from repolarization of the dendrites could be seen in the L2/3 population (Fig. 5C, gray curve); however, this effect was negligible in the total population response (Fig. 5C, 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. 5C, 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. 5C, 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 4A, 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 5A. 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. 4B) (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).

View larger version (27K): [in this window] [in a new window]   Figure 6. Simulated SI evoked responses for perceived and nonperceived trials (compare with Fig. 4A). The statistically significant differences in the SI evoked responses for perceived (dark blue) versus nonperceived (light blue) trials were reproduced in the model by decreasing the mean latency and increasing the synaptic input strength of the SGR input and late GR output. For perceived trials, the mean SGR input and late GR output latencies were each decreased by 5 ms (to 65 and 130 ms, respectively) and their strengths were increased by 5% and 30%, respectively (Table 3) (average of 100 trials). The initial GR output was fixed. The green arrows and schematic Gaussians mark the distributions and earliest input times for perceived trials. These manipulations reproduced the increase in M70 magnitude, the larger onset slope to the M100 peak, the larger magnitude of the M100 and M135 peaks, and the greater mean area under the curve between them. Nonperceived trials were simulated with default parameters that produced the threshold-level response in Figure 5A.

 

	Discussion

Top Abstract Introduction Materials and Methods Results Discussion References

 
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