Neuronal Mechanisms of Cortical Alpha Oscillations in Awake-Behaving Macaques

Surgical preparation.

Monkeys were prepared for chronic awake intracortical recording. Surgery was performed using aseptic techniques, under general anesthesia, as described previously (Givre et al., 1994; Schroeder et al., 1995, 1998; Chen et al., 2007; Lakatos et al., 2008). Briefly, the tissue overlying the calvarium was resected, and appropriate portions of the cranium were removed. The neocortex and overlying dura were left intact. To provide access to the brain and to promote an orderly pattern of sampling across the surface of visual areas, recording chambers incorporating parallel guide tube grids (Crist Instrument) were positioned normal to the brain surface for orthogonal penetration of (1) portions of area V2 underlying the foveal representation of area V1 in the lateral striate operculum, (2) foveal/parafoveal representations in area V4 lying on the lateral prelunate gyrus, and (3) foveal/parafoveal regions of IT cortex lying on the lower bank of the superior temporal sulcus, anterior to the MT/MST complex (medial temporal/medial superior temporal cortex) and lateral to FST (fundus of the superior temporal sulcus). Choice of these specific portions of the target regions enabled us to penetrate the target regions orthogonal to their lamination patterns to optimize the implementation of one-dimensional CSD analyses (see below). Implantation was guided by stereotaxic transformation of magnetic resonance imaging data, which delineated the cortical gyral pattern. Individual epidural guide tubes were positioned over central and frontal sites to serve as ground and reference electrodes. Together with socketed Plexiglas bars (to permit painless head restraint), they were secured to the skull with orthopedic screws and embedded in dental acrylic. A recovery time of 2 weeks was allowed before the beginning of data collection.

Paradigm.

As part of an experiment involving switching attention between auditory to visual input streams, two male macaque monkeys (Macaca fascicularis), referred to as B and V, were trained to perform an auditory oddball discrimination task (Lakatos et al., 2008). In this task, 100-ms-duration pure tones were presented at 1.5 Hz, with the stream of “standard” stimuli randomly interrupted by tones that differed in frequency (deviants). The monkey was required to respond to the deviants after the onset of an oddball stimulus. A liquid reward was given to correct responses. The reason for analyzing activity in visual cortices during auditory discrimination was that the discrimination kept the monkey verifiably alert without using visual stimuli, so that we could study spontaneous neural activity in V2, V4, and IT cortex. In addition, the behavioral performance measured by the reaction time provided an opportunity of examining the functional significance of alpha activity in various visual cortices. Data were collected during periods of adequate task performance (i.e., >80% target detection).

Recordings.

Animals sat in a primate chair in a completely dark, electrically shielded, sound-attenuated chamber with head fixed in position, and were monitored with infrared videography. LFP and MUA were sampled (2 kHz) with a linear array electrode with 14 contacts spanning all six cortical layers (see Fig. 1A for a schematic). The electrodes were constructed by threading 25-μm-diameter, Teflon-coated wires through a stainless-steel tube, aligning and spacing the wires in a linear array at the distal tip of the tube (beveled), embedding the wires in fast-curing epoxy, and cutting the wires flush with the epoxy after it dried. The resulting line of recording contacts was thus 25 μm in cross section, and had impedances between 0.3 and 0.5 MΩ at 2 kHz. The intercontact spacing used was 150 μm in V2 and 200 μm in V4 and IT. Multiple penetrations were made in each of the three visual cortical areas. The dataset analyzed here consists of two penetrations in V2, two penetrations in V4, and three penetrations in IT for subject B; and two penetrations in V2, two penetrations in V4, and one penetration in IT for subject V.

Functional position of multielectrodes.

For each experiment, electrodes were lowered into chambers targeting one or more of the aforementioned visual areas. Using the intracortical profiles of LFPs and MUA generated in response to visual stimulation, the electrodes were guided to the target region and then positioned so as to straddle its laminar expanse as described in detail previously (Schroeder et al., 1998) [see also more recent studies that have used this approach (Lipton et al., 2006; Chen et al., 2007; Lakatos et al., 2007)]. The positioning method is based on the fact that all visual areas recorded to date generate an LFP distribution with a constant amplitude decay (i.e., second derivative of approximately zero) over distance (see also Schroeder et al., 1995). Visually evoked LFP and MUA profiles are collected at successive depths beginning well in advance of the expected depth of the target region. Linear amplitude increase (near-zero second derivative) indicates approach to an active visual area. Arrival at the depth of the active region is indicated by steep acceleration in the LFP gradients and/or component polarity inversions (nonzero second derivative). MUA is also useful, because in association with a relatively flat CSD profile, active MUA indicates that the electrode array is traversing a span of white matter. Sulcal cortical visual areas are distinguished from intervening white matter regions by the conjunction of large CSD components with multiunit activity. Because the CSD indexes the net local transmembrane current flow in a cortical region, the large, current-balanced CSD profiles characteristic of cortical regions and laminated nuclear structures are not found in association with the MUA recorded in white matter. After localizing an active region, further depth adjustments are made based on the CSD profile, so as to straddle its laminar expanse with the electrode array. Because the uppermost laminae in most cortical areas are relatively quiescent, we can distinguish between sulcal areas whose pial surfaces were opposed. In this case, the relevant distinction is between IT cortex (on the lower bank of the superior temporal sulcus) and the overlying cortical regions. In any case, the outer boundaries of an active cortical region, both pial and white matter regions, are indicated by drops to near-zero CSD values. Then, after the direction from which the tissue is penetrated (pial or white matter) is determined (by track reconstruction), approximate internal laminar assignments can be made, based on the known internal geometry of the region (see below).

