Variations in cortical oscillations in the alpha (7–14 Hz) and beta (15–29 Hz) range have been correlated with attention, working memory, and stimulus detection. The mu rhythm recorded with magnetoencephalography (MEG) is a prominent oscillation generated by Rolandic cortex containing alpha and beta bands. Despite its prominence, the neural mechanisms regulating mu are unknown. We characterized the ongoing MEG mu rhythm from a localized source in the finger representation of primary somatosensory (SI) cortex. Subjects showed variation in the relative expression of mu-alpha or mu-beta, which were nonoverlapping for roughly 50% of their respective durations on single trials. To delineate the origins of this rhythm, a biophysically principled computational neural model of SI was developed, with distinct laminae, inhibitory and excitatory neurons, and feedforward (FF, representative of lemniscal thalamic drive) and feedback (FB, representative of higher-order cortical drive or input from nonlemniscal thalamic nuclei) inputs defined by the laminar location of their postsynaptic effects. The mu-alpha component was accurately modeled by rhythmic FF input at approximately 10-Hz. The mu-beta component was accurately modeled by the addition of approximately 10-Hz FB input that was nearly synchronous with the FF input. The relative dominance of these two frequencies depended on the delay between FF and FB drives, their relative input strengths, and stochastic changes in these variables. The model also reproduced key features of the impact of high prestimulus mu power on peaks in SI-evoked activity. For stimuli presented during high mu power, the model predicted enhancement in an initial evoked peak and decreased subsequent deflections. In agreement, the MEG-evoked responses showed an enhanced initial peak and a trend to smaller subsequent peaks. These data provide new information on the dynamics of the mu rhythm in humans and the model provides a novel mechanistic interpretation of this rhythm and its functional significance.
Two predominant rhythms are expressed in the neocortex in the frequency range from 7 to 30 Hz: alpha (7–14 Hz) and beta (15–29 Hz). Modulation of alpha and beta activity is correlated with successful perception in humans and awake monkeys (Bauer et al. 2006; Donner et al. 2007; Hanslmayr et al. 2007; Linkenkaer-Hansen et al. 2004; Mathewson et al. 2009; Mazaheri et al. 2009; Palva et al. 2005b; Pineda 2005; Schroeder and Lakatos 2009a; Schubert et al. 2008; van Wijk et al. 2009; Wilke et al. 2006; Worden et al. 2000; Zhang and Ding 2009). Recent studies have emphasized a potential role for the active deployment of these rhythms in the suppression of “distracting” sensory input (Jensen et al. 2002; Kelly et al. 2006; Mazaheri et al. 2009; Worden et al. 2000), presumably by suppression of evoked responses in early sensory cortices.
The mu rhythm measured with magnetoencephalography (MEG) over Rolandic cortex shows alpha and beta components (Hari and Salmelin 1997; Tiihonen et al. 1989). This finding is in contrast to the Rolandic mu rhythm measured with electroencephalography (EEG), in which only a dominant alpha component is typically observed (Kuhlman 1978; Zhang and Ding 2009). This historical distinction is likely attributable to differences in the recording techniques and has led to mixed usage of the term “mu” in the literature. This ambiguity in naming is indicative of the ongoing ambiguity with respect to the statistical characteristics and neural origins of the mu rhythm. Despite the fact that much research has been devoted to localizing the source of this rhythm in the brain—and to understanding the cellular-level neural mechanisms creating alpha and beta rhythms independently—the neural origin of the MEG mu complex remains unknown. In the present report, we investigated the two-component mu rhythm measured with MEG using experimental and modeling approaches. We refer to these components throughout as mu-alpha and mu-beta.
One prominent view of the origin of the MEG mu rhythm, based on source localization of sensor data from human studies, is that the mu-beta component is produced by the precentral motor cortex, whereas the mu-alpha component originates from the postcentral somatosensory cortex (Hari and Salmelin 1997; Salmelin and Hari 1994; Salmelin et al. 1995). These studies focused on localizing late event-related desynchronization (ERD) of the rhythm after movement. More recent work, focused on spontaneous activity and early ERD, has shown that both components can be expressed in a single area (Brovelli et al. 2004; Gaetz and Cheyne 2006; Kopell et al. 2000; Pinto et al. 2003; Szurhaj et al. 2003), with intracerebral recordings in humans suggesting a common source in primary somatosensory (SI) cortex (Szurhaj et al. 2003).
Studies that have investigated the cellular-level neural mechanisms inducing ongoing cortical alpha and beta rhythms have focused on the origin of the two frequency bands separately. A large body of experimental and computational work suggests neocortical alpha emerges from an approximately 10-Hz thalamocortical rhythm (Andersen and Andersson 1968; Contreras and Steriade 1995; Hughes and Crunelli 2005; Suffczynski et al. 2001; Traub et al. 2005). Other evidence suggests that neocortical alpha also depends on, or could emerge independently from, intrinsic properties in large layer V pyramidal neurons (Bollimunta et al. 2008; Jones et al. 2000; Pinto et al. 2003; Silva et al. 1991) and/or the local activity of low-threshold spiking interneurons (Fanselow et al. 2008).
Fewer studies of the neural origins of the cortical beta rhythm have been conducted. Slice recordings (Roopun et al. 2006; Whittington et al. 2000) and computational models (Jensen et al. 2002; Kopell et al. 2000; Pinto et al. 2003; Roopun et al. 2006) have shown a beta-frequency range oscillation in isolated cortex that depends on the kinetics of the M-type potassium current in excitatory neurons, combined with GABAergic inhibition. Roopun et al. (2006) further found that axonal gap junctions were critical to the maintenance of pharmacologically induced 20- to 30-Hz rhythms in slices from somatosensory cortex (referred to as a “beta 2” rhythm) and that lower-frequency beta rhythms (13–17 Hz; “beta 1”) could be produced by period concatenation of higher-frequency beta (20–30 Hz) and gamma (30–50 Hz) oscillations (Kramer et al. 2008; Roopun et al. 2008). In all of these studies, the beta rhythm is argued to emerge from the activity of localized cells in a single cortical circuit that are spiking nearly synchronously and that the frequency of these rhythms depends on intrinsic membrane time constants.
Another body of research has suggested that beta activity mediates long-range communication between cortical areas, suggesting that intracortical projections may play a role in local beta expression (Buschman and Miller 2007; Hanslmayr et al. 2007; Roelfsema et al. 1997; Schubert et al. 2008; von Stein et al. 2000; Witham et al. 2007). For example, Von Stein et al. (2000) found coherence in beta-band activity between temporal and parietal EEG sensors during multimodel object representation in humans, and Buschman and Miller (2007) between parietal and frontal indwelling electrodes during selective attention in monkeys. In support for a role of interareal projections in the neocortex, Whitham et al. (2007) found interactions between 20-Hz activity in somatosensory and motor areas that did not depend on intrinsic spiking (Witham and Baker 2007) and claimed there was oscillatory coupling across the central sulcus. Brovelli et al. (2004) applied Granger causality analysis to 20-Hz oscillatory activity measured intracranially from multiple cortical sites in the monkey and claimed that the beta activity propagated from the primary somatosensory (SI) cortex to primary motor and parietal cortices. Nonlemniscal thalamic nuclei have also been proposed as relays for signal transmission between neocortical areas (Sherman 2005). Such thalamic nuclei are also ideally poised to generate oscillatory coherence across multiple cortical areas (Llinás and Ribary 2001). As such, beta activity related to long-range communication could also arise, at least in part, from these projections.
Difficulty establishing the neural mechanisms generating the MEG mu rhythm comes in part from the fact that the shared temporal dynamics of the mu-alpha and mu-beta components are not well characterized. Tiihonen et al. (1989) observed qualitatively that the mu-alpha and mu-beta components do not always overlap in time, suggesting they are not harmonics, but did not quantify this assertion (Palva et al. 2005b; Tiihonen et al. 1989). Evidence for cross-frequency coupling of these components has also been reported (Palva et al. 2005b). Quantifying the degree to which mu-alpha and mu-beta co-occur and covary in amplitude is crucial to understanding their relative interdependence and to constraining computational models of their origins.
An important related topic is delineation of the mechanistic underpinnings of the impact of ongoing mu power on evoked sensory responses. This relation is particularly important, given studies showing that both ongoing mu power and SI tactile-evoked responses predict perception (Jones et al. 2007; Kulics 1982; Linkenkaer-Hansen et al. 2004; Zhang and Ding 2009). In recent studies of mu-alpha, prestimulus power shows an inverted U-shaped relation to tactile detection (Linkenkaer-Hansen et al. 2004; Zhang and Ding 2009) and a similar relationship was seen for mu-beta (Linkenkaer-Hansen et al. 2004). In studies of the relationship between time-domain peaks in tactile-evoked responses and detection, Kulics (1982) observed that greater evoked activity in the local field potential (LFP) near 70 ms predicted detection of tactile input by monkeys with electrodes implanted in SI and that differences in later activity (105–130 ms) were correlated with reaction time. Similarly, in humans, Jones et al. (2007) showed that greater amplitude of components of the SI-evoked response 70–130 ms poststimulus was a key predictor of detection. Zhang and Ding (2009) recorded greater late evoked components (∼140 ms) in sensorimotor EEG electrodes on successfully detected trials and Palva et al. (2005) showed greater response 30–150 ms poststimulus for detection in sensorimotor MEG sensors.
Despite the links between prestimulus mu expression or evoked response amplitude and detection, the impact of spontaneous mu on peaks in the time-domain tactile-evoked response (ER) has received limited attention. Nikouline et al. (2000) reported that large variation in the prestimulus mu-alpha band activity was related to “relatively stable” early parts (≤60 ms; referred to as P35 and P60) of median nerve ERs recorded using MEG, which showed a small positive correlation with greater alpha power predicting a larger early component (Nikouline et al. 2000). Zhang and Ding (2009) reported that for late components (∼140 ms) of the tactile ER measured with EEG, there was a complex and largely parabolic relationship between mu-alpha power and evoked amplitude (Zhang and Ding 2009). Further, although mechanisms for changes in later (>250 ms) components of median nerve (Nikulin et al. 2007) and visual (Mazaheri and Jensen 2008) evoked responses have been associated with baseline shifts and phase resetting of alpha activity (Hanslmayr et al. 2007; Makeig et al. 2002), to our knowledge no systematic study has been performed to delineate the specific impact of MEG mu power (containing alpha and beta components) on different aspects of the earlier tactile responses that have been linked to perceptual success.
To investigate these questions, we recorded whole-head 306-channel MEG data and applied an equivalent current dipole inverse solution technique to look at activity from the hand area of SI, a technique that has been found to consistently localize signals discretely to the anterior bank of the postcentral gyrus, area 3b (Jones et al. 2007). We studied the prestimulus mu rhythm generated by this somatosensory source and its connection to peaks in the early (<175 ms) time-domain–evoked responses. On single trials, we observed that during the prestimulus time period high-power epochs of mu-alpha and mu-beta were not simultaneous in their emergence, but did co-occur at rates greater than chance. These findings suggest that mu-alpha and mu-beta emerge from separable generators within SI that share common elements. We did not explore evoked poststimulus oscillations because our analysis of this period was focused on the region in which our model had predictive power in the evoked response (0–175 ms).
We then implemented a biophysically principled laminar cortical model of SI to make specific predictions with respect to the underlying neural mechanisms inducing the mu rhythm, its modulations in prestimulus power, and their influence on early peaks (<175 ms) in the time-domain tactile-evoked response. This model was expanded from our previous SI model that accurately reproduced the different time-domain–evoked response components in SI measured with MEG. The model included pyramidal neurons and interneurons in the supra- and infragranular layers. Importantly, this model also included thalamocortical feedforward (FF) input and inputs to distal dendrites that could arise from intracortical feedback (FB) or from nonlemniscal thalamic input.
Our results showed that the mutual generation of the mu-alpha and mu-beta components of the prestimulus rhythm, as well as the relative separation of these components in time, can be reproduced by the SI model. The data and model indicated these rhythms are not simple harmonics and are not regulated by time constants of intrinsic membrane properties. The mu-alpha component was reproduced with an approximately 10-Hz stochastic thalamocortical FF input. Accurately capturing mu-beta and its interactions with mu-alpha, in contrast to previous modeling studies emphasizing local cortical circuits and intrinsic properties, required rhythmic supragranular input, consistent with FB from higher cortical areas and/or input from nonlemniscal thalamic sources. To accurately reproduce the observed mu-alpha and mu-beta emergence and statistical interdependence, this FB was simulated not as a 20-Hz signal, but rather as a stochastic approximately 10-Hz signal. Further, mu-beta was not achieved when the approximately 10-Hz supragranular inputs were perfectly out of phase, as one might expect for generation of a 20-Hz oscillation, but rather almost perfectly aligned with the approximately 10-Hz FF input (<10-ms mean delay).
Investigation of the influence of prestimulus mu rhythm on the gain of the early-evoked sensory response (0–175 ms) in the model showed that mu power had a dominant effect on the initial peaks of the time-domain signal (<70 ms). This prediction was confirmed in the MEG data. Further, the model predicted that a key impact of high mu on peaks of the early evoked response is enhanced recruitment of excitatory and inhibitory interneurons and that recruitment of the interneurons caused suppression of subsequent (70–100 ms) components.
MEG data were collected from 10 neurologically healthy, right-handed, 18- to 45-yr-old adults during performance of a tactile-detection paradigm. The experimental protocol was approved by the Massachusetts General Hospital Internal Review Board and each subject gave informed consent prior to data acquisition. The stimulus paradigm and data acquisition are described in detail in Jones et al. (2007), which reported on tactile-evoked response components for 7 of these subjects; data from an additional 3 subjects were collected for the present report. Here, we outline key aspects of the experimental paradigm and current data analysis methods.
Brief taps were delivered to the subject's right hand in the form of a single cycle of a 100-Hz sine wave (10-ms duration) via a custom piezoelectric device. Subjects rested their hand on a Delrin frame that held a piezoelectric transducer parallel to the finger (Noliac ceramic multilayer bender plate: 32 × 7.8 × 1.88 mm). A deflection stroke drove a 7-mm-diameter Delrin contractor through a 1-cm circular rigid surround and into the fingertip.
The detection threshold of each subject was obtained prior to imaging using a parameter estimation by sequential testing (PEST) convergence procedure (Dai 1995; Leek 2001). During MEG recordings, for 70% of presented trials, stimulus strength was maintained at a perceptual threshold level (50% detection) using a dynamic algorithm (see Jones et al. 2007). Suprathreshold stimuli (10% of all trials; 350-μm deflection; 100% detection) and null trials (20%) were randomly interleaved with the threshold stimuli. 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 2 s. During the auditory cue, the 10-ms finger-tap stimulus was delivered between 500 and 1,500 ms, in 100-ms intervals, from trial onset. The number of trials of a given latency to tap was randomly distributed during each run. Following the cessation of the auditory cue, subjects reported detection or nondetection of the stimulus with button presses, using the second or third digit of the left hand, respectively. The auditory cue ended ≥500 ms after the tactile stimulus and 1,000 ms before the next trial began.
MEG DATA ACQUISITION.
By use of a 306-channel MEG (Elekta-Neuromag Vectorview), 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) signals were recorded with electrodes placed close to the left eye. Four head-position indicator (HPI) coils were placed on the subject's head to coregister the subject's anatomical magnetic resonance image (MRI) and the MEG sensors. The data were sampled at 600 Hz with the band-pass 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 band-pass of 0.1–200 Hz. The chosen high-pass filter corner frequency was low enough to retain possible slow variations in the dc level of the neural signals while eliminating low-frequency environmental noise. Inspection of the localized SI spontaneous activity (see description of source analysis in the following text) showed that it was stable through ≥10-s time windows during which there was modulation of the mu rhythm without dc offset. Epochs with EOG peak-to-peak amplitude exceeding 150 μV were excluded from the analysis.
MEG SOURCE ANALYSIS.
Source analysis was used to locate the primary current-dipole source to contralateral SI and to find the time course of this source, taking into account the presence of other active areas. This method was motivated by and described in greater detail in Jones et al. (2007), but is also described here. The SI contribution to the somatosensory-evoked field was isolated using the following approach. Because we did not observe consistent activity over the ipsilateral secondary somatosensory cortex (SII), or in other brain areas, we modeled the data with two dipoles (contralateral SI and SII). This fit was then optimized by use of signal-space projection (SSP) (Tesche et al. 1995; Uusitalo and Ilmoniemi 1997). A least-squares fit with the dipole forward solution was calculated through the use of a spherically symmetric conductor model (Hamalainen and Sarvas 1989; Sarvas 1987). At the peak activity in the suprathreshold stimulus signals from one data run (average of 12 trials; mean = 68 ms, SD = 8 ms), we observed an initial equivalent current dipole (ECD). An anatomical MRI was then coregistered with the MEG, confirming this source localized to SI in 7 of 10 subjects. We then removed the contribution of the SI ECD using the SSP method and the residual data fit with a second ECD (Nishitani and Hari 2000; Tesche et al. 1995; Uusitalo and Ilmoniemi 1997). In all fit data during peak responses, the goodness-of-fit of the two-dipole model was >70%. We then removed the effect of the second ECD from the data using SSP and refitted the SI ECD to the residual. The ECD localizations fitted to the suprathreshold data were used to model all responses. Only the last 100 trials for a given response were considered for analysis, to account for adaptation and learning effects with training.
Large baseline rhythmic activity interfered with the localization of peak responses for 2 of the 10 subjects. In these cases, the SI source was placed in the finger representation of area 3b within the anterior bank of the postcentral gyrus (Moore et al. 2000; Penfield and Rasmussen 1950; Sastre-Janer et al. 1998; Uematsu et al. 1992; White et al. 1997; Yousry et al. 1997) and the SII dipole source was placed in the parietal operculum. For one of the 10 subjects, anatomical MRI data could not be obtained. In this case, SI dipole localization was determined by field contours on the spherical head model, consistent with the predicted position of contralateral SI. Removal of these 3 subjects from our analysis did not have a significant impact on our results (see Fig. 3).
FREQUENCY DOMAIN ANALYSIS.
Time–frequency representations (TFRs) or spectrograms of the data were calculated from 1 to 40 Hz on the SI ECD time courses by convolving the signals with a complex Morlet wavelet of the form w(t, f0) = A exp(−t2/2σt2) exp(2iπf0t), for each frequency of interest f0, where σt = m/2πf0 and i is the imaginary unit. The normalization factor was and the constant m, defining the compromise between time and frequency resolution, was 7. TFRs of power were calculated as the squared magnitude of the complex wavelet-transformed data. The normalization factor used is such that the sum of the magnitude of the wavelet coefficients over all frequencies is one, unlike that preserving the sum of squared magnitudes of the wavelet coefficients used in, e.g., Tallon-Baudry et al. (1997). Our normalization factor emphasizes higher frequencies with the distinct benefit of allowing us to more clearly visualize the time course of 15- to 29-Hz mu-beta activity (e.g., compare Fig. 3, A and C). Therefore the low-frequency activity (<7 Hz) is much less pronounced than that in the traditional PSD plots.
POWER SORTING ALGORITHM.
When sorting trials over mu (7–29 Hz), mu-alpha (7–14 Hz), or mu-beta (15–29 Hz) power for analysis in Figs. 5, 7, and 9B, spectrograms (described earlier) calculated for each trial were averaged over the frequency band of interest and then sorted from high to low power. The top and bottom 10% of the sorted trials were excluded from further analysis.
CALCULATION OF SYMMETRY INDEX.
The symmetry index of the oscillation waveform around zero, as shown in Fig. 6, was calculated as follows. We first applied a band-pass of the entire mu frequency range (7–29 Hz) to the signal. We then found the local maxima (peaks) and minima (troughs) in the filtered data and identified the time points corresponding to each peak and trough. These time points were used to obtain signal peak and trough values from the unfiltered data. The symmetry index (SInd) was calculated as [abs (peak) − abs (trough)]/[abs (peak) + abs (trough)] (Galaburda et al. 1990). A positive symmetry index indicates greater-amplitude peaks, a negative value indicates greater-amplitude troughs, and a zero value indicates that an oscillation was symmetric around zero.
Computational neural model
The computational neural modeling presented is expanded from our previous model of a laminar SI network, as described in detail in Jones et al. (2007) and the code is available to the public on the NEURON ModelDB website (http://senselab.med.yale.edu/modeldb/ShowModel.asp?model=113732). Here, we describe key features of the model and its expansion.
SI CORTICAL COLUMN MODEL.
Our simulated SI cortical column network consisted of 100 multiple-compartment excitatory pyramidal neurons (PNs) and 35 single-compartment inhibitory interneurons (INs) per layer (Thomson et al. 2002) in layers II/III and V, expanded from 10 PNs and 3 PNs in each layer in our previous model (Jones et al. 2007). Postsynaptic dendritic contact points of the local excitatory and inhibitory synapses are depicted in Fig. 1 A (see Jones et al. 2007 for supporting literature). Connection lines are schematic representations of axonal-to-dendritic input. Axons were not explicitly modeled. The PNs were arranged in a two-dimensional (2D) grid as shown in Fig. 1D. INs were interleaved evenly between every 2 PNs (not shown in Fig. 1D). Fast and slow excitatory (α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid/N-methyl-d-aspartate [AMPA/NMDA]) and inhibitory (γ-aminobutyric acid type A/type B [GABAA/GABAB]) synapses were simulated using an alpha function that was “turned on” by the soma of the presynaptic cell crossing a voltage threshold (mV = 0). 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 conductance of the synaptic connections within the local network grid were defined with a symmetric 2D 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 (WSCs, number of cells from center) are listed in Table 1 along with the minimum synaptic delay and corresponding Gaussian delay space constant (DSC).
SINGLE-CELL MORPHOLOGY AND PHYSIOLOGY.
The morphology and physiology of the INs in each layer were simulated with single compartments and contained fast sodium (INa) and potassium currents (IKdr) to create spiking activity. The PNs in layers II/III and V were simulated with eight and nine segments, respectively, based on the reduction by Bush and Sejnowski (1993), and used a compartment length of 50 microns. The PNs in layer II/III produced an adapting spike trains to injected current created by active currents in the somatic and dendritic compartments, including a fast sodium current (INa), a delayed rectifier potassium current (IKdr), an adapting potassium current (IM), and a leak current (IL). The layer V PNs produced bursting responses to injected current (Fig. 1B) and contained the same currents as those of the PNs in layer II/III, with the addition of a calcium current (ICa), a potassium-activated calcium current (IKCa), and h- and T-currents in the somatic and dendritic compartments. See Fig. 1B in Jones et al. (2007) for examples of spiking behavior. A low-threshold calcium current (IT) and a hyperpolarization-activated mixed cation current (Ih) are additions to our previous model. The IT channel density was constant in all segments at 2 × 10−4 S/cm2. The Ih channel density was increased exponentially from the soma (where it was set at 1 × 10−6 S/cm2) to apical dendrite, with a space constant of 3 × 10−3. This matches Ih densities observed in rat somatosensory cortex (Kole et al. 2006). For further details, supporting literature, and specific parameters see Jones et al. (2007) and the web-available code at http://senselab.med.yale.edu/modeldb/ShowModel.asp?model=113732.
SYNAPTIC ARCHITECTURE OF EXOGENOUS DRIVE TO SI.
Exogenous drives to the local network were simulated to reproduce 1) ongoing mu rhythms and 2) evoked responses in the model. In each case, the exogenous drive 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 (Douglas and Martin 2004; Felleman and Van Essen 1991; Friedman and Jones 1980; Jones 2001; Rockland and Pandya 1979). The sources of SI drive were modeled as spike generators with predefined temporal profiles and postsynaptic conductances that were distinct for the activity of ongoing rhythmic and evoked responses (described in the following text).
Feedforward (FF) drive emerged from the granular layer, layer IV, and contacted the supragranular L2/3 neurons, with a delayed and weaker connection to the infragranular layer V neurons (see Fig. 1B for specific poststimulus dendritic compartments). Activity in layer IV is modeled to reflect drive from the thalamus based on several studies of intracranial laminar electrophysiological recordings of evoked responses in SI, including responses to vibrissa and thalamic stimuli in rodents (Barth and Di 1991; Castro-Alamancos and Connors 1996; Di et al. 1990; Douglas and Martin 2004; Kandel and Buzsáki 1997), trigeminal stimulation in piglets (Ikeda et al. 2005), and to tactile (Cauller and Kulics 1991; Kulics and Cauller 1986) and median nerve stimuli in awake monkeys (Lipton et al. 2006; Peterson et al. 1995). The maximal conductances onto INs were always twice as strong as those onto PNs for all FF inputs (Giove et al. 2003). In our previous study, this FF drive was referred to as “granular layer output.” For descriptive purposes, we have changed the nomenclature in the current study.
Feedback (FB) drive to the SI network contacted the distal apical dendrites in the supragranular layers of each neural population (Fig. 1C). This connection was representative of input from higher-order cortical areas or nonspecific thalamic sources (Douglas and Martin 2004; Felleman and Van Essen 1991; Friedman et al. 1980; Jackson and Cauller 1998; Jones 2001; Rockland and Pandya 1979). In our previous study, this FB drive was referred to as “supragranular layer input.”
Temporal dynamics and conductances of ongoing and evoked exogenous drive to SI
Stochastic ongoing rhythmic drive (Fig. 2) was generated by delivering 10 “burst” spike trains, each consisting of two spikes with an interstimulus interval (ISI) of 10 ms, which was set based on experimental evidence (Hughes and Crunelli 2005), to the SI network in an FF synaptic activation pattern (Fig. 1C). To reproduce approximately 10-Hz input, all 10 bursts arrived nominally every 100 ms, with a Gaussian random distribution in arrival time for each burst (mean ISI of 100 ms; default variance of 400 ms). On every cycle, a similar input pattern (same number of bursts and arrival time statistics) arrived to the SI network in an FB synaptic activation pattern (Fig. 1D), delayed from the FF input by a fixed amount. This driving sequence is shown schematically in Fig. 2. The mean delay between the approximately 10-Hz FF and FB ongoing inputs, the number of prestimulus “bursts”/spike trains on each cycle (i.e., the “amplitude” of the input), the SD of the arrival time of each burst to the SI network (i.e., the “variance” of the input), and the maximal postsynaptic conductance of the inputs were varied parametrically to investigate their separate influence on an expressed mu rhythm (Fig. 8).
Stochastic events in each simulation were regulated by the arrival time of each of the input bursts, for both FF and FB inputs. On each cycle of the rhythmic input, the timing of each two-spike burst event (Fig. 2) was chosen from a Gaussian distribution, with mean fixed at 100-ms intervals (e.g., 0, 100, 200 ms, etc.) and variance of 400 ms. With this mechanism, on every cycle the relative net postsynaptic conductance of the FF and FB drives changed. All other parameters were fixed for each simulation with default values of delay = 5 ms, number of input bursts = 10, variance of input bursts = 400 ms. The default conductances for FF and FB inputs were as follows: low mu, maximal weight of 0.4 picosiemens (pS) onto PNs and 0.8 pS onto INs; high mu, maximal weight of 0.6 pS onto PNs and 1.2 pS onto INs. Synaptic weights and delays were distributed by a symmetric 2D Gaussian spatial profile with maximal weight and minimal delay in the center and WSC of 100 and DSC of 100 for all connections. The minimal delay between the FF input and the layer II/III PNs and INs was 0 ms and the layer V PNs and INs was 1 ms. The driving parameters for the ongoing rhythmic input in our model were all chosen so that the oscillations in the pyramidal neurons remained subthreshold (Zhu et al. 2009).
EVOKED RESPONSE INPUTS.
These were simulated as in Jones et al. (2007), with a sequence of FF input, followed by FB, followed by a reemergent late feedforward (LFF) input (schematically drawn in Fig. 9A). The timing of the input sequence was fixed as in Table 2. Each driving spike train consisted of a single presynaptic spike on each trial and the weights were distributed uniformly across the SI networks as in Table 2. The stochastic ongoing drive was unchanged during the evoked response sequence. Evoked responses were simulated to begin at various phases in simulated mu-alpha and mu-beta cycles as described in results.
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.
All simulations were performed using the shareware software program NEURON (available http://www.neuron.yale.edu/neuron/). A fixed time-step implicit Euler integration method was used with a time increment dt = 0.025 ms. Frequency analysis of all simulated rhythms was identical to that performed on the MEG data. Simulated evoked responses were smoothed by convolution with a 15-ms box filter. On publication, the code that produced all simulated data herein will be available on the ModelDB website http://senselab.med.yale.edu/senselab/modeldb/.
Somatosensory mu rhythm in trial averages.
Consistent with previous studies, the ongoing mu rhythm recorded here expressed peaks of activity in the mu-alpha (7–14 Hz) and mu-beta (15–29 Hz) frequency bands. The dual peak in mu can be observed in the mean and SE across subjects of the PSD calculated from the frequency range 1–60 Hz (Fig. 3A, blue curve; n = 10 Ss, mean of 200 1-s prestimulus trials per subject, calculated with Welch's periodogram; see methods). Peaks in both the mu-alpha and mu-beta range are evident in this analysis and the relative mu-alpha and mu-beta expressions vary across subjects (Fig. 3B). The gray curve in Fig. 3A shows the corresponding results when the data from the three subjects—whose SI dipole localizations were not determined by standard methods—were removed (see methods) and their individual subject data are shown in black in Fig. 3B. Although the overall power of the oscillations is smaller without these subjects, the overall shape is preserved and the dual-peak mu phenomenon is still present.
The two distinct components, mu-alpha and mu-beta, can be visualized more clearly when looking at time–frequency representations (TFRs; spectrograms) of the data from two example subjects (Fig. 3C, average n = 100 trials of 1-s prestimulus data, each calculated with Mortlet wavelets; see methods). The spectrograms from these subjects show clear bands of mu-alpha and mu-beta activity in the average data. We subsequently show that, although the relative mu-alpha and mu-beta expressions across subjects are not the same, the existence of separable activity in both frequency ranges persists in all subjects and in single trials.
Somatosensory mu rhythm on single trials.
These initial analyses suggested that there was not a fixed relationship between mu-alpha and mu-beta expressions across subjects, indicating different neural generators. To further investigate the relation between these bands, we analyzed the spontaneous SI signal in individual trials. Figure 4 shows two examples of SI frequency distributions and corresponding waveforms during 1-s prestimulus epochs for four subjects. These spectrograms are different from the average responses (Fig. 3C). In single trials, periods of high mu-alpha and mu-beta power often have nonoverlapping time courses. We note that subjects 3 and 4 had nonstandard dipole localizations and larger overall power (Fig. 3B). However, these subjects still showed the characteristic periods of nonoverlapping high mu-alpha and mu-beta power. This feature implies that these rhythms are not simple harmonics, further indicating separable neural generators.
To quantify this observation, we calculated the probability that high power in mu-alpha and mu-beta was expressed simultaneously in the spontaneous mu rhythm. We sorted the relative power of mu-alpha and mu-beta in each 100-ms window. We then calculated the probability that bins in the top third of the mu-alpha and mu-beta power distributions were identical. If these rhythms were generated from an identical neural process, they would be predicted to correlate perfectly. A completely random association would predict overlap on roughly 10% of bins. We found overlap on 50% of bins (Fig. 5, inset, 50.31% mean, SD = 0.056% n = 10 Ss, calculated from 2,000 100-ms time windows per subject). This finding implies a degree of independence for the two rhythms, but also that they co-occurred at greater than chance levels, suggesting they share components of neural mechanism.
To further characterize the relation between alpha and beta occurrence, we plotted a histogram of the ratio of alpha to beta power calculated in the same 100-ms time windows (Fig. 5, n = 2,000 trials, 10 Ss, 1,000 bins). The mean and median of the histogram were 2.4 and 1.3, respectively. The normalization applied to the TFR data in our analysis favors the observation of higher-frequency bands (see methods). As such, the dominance of alpha power in this analysis, despite the normalization applied, underscores the relative prevalence of this oscillation in the mu signal.
Computational neural model
Simulating MEG SI activity.
We have previously developed a laminar SI model that predicted the neural origin of the tactile-evoked response in MEG and provided insight into the changes in this response that predicted successful detection of a threshold-level stimulus (Jones et al. 2007). Here, we expanded this model from 10 pyramidal neurons (PNs) and 3 interneurons (INs) per layer, to a grid of 100 PNs and 35 INs per layer (see methods, Fig. 1) and used it to investigate the neural origin of the spontaneous SI mu rhythm, its modulations in power, and how these modulations influence cortical excitability and evoked response gain.
As in our previous study (Jones et al. 2007), the model consisted of a biophysically realistic laminar cortical model of a local SI network (Fig. 1). Two types of exogenous drive provided excitatory synaptic input to the SI column, one in an FF connection pattern representing input from the thalamus to layer IV and subsequently to layer II/III (Fig. 1B), and one in an FB connection emulating higher-order cortical sources to the superficial layer II/III (Fig. 1C).
The MEG SI ECD is calculated by the net longitudinal current flow within PNs in the model (Hamalainen et al. 1993; Hari and Salmelin 1997; Hari et al. 1980; Ikeda et al. 2005; Jones et al. 2007; Okada et al. 1997). Therefore the polarity of the waveform of a given oscillation provides specific insight into the direction of current flow within the pyramidal neurons across the cortical lamina. In a previous study (Jones et al. 2007), we determined that positive polarity corresponded to net current flow anterior, which in the case of area 3b corresponds to flow “up” the PN dendrites toward the cortical surface and negative polarity corresponded to net current flow down the PN dendrites toward the cell bodies/white matter (see methods). The example single-trial waveforms of the spontaneous SI mu rhythm in Fig. 4 (bottom panels) oscillated around zero for each subject, suggesting current flow that oscillates up and down the cortical laminae.
Nearly synchronous stochastic FF input followed by FB input creates a physiologically realistic mu rhythm.
Our strategy for simulating single-trial mu rhythms in the model was based on current theories as to the origin of 10- and 20-Hz activity in the cortex (see introduction and methods). In expanding our cortical model to study the mu rhythm, we included several intrinsic currents in our PNs that have previously been proposed to be important in the generation of alpha and beta activity (h-currents, T-currents, M-currents, Ca-currents). We initially thought that these currents would be essential for the generation of the mu-alpha and mu-beta components of the SI mu rhythm. Because our network was trying to reproduce in vivo data from awake humans, we also included exogenous driving inputs to the SI network in an FF and FB manner (Fig. 1). In the following text, we show that, whereas the intrinsic currents necessarily define the time constant of integration along the dendrite, the timing and strength of the excitatory FF and FB inputs determine the mu-alpha and mu-beta expression.
To simulate the single-trial spontaneous SI MEG mu rhythms, as shown for example in Fig. 6A, we first drove the SI network with a stochastic approximately 10-Hz burst of FF rhythmic drive, to see whether the mu rhythm would emerge from this simulated rhythmic lemniscal thalamic input combining with intrinsic cellular and local network properties. Details of the driving sequence are shown in Fig. 2 and described in methods.
The spectrogram from the simulated data in Fig. 6B shows that approximately 10-Hz FF rhythmic input alone created a primarily 10-Hz response in the model. The corresponding waveform shows that each cycle of the FF input created an upward current flow in the pyramidal neuron dendrites, resulting in a solely positive-polarity oscillation that lasted about 100 ms. This waveform contrasts with MEG data, where the polarity of the signal oscillates around zero (compare waveforms in Fig. 6, A and B; see also Fig. 4), suggesting that the current flow oscillates up and down the pyramidal neurons. Importantly, consistent with the notion that large-scale MEG/EEG ongoing oscillations represent the subthreshold activity of large numbers of synchronous neurons (Zhu et al. 2009), the driving parameters in our model were all chosen so that the oscillations in the pyramidal neurons remained subthreshold.
Our previous modeling study (Jones et al. 2007) showed that negative-polarity signals could be produced by FB excitatory input to the apical dendrites of pyramidal neurons lying in the superficial layers that induce a downward current flow. Therefore emulating approximately 10-Hz FB, a second stochastic 10-Hz burst of rhythmic excitatory synaptic input was provided to the apical dendrites of the layers II/III and V pyramidal neurons. The statistics of this FB drive were analogous to the FF drive, as described in methods and shown schematically in Fig. 2. Approximately every 100 ms, an FF burst was followed by an FB burst. In the example shown in Fig. 6C, on each cycle, the FF and FB inputs were nearly synchronous (Gaussian mean difference, 5 ms; FF and FB weights, 0.4 pS; n = 10; variance, 400). See the following text for further details on the importance of phasic latency and other input statistics (Fig. 8). Under these conditions, the 10- and 20-Hz components of the mu rhythm were robustly reproduced and the waveform of the rhythm oscillated around zero, with a subset of nonsimultaneous intervals of strong mu-alpha and mu-beta power, as observed in our MEG data (compare Fig. 6C with Fig. 6A and Fig. 4). The parameters regulating the mu-alpha and mu-beta components are described in detail in the following section (Fig. 7).
To test whether the waveforms are symmetric around zero, we calculated a symmetry index (SInd) for the MEG (Fig. 6D) and model data (Fig. 6E). The SInd was derived as described in methods and is such that a positive SInd value indicates greater-amplitude peaks, a negative value indicates greater-amplitude troughs, and a zero value indicates that an oscillation was symmetric around zero. Figure 6D (top panel) depicts the mean and SD of the SInd for each of the 10 subjects (n = 200, 1-s prestimulus trials); the bottom panel shows a histogram of the SInd across all subjects and trials, with the inset showing the mean and SD of the histogram. The SD of each individual subject overlaps zero, with four subjects >0 and six <0. The SInd was not significantly different from zero across the subjects (P < 0.001, t-test). Similarly, Fig. 6E shows a histogram of the SInd across trials in the model (n = 40, 1-s trials, parameters as in Fig. 6C) and inset with mean and SD of the histogram. As in the MEG data, the SInd index in the model was not significantly different from zero (P < 0.001).
In the model, prominent mu-beta cycles emerged when the stochastic FB input came in at the proper moment and strength to cut the FF alpha cycle in half (compare red boxes in Fig. 6, B and C), creating an oscillation that looks similar to that in the MEG data (Fig. 6A). As described historically in recordings of the MEG mu rhythm, our SI signal waveforms exhibit a “comb-like” or “arch-shaped” signal (Tiihonen et al. 1989), where the upward deflections appear more rounded, whereas the downward deflections have sharp edges forming what looks like several connected arches, as shown in Figs. 4 and 6A. This feature is also qualitatively reproduced by our model data (Fig. 6C) and arises as follows in the model. The 10-Hz FF input produces positive upward deflections that last about 100 ms (i.e., the rounded portion of the “arch”; Fig. 6B, right). The additional nearly 10-Hz FB input acts to push current flow in the opposite direction. In the example shown in Fig. 6C, the 10-Hz FF and FB inputs arrive nearly synchronously (5-ms delay); thus for part of the upward deflection the FB input works directly against the FF input, whereas at the end of the upward deflection, the FB input is able to abruptly push current flow down the dendrites, creating the “sharp” downward edges of the arch.
To create oscillations that were the same magnitude as that of oscillations seen experimentally, on the order of 100 nAm (Fig. 6A, right), the modeling data were multiplied by 30,000 (Fig. 6. B and C, right), which suggests that the subthreshold activity of nearly 200 × 30,000 = 6 million layers II/III and V pyramidal neurons contributed to the observed mu rhythm. Based on the general estimate of about 75,000 PNs per about 1 mm2 of the cortical sheet (Cheung et al. 2007), these findings suggest that in this study, about 80 mm2 (0.8 cm2) of cortical space was synchronized to generate the observed rhythms. This volume estimate would lie within the volume of the human hand representation, defined as the omega-shaped bend in the postcentral gyrus (Moore et al. 2000). This prediction is significantly larger than the number of neurons predicted to contribute to the tactile-evoked response in our previous and current modeling results shown in Fig. 9A, which was on order of 60,000 spiking neurons (200 PNs × 300 scaling factor = 60,000 PNs).
Statistical properties of the mu rhythm on single trials.
To further validate the model prediction that the SI mu rhythm is created by nearly synchronous, stochastic, alternating 10-Hz FF and 10-Hz FB inputs, we quantified other properties of the simulated mu rhythm in a manner analogous to the MEG data. In each case, we found consistency between the MEG and model data.
First, we calculated the mean and SE of the PSD for frequencies from 1 to 60 Hz across 100 1-s trials of the stochastically simulated spontaneous mu rhythm (Fig. 7A, n = 50 trials; all parameters fixed as in Fig. 6C). Peaks of activity emerge in the mu-alpha and mu-beta range. As in the grand average MEG data, the mu-alpha power has a higher peak in this PSD analysis, but a second smaller peak in the mu-beta range is also present (compare Fig. 7A with Fig. 3A).
Second, we calculated the probability that, for a given 100-ms time window, mu-alpha and mu-beta were in the top third of trials sorted from low to high power in the model. We found that they co-occurred roughly 50% of the time (Fig. 7B, inset, 48.7% calculated from 100 100-ms time windows), analogous to the mean of our grand-average MEG data (compare with Fig. 5, inset).
Third, we plotted a histogram of the ratio of alpha to beta power calculated in 100-ms time windows (Fig. 7B, grand total 1,000 prestimulus time windows, 100 bins). We found that the distribution of this histogram had the same qualitative shape as that calculated from the MEG data (compare Figs. 7B and 5). Further, as in the MEG data, the mean (1.4) and median (1.1) of the histogram were >1. A further examination of the model parameters controlling the relative dominance of alpha and beta power is given in the following text (Fig. 8).
Parameters controlling relative mu-alpha and mu-beta power dominance.
We varied the statistics of the input patterns shown schematically in Fig. 2 and detailed in methods, to examine their effect on the dominance of mu-alpha or mu-beta power expressed in the spectrograms generated from the model data. We began by parametrically varying the mean delay between the approximately 10-Hz FF and FB input, from 0 to 95 ms at 5-ms intervals (Fig. 8A; all other parameters were fixed as in Fig. 6C; also see methods). We calculated the mean ratio of mu-beta/mu-alpha power over a 1-s simulation and found that when inputs arrived nearly synchronously with FB following FF (delay <10 ms), or vice versa, the relative dominance of mu-alpha and mu-beta power was approximately equal (mu-beta/mu-alpha power ≃ 1; see endpoints of parabola in Fig. 8A). In contrast, when the alternating approximately 10-Hz FF and FB inputs were asynchronous (delay = 50 ms), the mu-alpha rhythm dominated the signal (mu-beta/mu-alpha ≃ 0). In this case, the FF input created a positive-polarity signal for half of a 10-Hz cycle and the FB input created a negative-polarity signal for the remaining half of a 10-Hz cycle, where the peaks of the two halves are 50 ms apart, defined by the 50-ms delay (data not shown).
We then investigated the influence of the variance (relative synchrony in presynaptic spiking, set by the SD of the “input bursts”), amplitude (measured as number of “input bursts”), and efficacy (differences in the postsynaptic conductance evoked by an identical presynaptic input) of the separate FF and FB inputs on the relative mu-beta to mu-alpha power (Fig. 8, B, C, and D, respectively; each parameter was tested at the three values on the x-axis; see methods and Fig. 2 for definition and illustration of each parameter). In each case, manipulations that enhanced the FF drive [increased synchrony (decreased variance), amplitude, and postsynaptic conductance] increased mu-alpha power (red curves), and manipulations that enhanced FB drive increased mu-beta power (blue curves).
Each of these mechanisms for increasing drive to the network [increased synchrony (decreased variance), amplitude, and efficacy of postsynaptic input] had the net effect of increasing the size of the net postsynaptic conductance. Based on this parameter search, in the remainder of our modeling investigation we simulated high mu power by simultaneously increasing the postsynaptic conductance of the FF and FB inputs, thus producing an increase in mu-alpha and mu-beta power (Fig. 8E). Low mu power was simulated with smaller postsynaptic conductance (high mu 0.6 pS to PNs and 1.2 pS to Ins; low mu 0.4 pS to PNs and 0.8 pS to INs).
Impact of mu power on tactile-evoked responses.
Our previous modeling results indicated that the SI-evoked response measured with MEG was created by a specific temporal sequence of exogenous excitatory inputs to the SI network (Jones et al. 2007). This sequence consisted of FF input at ≃25 ms poststimulus, followed by FB input at ≃70 ms, and by a late feedforwad (LFF) input at ≃135 ms (shown schematically in Fig. 9A, where each input consisted of single presynaptic spikes per trial; see methods). The sequence that induced the SI-evoked signal may be interpreted as initial FF input from the periphery through the lemniscal thalamus to granular layers, followed by FB input from higher-order cortex or nonspecific thalamic sources to the supragranular layers, followed by a second wave of lemniscal thalamic input to the granular layers.
Here, we used our expanded model to investigate the effects of prestimulus mu power on SI-evoked response magnitude and timing. The evoked response input sequence was delivered at several starting phases (SPs) within simulated high and low prestimulus mu-alpha and mu-beta cycles (marked with black bars in Fig. 8E) and results were averaged (Fig. 9A). The SP is defined as the mean time of the initial lemniscal thalamic FF input to the SI network, which was followed by FB input 45 ms later (mean SP +45 ms) and a subsequent LFF input 65 ms later (mean SP +65 ms). This sequence was delivered at 20 different equally spaced SPs in the marked mu-alpha cycle and 10 different equally spaced SPs in the marked mu-beta cycle; the results were averaged over all 30 trials (see Table 2 for default postsynaptic conductances and input times of evoked response inputs).
When given during either high or low prestimulus mu power, the average simulated evoked response created a waveform with a negative-polarity peak at about 70 ms (M70) and two positive-polarity peaks at about 100 ms (M100) and about 135 ms (M135), respectively (Fig. 9A). The timings and polarities of the peaks are similar to those generated in our previous study, in which prestimulus rhythms were not considered (Jones et al. 2007).
Several significant differences emerged in the evoked response waveforms simulated during low and high prestimulus mu conditions (compare dark and light blue curves in Fig. 9A; red stars indicate time points where the difference was P ≤ 0.05; purple stars: P ≤ 0.01, paired t-test, n = 30 trials). The greatest difference was the emergence of a positive peak near approximately 50 ms, labeled M50, under high prestimulus mu that was negligible under low mu (M50 = max 40–60 ms, P < 0.01). In addition, in the model, the magnitude of the M70 was smaller (M70 = min 50–100 ms, P < 0.05). There were also slight but significant differences between the M100 and M135 peaks, and later [150–175] ms activity such that the magnitude was larger under low mu conditions (M100 = max 75–125; M135 = max 125–175, mean 150–175 ms, P < 0.05; see Fig. 9A). Predictions as to the underlying neural activity creating these differences are discussed in the following text (Fig. 9C).
We next related these predicted differences to the SI-evoked response waveforms in the MEG signal and found consistency between the model and experimental data. Figure 9B shows the MEG SI-evoked response sorted over high and low prestimulus mu trials for an equal number of hit-and-miss trials comprising a 50% detection rate (n = 10 Ss, 30 high and low mu trials per subject). We found that the prestimulus mu power significantly influenced the early evoked response (<70 ms), and showed trends to a difference in the later components. As in the model data, a significant positive peak near 50 ms, labeled M50, emerged under high prestimulus mu (Fig. 9B, P < 0.05 marked with red asterisks). This peak was followed by a trend that approached significance for a decreased M70 peak under high mu (M70 = min 60–75 ms, n = 7/10 Ss, P = 0.1 sign-test), followed by a late trend toward larger magnitude under high mu (mean 150–175 ms, n = 8/10 Ss, P = 0.055 sign-test).
The late trend to higher-magnitude response under high mu conditions and the nonreturn to baseline between the M100 and M135 peaks, in the MEG data are discrepant with the model results. The nonreturn to baseline was also observed in our previous study (Jones et al. 2007) where the strengths of the FF and FB inputs in the current study were derived. One way that this difference could be rectified in the model is by assuming that the late FF input at ≃135 ms synchronously drives a greater number of neurons. This effect can be achieved in the model by multiplying by a separate increased scaling factor at the late time points (>100 ms; data not shown). A large late FF input that occurs ≃135 ms after the stimulus may arise as part of an induced poststimulus approximately 10-Hz thalamocortical oscillation, which may be stronger with high prestimulus mu. A strong induced 10-Hz oscillation is consistent with the notion that the stimulus “phase resets” the 10-Hz component of the oscillation (Hanslamayer et al. 2007; Makeig et al. 2002).
Increased excitation and inhibition during high mu influences the M50 and M70 and evoked response peaks.
To better understand how increased prestimulus drive during high-prestimulus mu rhythms created an initially larger M50 response followed by a smaller M70 and later response between the M100 and M135 peaks, we examined the evoked spiking in the network model. We plotted the mean number of spikes (smoothed over 5-ms bins) for the pyramidal and inhibitory neuron populations in our network during high and low prestimulus mu power (Fig. 9C, pooled responses, 30 trials). We found that when the evoked response was simulated during high mu power (light blue curves), the FF input at ≃25 ms induced more spiking in pyramidal and inhibitory neurons in layers II/III and V compared with low mu (dark blue curves). Firing in the pyramidal cells to the initial FF approximately 25 ms input peaked near 50 ms, creating the initial M50 positive peak seen in the evoked response following high prestimulus mu (light blue curve, Fig. 9A). The positive polarity of the M50 peak comes from back-propagation of action potentials up the apical dendrites of the pyramidal neurons (Jones et al. 2007). The strong early evoked inhibition suppresses the pyramidal response to subsequent FB input at about 70 ms, creating a decreased response and less spiking under high mu conditions (compare <70 ms dark and light traces for L5 and L2/3 excitatory cells in Fig. 9C). In turn, the decreased activity in the E cells created a smaller M70 and later responses in the model (compare dark blue and light blue curves in Fig. 9A).
We used MEG imaging and biophysically principled computational neural modeling to provide a detailed characterization of the mu-alpha and mu-beta components of the spontaneous SI mu rhythm, investigate its neural origin, and investigate its impact on early (0–175 ms) evoked responses. The MEG recordings revealed an SI mu rhythm with mu-alpha and mu-beta components, whose relative expression was subject dependent. On single trials, the mu-alpha and mu-beta components were not simultaneous in their emergence, but did co-occur at rates greater than chance.
Our model results led to the novel prediction that the SI mu rhythm arises from a stochastic nearly synchronous alternating sequence of approximately 10-Hz FF followed by FB excitatory input to SI. Simulating mu in this manner reproduced many key features of the MEG data, including: a mu rhythm with both mu-alpha and mu-beta components; a rhythm that oscillates symmetrically around zero polarity; nonsimultaneity of mu-alpha and mu-beta components in time; and an enhanced early (M50) peak in the tactile-evoked response under high mu.
The predicted enhancement in the early response when stimulus presentation coincided with epochs of high mu power was observed in our data. We note that the model was completed prior to running our analyses of rhythmogenesis or evoked response modulation and that the only model parameters available to manipulation during analysis were those described in the text (primarily the strength and arrival synchrony of FB and FF signals). The convergence of experimental and model analysis given this principled modeling approach supports the predictive value of the neural interpretation of model behavior.
To our knowledge, our results provide the first evidence that the mu rhythm arises from the combination of two stochastic approximately 10-Hz rhythms that arrive to the SI network at different laminar locations corresponding to thalamic FF input and intracortical FB input. Because the rhythmic approximately 10-Hz FF and FB inputs and SI cortical structure in our model are defined by general principles of cortical circuitry (Felleman and Van Essen 1991; Hughes and Crunelli 2005), the predicted mechanisms may be applicable to the generation of mu-alpha and/or mu-beta frequency rhythms observed in other cortical areas. Table 3 summarizes the results of our MEG findings, the corresponding predictions of our biophysical model, and the mechanistic interpretation of each finding from the model.
Importance of feedforward and feedback input
In expanding our previous SI model (Jones et al. 2007) to study mu rhythms, we included several pyramidal neuron intrinsic currents (see introduction and methods) and local excitatory and inhibitory interactions, hypothesizing that these factors would prove important for the generation of the mu-alpha and mu-beta components of the SI mu rhythm. Because our network was trying to reproduce in vivo data from awake humans, we also included exogenous driving inputs to SI in an FF and FB manner (Fig. 1). In contrast to our initial assumptions, we found that although the intrinsic current kinetics in the pyramidal neurons were essential to the time constant of the subthreshold integration of current flow along the dendrites, the dominant time constants regulating mu-alpha and mu-beta frequency in the model were determined by the relative strength and timing of the exogenous rhythmic excitatory inputs that drive current flow up and down the pyramidal neuron dendrites (Fig. 8). These exogenous inputs were necessary to accurately reproduce 10- and 20-Hz components and distribution of the awake human SI mu rhythm.
Intuitively, one might predict that alternating 10-Hz inputs at a 50-ms delay would ideally combine to generate a 20-Hz MEG signal. However, our results suggest that the inputs must arrive nearly simultaneously (<10 ms) for 20-Hz cycles to emerge. In this case, approximately 10-Hz current flow that is driven up the pyramidal neuron dendrites (via FF input) is cut in half by current flow driven down the dendrites (via FB input) to create a 20-Hz cycle. With perfectly asynchronous (50-ms delay) FF and FB inputs, a dominant mu-alpha rhythm is produced. This result, and the fact that stronger FF inputs increase mu-alpha, whereas stronger FB inputs increase mu-beta (Fig. 8, B–D), leads to direct predictions as to the source of variability in mu-alpha or mu-beta expression across subjects (Fig. 3B). We may surmise that subjects with dominant mu-alpha components possess stronger rhythmic FF input and/or that the intracortical FB input arrives with a phase delay near 50 ms. In contrast, we predict that subjects with dominant mu-beta components possess stronger FB input and/or nearly perfectly aligned (<10-ms delay) FF and FB inputs.
Subthreshold oscillations and their relation to higher-frequency gamma rhythms
Several studies report that low-frequency alpha and/or beta range rhythms are directly coupled to higher-frequency gamma (40–80 Hz) rhythms. This effect has been observed in vivo in humans in electrocorticographic data (Canolty et al. 2006), MEG data (Palva et al. 2005a), and intracranial cortical LFP recordings (Lakatos et al. 2005; Schroeder and Lakatos 2009b). In our data, we do not see prominent gamma activity in the prestimulus time period (Fig. 3), likely due to the fact that our MEG signal measures the subthreshold synchronous activity of a large number of pyramidal neurons (Hamalainen et al. 1993; Zhu et al. 2009), estimated to be on the order of 6 million. Gamma rhythms are likely produced by a smaller subnetwork of spiking excitatory and inhibitory neurons (Cardin et al. 2009; Kopell et al. 2000; Pinto et al. 2003; Vierling-Claassen et al. 2008; Whittington et al. 2000) and are thus not recorded in our signal. Further, due to the limited size of our network, the model cannot currently simultaneously reproduce both large- and small-scale phenomena.
Coordination of ongoing 10- and 20-Hz activity between brain areas
The prediction that the SI mu rhythm is driven by approximately 10-Hz rhythmic FF and FB in the model implies that the SI mu rhythm is correlated with 10-Hz activity in lemniscal thalamic sources (providing FF inputs), as well as in higher-order cortical sources (providing FB inputs) and/or nonlemniscal thalamic input. These hypotheses are consistent with several studies that report the existence of coherence in power (Schubert et al. 2008; Zhang and Ding 2009) and phase-locking (Hanslmayr et al. 2007; Palva et al. 2005a) between SI mu rhythms and alpha and beta frequency activity in prefrontal cortex and other higher-order cortical areas. Most significantly, Zhang and Ding (2009) recently applied Granger causality analysis to 10-Hz activity measured from EEG electrodes above somatosensory and prefrontal cortices and found that 10-Hz activity propagates from prefrontal cortex to SI. Palva et al. (2005) observed local cross-frequency phase synchrony between alpha- and beta-band oscillations widely across the cortex, particularly over the somatomotor regions, during rest conditions that increased with execution of mental arithmetic tasks (see also Nikulin and Brismar 2006). Interareal 10- to 10-Hz phase synchrony (and even more so 10- to 20-Hz synchrony) increased with task demands with the strongest effects at long distances. Hanslmayr et al. (2007) examined interareal within-frequency phase synchrony between EEG electrodes during a visual-perception task and found prestimulus alpha (8–12 Hz) phase synchrony mainly between frontal and parietal electrode sites that decreased with perception. Beta (20–30 Hz) phase synchrony was also observed and increased with perception. Schubert et al. (2008) used EEG during a somatosensory perceptual masking task and found that prestimulus beta activity (18–26 Hz) in prefrontal and sensorimotor electrodes covaried such that they were both lower during perceived trials. These results emphasize the cooperation of alpha- and beta-band activity across the brain and support the prediction that the SI mu rhythm may arise from rhythmic interactions of 10-Hz oscillations projected to SI in an FB manner from higher-order areas. However, our FB inputs could also arise from nonspecific thalamic inputs oscillating with a slight phase shift from that of the lemniscal driver nucleus (Guillery and Sherman 2002; Hughes and Crunelli 2005; Jones 2001). The prediction of near synchrony in inputs necessary to generate beta may be better explained by a common thalamic source, which may represent a different form of “feedback” or may represent a distinct form of input.
Symmetric oscillations and their relation to evoked responses
Our finding that the prestimulus SI mu rhythm oscillates symmetrically around zero is in contrast to previous studies that have observed nonsymmetric alpha-frequency oscillations that contribute to the generation of late (>250 ms) components of sensory-evoked responses. Nikulin et al. (2007) reported baseline shifts in mu-alpha frequency activity in sensorimotor MEG sensors on the scalp. These shifts in the mean of the oscillation over time were reported to be responsible for late (>250 ms) components of median nerve-evoked response. In another study, Mazaheri and Jensen (2008) observed asymmetric amplitude modulations in posterior alpha activity that was proposed to contribute to slow components (>300 ms) of visual-evoked responses over occipital cortex. Both studies conjectured that the baseline shifts and asymmetric oscillations were likely due to differences in outward and inward neural currents producing the magnetic field, which they posit are unlikely to be equal.
Our data and modeling results suggest that the net inward and outward current flow, within pyramidal neurons across the SI network, are approximately balanced, creating a symmetric index near zero in our signal (Fig. 6). The discrepancy between our results and those of Nikulin et al. (2007) and Mazaheri and Jensen (2008) may arise from several areas. These studies showed that differences in the later components of evoked sensory responses (>250 ms) were tied to asymmetry in the prestimulus alpha oscillation. In our study, we investigated only the evoked response components <175 ms. Notably, in our data there is a late trend (135–175 ms) to a larger MEG SI-evoked response under high prestimulus mu that was not reproduced by the model. We conjecture that this trend may arise from stronger poststimulus approximately 135-ms FF inputs to the network during high prestimulus mu conditions. This timing (∼100 ms after the first evoked FF input arrives to the network) is consistent with the idea that there is also a stronger poststimulus evoked 10-Hz oscillation that is perhaps phase-locked (Hanslamayer et al. 2007; Makeig et al. 2002). Investigation of evoked poststimulus oscillations is not a focus of our study and we do not make any strong claims regarding these phenomena here.
Further, we studied activity localized to a primary equivalent current dipole in SI, whereas the other studies investigated sensor data. Both approaches have strengths. A benefit of dipole localization is that it provides a more concrete localization of a single signal source. Mazaheri and Jensen (2008) described the presence of a bipolar field pattern in their asymmetry index measure that was indicative of the existence of a current dipole located between the positive and negative signals. However, they did not investigate the location of the proposed dipole in the brain or find a consistent direction of the dipole current across subjects, from which inward and outward could be defined.
Last, we calculated the symmetry index of our signal across the entire frequency band of the SI mu rhythm complex (7–29 Hz), whereas the previous studies investigated asymmetries in activity from signals that had been band-passed in the alpha or beta range. Nikulin et al. (2007) used an independent component analysis to filter sensor data in the range of 8–13 Hz and investigated baseline shifts in the mean amplitude of the filtered signal. Inspection of our signal over long time periods (10 s) showed modulation of the mu rhythm without such a dc-offset effect (see methods). Mazaheri and Jensen (2008) band-passed their data using narrow width intervals from 5 to 40 Hz (e.g., 8–12 Hz for alpha activity) and found peaks and troughs in the band-passed data. The time points of the peak and trough values in the band-passed signal were used to calculate the peak and trough values in the original signal from which the symmetry index was calculated.
Our analysis methods for finding peaks and troughs in the data are analogous to those of Mazaheri and Jensen (2008). However, in Fig. 6 we quantified the symmetry index from peaks and troughs calculated from a broader band filter of our signal that contained the entire mu range (5–30 Hz) (see methods). In an additional analysis, we quantified the symmetry index from calculated peaks and troughs in our signal by first band-passing in narrow width bands (3-Hz widths) from 5 to 40 Hz, precisely as in Mazaheri and Jensen (2008). We found the symmetry index to be significantly different from zero at only 5–7 Hz in that analysis of the data (P = 0.03, t-test; data not shown). However, closer inspection of the calculation of peaks and troughs from the data filtered in this small band showed that this is not the most accurate way to calculate these values in our signal. This method produced peaks and troughs that did not align with the true values in the raw data. In contrast, a band-pass that contained the entire mu range (5–30 Hz) produced precise alignment of peaks and troughs in the band-passed and raw data, thus providing a more accurate measure of symmetry in our signal. Supplemental Fig. S1 shows examples of single-trial (1-s) raw-data waveforms from the MEG and model SI signal, along with the band-passed versions of these signals in a narrow low-frequency band from 5 to 7 Hz (Supplemental Fig. S1, A and C) and in a band containing the entire mu range, 5–30 Hz (Supplemental Fig. S1, B and D).1 Calculated peaks and troughs are shown on the band-passed and raw data with red and green asterisks, respectively. In the case of calculating the peaks and troughs from the 5- to 7-Hz filtered signals, the calculated peak and trough values do not line up with those apparent in the raw data and, in several instances, are flipped such that troughs are higher than peaks (Supplemental Fig. S1, A and C). In contrast, when calculating the peaks and troughs from the broader-band 5- to 30-Hz filtered signal, there is clear alignment between the calculated values and those in the raw data (Supplemental Fig. S1, B and D). The misalignment in peaks and troughs when using only a narrow low-frequency band is likely due to the fact that our SI dipole signal is a true two-component signal that contains prominent alpha and beta activity and band-passing the signal in a small frequency range distorts this fact. The observed misalignment in peaks and troughs would be less pronounced in signals that contained a predominantly 10-Hz oscillation, as in Mazaheri and Jensen (2008).
Motor gating of sensorimotor mu rhythms and somatosensory-evoked responses
Our study focuses on the sensorimotor mu rhythm and tactile-evoked responses exhibited by a primary current dipole localized to the hand representation of SI, without considering motor activity. Nevertheless, there is synergy between our results and the literature examining movement-induced decreases in sensorimotor mu rhythms and gating of somatosensory-evoked responses. Movement desynchronizes sensorimotor mu rhythms (Neuper et al. 2006; Pfurtscheller et al. 1997; Salenius et al. 1997). Several studies have also shown that limb and digit movements attenuate early components (≤45 ms) of electrically induced (typically median nerve) somatosensory-evoked potentials (Cheron and Borenstein 1987; Cohen and Starr 1987; Nishihira et al. 1997; Rossini et al. 1999; Rushton et al. 1981; Tapia et al. 1987). Direct electrophysiological recordings in animal models suggest that this attenuation is due to decreases in neuronal spiking activity, most probably of pyramidal neurons (Chapin and Woodward 1982; Ro et al. 2000). Rossini et al. (1999) used high-resolution EEG to localize this effect in human somatosensory, motor, and supplementary motor areas and found that in each area the gating effect was strongest at peaks between 30- and 45-ms post-nerve stimulation.
Our results show that decreased mu in SI predicts a decreased M50 peak, occurring with a maximum between 40- and 60-ms post-tactile stimulation. Our model predicts that this decrease is due to lower neuronal firing rates in SI pyramidal neurons (Fig. 8). The decreases in mu in our study are spontaneous and not directly related to motor movement; the attenuated peak occurs at a slightly longer latency than that in the described studies, likely due to the slower propagation time to the cortex from a brief tactile stimulation versus a strong electric stimulation to the nerve. However, the parallels in the two phenomena—decreased mu power predicting decreased evoked responses—suggest that similar circuits and mechanisms are common. In the following text we discuss how prestimulus mu-related changes in evoked activity may also be connected to facilitating detection of tactile-evoked responses.
Influence of the prestimulus mu rhythm on evoked response gain and implications for detection
In the present study, we intentionally excluded analysis of the impact of mu oscillations on detection. The focus of the present study was on mu—its rhythmogenesis and its impact on the evoked response—and the additional domain of relating these findings to perceptual success requires its own extensive treatment. Specifically, to determine the effects of mu power on the SI-evoked response, independent of perceptual success, trials were chosen to comprise a 50% detection rate for high and low mu conditions. That said, our results show that the changes in the evoked response driven by high or low mu conditions are different from those we showed in a previous study to predict detections. Our modeling results predicted that the ongoing mu rhythm influenced early components of the SI-evoked response (≤70 ms), with a weaker effect on later evoked activity ≤175 ms (Fig. 9A). This hypothesis was confirmed in the MEG data when sorting trials over high and low prestimulus mu power. Under high prestimulus mu conditions, a significant early M50 positive peak emerged, followed by a trend toward a weaker M70 response (Fig. 9B). In our previous study (Jones et al. 2007), we found that the later components of the SI-evoked response, beginning at about 70 ms (M70), were correlated with detection of a threshold-level tactile stimulus, such that the response magnitude was greater on detected trials. The model predicted that the late differences with detection reflected poststimulus changes in higher-order cortex or nonspecific thalamic projections that created earlier and stronger evoked FB at ≃70 ms and LFF input at ≃135 ms to SI on detected trials (see Jones et al. 2007). Several studies have reported a nearly inverse relationship between mu power and detection (Hanslmayr et al. 2007; Linkenkaer-Hansen et al. 2004; Schubert et al. 2008; Zhang and Ding 2009). In light of these findings, we may conjecture a possible relationship between low prestimulus rhythms in SI and “higher-order areas” on detected trials. Low-power ongoing rhythms in the higher-order area, producing weak ongoing approximately 10-Hz FB to SI, have greater early (<70 ms) poststimulus-evoked responses that in turn create greater FB to SI at nearly 70 ms on detected trials.
In summary, our study is the first to characterize the interdependence of the mu-alpha and mu-beta components of the SI MEG measured mu rhythm on single trials. Our computational model provides a novel mechanistic interpretation of the mu rhythm and its functional significance and presents the first evidence that the SI mu rhythm may arise from a stochastic nearly synchronous alternating sequence of approximately 10-Hz FF followed by FB input to the SI network. Importantly, specific predictions of our computational model are in agreement with the characteristics of the SI activity based on source modeling of our MEG data.
This work was supported by National Institutes of Health Grants P41-RR-14075, K25-MH-072941, 1RO1-NS-045130-01, and T32-GM-007484; National Science Foundation Grant 0316933; the Athinoula A. Martinos Center for Biomedical Imaging; and the McGovern Institute for Brain Research.
We thank D. Ziegler and P. Hosseini-Varnamkhasti in the Brain and Cognitive Science Department at Massachusetts Institute of Technology for excellent technical support in implementation of data analysis methods, D. Vierling-Claassen and Q. Wan for thorough remarks on the manuscript, and M. Hines for continued support in NEURON software coding.
↵1 The online version of this article contains supplemental data.
- Copyright © 2009 the American Physiological Society