Localization of V2, V4, and IT.

In both monkeys, visual areas were defined by anatomic location, subject to histologic verification [for detailed illustration, see Schroeder et al. (1998) and Chen et al. (2007)]. V2 recordings were made in the first belt of cortical tissue underlying the operculum, immediately posterior to the lunate sulcus, by penetrating through V1. V4 recordings were in the crown of the prelunate gyrus. IT recordings were made in the lower bank of the superior temporal sulcus.

Identification of laminae.

Once the electrode is positioned so that the contacts straddle the active current flow pattern of an area, gross laminar positions can be gauged from the known anatomical dimensions of each region. In extrastriate visual cortices, for example, layer 4 lies just below the midpoint of the laminar depth of the area (Hendry et al., 1990; Jones, 1990). In addition, our previous studies have revealed some functional clues to laminar locations. In particular, we have noted that in V2, V4, and IT (other extrastriate areas as well), the largest current sink associated with the N95 component of the visual evoked potential (this corresponds to a human N1) is located at the depth of layer 4; in fact, this current sink reflects the main feedforward activation of the area from lower visual cortices (Givre et al., 1994; Schroeder et al., 1995, 1998; Chen et al., 2007).

The data.

The neocortex is a structured tissue with high specificity of connections between its six layers (Lund, 2002). Whereas extracellular LFP, generated by transmembrane current flow, reflects mainly dendritic processing, the MUA gauges the intensity of pyramidal cell firing. Together, LFP and MUA form a comprehensive measure of cortical neuronal ensemble activity at both the input and output level (Nicholson, 1973; Mitzdorf, 1985; Schroeder et al., 1995). Previous work has shown that simultaneous recordings of LFPs and MUAs from different layers of the cortex are critical to controlled analysis of the dynamics of cortical circuitry under different behavioral states (Schroeder et al., 1995; Mehta et al., 2000a,b). It is generally accepted that the largest contributions to LFPs come from pyramidal neurons, which comprise nearly 85% of the cortex (Lopes da Silva and Rotterdam, 2005). Under certain conditions, additional neuronal elements are also known to contribute (Tenke and Schroeder, 1990; Schroeder et al., 1991, 1992, 1995). Typically, LFPs are recorded against a distant reference, making them susceptible to volume conduction of potentials from other sites. The first derivative used for generating the bipolar LFPs and the second-derivative approximation used for the CSD analysis help to eliminate this problem (Mitzdorf, 1985). Moreover, the CSD depth profile allows the localization of transmembrane current generators to specific cortical layers. Concomitant MUA helps assess the physiological processes (i.e., excitability changes) that current generators are associated with.

Logic of the analysis protocol.

We seek to characterize the laminar organization of the cortical alpha rhythm by analyzing depth LFPs and MUAs. The analysis protocol has three main steps. First, transmembrane current generators that give rise to field oscillations at the alpha frequency are identified with the CSD method. Traditionally, CSD analysis is performed on field potentials averaged with respect to a repetitive sensory stimulus (Mitzdorf, 1985; Schroeder et al., 1991, 1995). The CSD of ongoing neural activity is more difficult to ascertain. Single-trial CSD estimates tend to be noisy (Shah et al., 2004; Lakatos et al., 2005, 2007), and because there is no stimulus-related trigger, LFP averaging requires an alternate procedure for the alignment of trials. The phase-realigned averaging technique (PRAT; see below) is such a procedure that was used here to identify the alpha current generators within the cortical column. Second, alpha current generators that have the potential of pacemaking are identified with CSD–MUA coherence. In the context of studying evoked potentials, a source or sink is considered active if simultaneously recorded MUA is depressed or enhanced, indexing net local hyperpolarization or depolarization, respectively. For ongoing oscillatory activity, the membrane undergoes rhythmic transition between hyperpolarization and depolarization. It is thus more appropriate to consider the oscillatory current generator in terms of its pacemaking capability. During the depolarizing phase of the oscillation, the pacemaker cells fire bursts of action potentials, which, via synaptic transmission, entrain neural activity in other laminae and cortical areas. For the present work, significant phase coherence between CSD and MUA is used to indicate that a current generator is accompanied by rhythmic firing and thus has the potential of being a pacemaker. Third, the primary local pacemaking generator is identified with the Granger causality analysis. For a cortical column with multiple alpha current generators distributed across different layers, the relationship among these generators needs to be further delineated. This is particularly so if the second step reveals that more than one generator has the potential of being the pacemaker. Granger causality analysis, based on bipolar LFPs, is used to further disambiguate the roles of different current generators, because the primary local pacemaking generator is expected to exert unidirectional causal influence on other neural ensembles (see below).

Current source density analysis of ongoing oscillatory activity.

The length of a contiguous segment of spontaneous ongoing activity was on average 30 s long, and there were five to seven such segments for a given penetration. After high-pass filtering (3 Hz, zero phase shift) and downsampling to 200 Hz, the LFP data were further divided into epochs of 200 ms in duration, which were considered trials or realizations of an underlying stochastic process. To calculate the average current source density, the following procedure was used. (1) The power spectrum of each recording contact for a given electrode was estimated, and the contact showing the highest power spectral density at the alpha frequency was chosen as the “phase index contact.” (2) Sinusoid of the same frequency as the spectral peak in 8–12 Hz was then fitted to the data from the phase index contact for each epoch to obtain the phase at that frequency with respect to the beginning of the epoch. (3) The LFP data from all the contacts were shifted according to this estimated phase to realign all the trials, and the realigned signals were then averaged across epochs (trials) to obtain the averaged LFP for each contact. (4) A second spatial derivative was calculated to produce the CSD profile (Schroeder et al., 1995), from which the alpha current generators were identified. Csicsvari et al. (2003) has used a similar method to study the oscillatory activity in the hippocampus. For ease of reference, in the rest of the paper, the above procedure will be called the PRAT, and the resulting averaged LFP and CSD are referred to as PRAT-LFP and PRAT-CSD, respectively.

CSD and MUA coherence.

The MUA data were epoched the same way as the LFP data and downsampled from 2 kHz by taking a temporal average in nonoverlapping windows of 5 ms duration to achieve effectively the same sampling resolution of 200 Hz as the downsampled LFPs. Penetrations and recording contacts with poorly recorded MUA activity were excluded from further analysis. The coherence between single-trial CSDs around alpha current generators identified by the PRAT-CSD method and the corresponding mean-centered single-trial MUAs was calculated by the multivariate autoregressive (MVAR) spectral analysis method (see below). The coherence estimates were cross-checked by standard Fourier-based routines in Matlab (MathWorks). In addition, the phase spectrum between CSD and MUA at the alpha oscillation frequency was obtained to help assess whether, during an alpha epoch, a sink (or source) corresponds with increase (or decrease) in MUA.

MVAR spectral analysis.

The interaction between the current generators in different laminae was studied by coherence and Granger causality analysis. For this purpose, we made appropriate bipolar derivations of the LFPs around the alpha current generators to yield bipolar LFPs. The bipolar LFP signals were subjected to MVAR spectral analysis from which power, coherence, and Granger causality spectral estimates were derived. The essential steps of MVAR can be summarized as follows. Let p channels (p = 2 or 3 in this study) of bipolar LFP recordings at time t be denoted by X_t = (x_1t, x_2t, … x_pt)^T, where T stands for matrix transposition. Assume that the data over an analysis window are described by an MVAR model: where E_t is a temporally uncorrelated residual error series with covariance matrix Σ, and A_k are p × p coefficient matrices to be estimated from data (Ding et al., 2000). The MVAR model order m was determined by the Akaike information criterion (Akaike, 1974). For the data analyzed in this study, m = 10 (50 ms) was chosen as a tradeoff between sufficient spectral resolution and overparameterization. An adequate parametric model fit of the data means that the residual noise process E_t in Equation 1 must be temporally uncorrelated (white) (Lutkepohl, 2005). The Durbin–Watson test (Durbin and Watson, 1950) was used to check the goodness of fit. The whiteness of the residuals was confirmed at the p = 0.05 significance level.

Once the model coefficients A_k and Σ are estimated, the spectral matrix can be evaluated as where the asterisk denotes matrix transposition and complex conjugation, and H(f) = (Σ_k=0^mA_ke^−2πikf)⁻¹ is the transfer function. The power spectrum of channel l is given by S_ll(f), which is the lth diagonal element of the spectral matrix S(f). The coherence spectrum between channel l and channel k is as follows: The value of coherence ranges from 1 to 0, with 1 indicating maximum interdependence between channel l and channel k at frequency f and 0 indicating no interdependence. For p = 2 in Equation 1, the Granger causality spectrum from x_2t to x_1t is defined as follows (Geweke, 1982; Brovelli et al., 2004): which can be interpreted as the proportion of x_2t's causal contribution to the power of the x_1t series at frequency f. The logarithm is taken to preserve certain favorable statistical properties. Similarly, the causality spectrum from x_1t and x_2t can be obtained by switching the indices 1 and 2 in Equation 4. In addition to the bivariate causality algorithm, we also applied a conditional Granger causality algorithm (Chen et al., 2006; Ding et al., 2006), which allows the examination of whether the causal influence between two time series is direct or is mediated by a third time series.

Interpretation of Granger causality.

For two simultaneously measured time series, one series can be called causal to the other if we can better predict the second series by incorporating past knowledge of the first one (Wiener, 1956). This concept was later adopted and formalized by Granger (1969) in the context of linear regression models of stochastic processes (see Eq. 1). Specifically, if the variance of the prediction error for the second time series at the present time is reduced by including past measurements from the first time series in the linear regression model, then the first time series can be said to have a causal (directional or driving) influence on the second time series. Reversing the roles of the two time series, one repeats the process to address the question of causal influence in the opposite direction. From this definition, it is clear that the flow of time plays an essential role in allowing inferences to be made about directional causal influences from time series data. Alternatively, improvement in prediction can also be viewed from the perspective of comparing relative estimates of conditional probability. Recent work has modeled the relation between multiple point processes along this view (Okatan et al., 2005; Truccolo et al., 2005). In our analysis, the bipolar LFP data constitute the time series, and Granger causal influence is equated with the direction of synaptic transmission between neuronal ensembles. Physiologically, rhythmic synaptic activity (LFP) in a neuronal ensemble gives rise to periodic bursts of action potential firing (MUA), which, through synaptic transmission, drives the rhythmic synaptic activity in the target neuronal ensemble. Thus, evaluation of the causal relation between two bipolar LFP time series forms the basis for identifying the layer containing the pacemaking cells in a column. In this study, Granger causality analysis, as a principled approach for inferring causal influence among neural time series (Freiwald et al., 1999; Bernasconi et al., 2000; Liang et al., 2000; Baccalá and Sameshima, 2001; Hesse et al., 2003; Brovelli et al., 2004; Ding et al., 2006), is used in lieu of the trisection method in in vitro studies to identify the laminar pacemakers of the alpha rhythm. We analyzed bipolar LFPs rather then simple LFPs because the latter are vulnerable to contamination by volume conduction. The use of bipolar LFPs rather than CSDs was based on the consideration that the latter are more noisy in single trials (Lakatos et al., 2005, 2007; Shah et al., 2004), and excessive noise can lead to spurious causality results (Nalatore et al., 2007).

Assessment of statistical significance.

For coherence and Granger causality spectra, we have adopted a random permutation approach (Brovelli et al., 2004) to build a baseline null-hypothesis distribution from which statistical significance can be derived. Consider two channels of recordings with many repeated trials (one trial being a 200 ms epoch for this work). We can reasonably assume that the data from different trials are approximately independent of one another. Randomly pairing data for channel 1 with data for channel 2 from a different trial leads to the creation of a synthetic ensemble of trials for which there is no interdependence between the two channels by construction. However, the temporal structure within a channel is preserved. Performing such random pairing with many different permutations will result in a distribution of coherence or causality spectra corresponding to the null hypothesis of no statistical interdependence. The calculated value for a given statistic from the actual data is compared with this baseline null-hypothesis distribution for the assessment of significance levels. In addition, the difference between two nonzero quantities was tested by a bootstrap resampling technique (Efron, 1982).

Visual alpha power and auditory reaction time.

The functional significance of alpha activity in visual cortices was assessed by correlating its magnitude with auditory behavioral performance. For each 30 s segment of contiguous data analyzed (five to seven such segments in a penetration), the monkey made on average six lever-pressing responses to the deviant auditory stimulus (see Experiment, Paradigm). With a target detection rate of ∼90%, there were 25–35 correct responses for each penetration. For each correct response, an 800 ms interval before the onset of a deviant stimulus was defined as the prestimulus interval. All the correct responses in a penetration were rank ordered by the reaction time and then sorted into groups of five, starting with the fastest reaction times and proceeding to the slowest, each group sharing four trials with the previous one. For each group, the alpha power was estimated for each electrode contact and averaged across all 14 electrode contacts. The relation between the average group reaction time and average group alpha power was assessed by Spearman's rank correlation. The same alpha power and reaction time relation was also examined on a contact-by-contact basis.