Abstract
Electrocorticography (ECoG) methodologically bridges basic neuroscience and understanding of human brains in health and disease. However, the localization of ECoG signals across the surface of the brain and the spatial distribution of their generating neuronal sources are poorly understood. To address this gap, we recorded from rat auditory cortex using customized μECoG, and simulated cortical surface electrical potentials with a full-scale, biophysically detailed cortical column model. Experimentally, μECoG-derived auditory representations were tonotopically organized and signals were anisotropically localized to less than or equal to ±200 μm, that is, a single cortical column. Biophysical simulations reproduce experimental findings and indicate that neurons in cortical layers V and VI contribute ∼85% of evoked high-gamma signal recorded at the surface. Cell number and synchrony were the primary biophysical properties determining laminar contributions to evoked μECoG signals, whereas distance was only a minimal factor. Thus, evoked μECoG signals primarily originate from neurons in the infragranular layers of a single cortical column.
SIGNIFICANCE STATEMENT ECoG methodologically bridges basic neuroscience and understanding of human brains in health and disease. However, the localization of ECoG signals across the surface of the brain and the spatial distribution of their generating neuronal sources are poorly understood. We investigated the localization and origins of sensory-evoked ECoG responses. We experimentally found that ECoG responses were anisotropically localized to a cortical column. Biophysically detailed simulations revealed that neurons in layers V and VI were the primary sources of evoked ECoG responses. These results indicate that evoked ECoG high-gamma responses are primarily generated by the population spike rate of pyramidal neurons in layers V and VI of single cortical columns and highlight the possibility of understanding how microscopic sources produce mesoscale signals.
Introduction
Brains are composed of neuronal microcircuits (e.g., cortical columns) that perform specific computations that are integrated into larger networks (Koch, 1994; Shepherd, 2004). Microscale measurements investigating the activity of individual neurons and small neuronal populations have yielded insight into microcircuit mechanisms of local computations. Macroscale measurements (e.g., fMRI) have revealed principles of global processing across the brain (Buzsáki et al., 2012). Much less is known about how local neural processing is organized and coordinated across distributed brain networks. Measuring distributed cortical function is critical to understanding complex perceptions and behaviors (Canolty et al., 2006; Bouchard et al., 2013; Khodagholy, et al., 2017). High-density micro-electrocorticography (μECoG) arrays record neuronal activity directly from the cortical surface, can have tight spacing of many thousands of channels, and are minimally invasive (Ledochowitsch et al., 2011; Viventi et al., 2011; Khodagholy et al., 2015). Furthermore, because ECoG/μECoG is used in humans, it is a critical methodological bridge between basic neuroscience findings and understanding human brains in health and disease (Khodagholy et al., 2015). However, use of (μ)ECoG for basic neuroscience is impeded by a lack of understanding of the spatial localization of recorded signals across the surface and the specific neuronal sources generating those signals.
ECoG records cortical surface electrical potentials (CSEPs), which reflect a weighted superposition of all electrical sources in the brain (Buzsáki et al., 2012; Einevoll et al., 2013). Many ECoG studies focus on lower frequencies (e.g., <60 Hz), although it is becoming common to use the high-gamma (Hγ; 65–170 Hz; Canolty et al., 2006; Bouchard and Chang, 2014) band. The motivation for using Hγ stems from the proposal that higher frequencies reflect more spatially localized signals (Buzsáki et al., 2012). However, experimental estimates of spatial spread of correlations range from several hundred micrometers to a few millimeters. Determining the spatial localization of evoked ECoG signals is critical for both interpreting data as well as guiding the design of future devices.
During sensory processing, neuronal populations are activated with spatial correlations that partly depend on cortical functional organization and the biophysics of signal propagation through cortical tissue (Lindén et al., 2011; Lin et al., 2015). In primary sensory cortices [e.g., primary auditory cortex (A1)], the functional organization of response properties varies smoothly across adjacent columns (Schreiner and Winer, 2007). Cortical columns in rodents have a radius of ∼350 μm and span ∼1800–2100 μm in depth, including all six cortical layers, which are composed of different neuron types (Mountcastle et al., 1997; Tischbirek et al., 2019). Although different excitatory and inhibitory cell types within columnar microcircuits play different roles in sensory computation, in general, neurons within a column share similar sensory tuning properties (Atencio and Schreiner, 2010a, 2012). At the same time, electrical fields spread passively through cortex, diffusing the signal (Lindén et al., 2011). Thus, both passive spread through the tissue and functional organization of neuronal populations are potential mechanisms dictating localization of evoked ECoG signals across the surface. However, which of the two is the dominant mechanism is unknown.
We lack biophysical understanding of the sources generating CSEPs. The spatial spread of population local field potential (LFP) is determined by single-neuron morphology, temporal correlations between sources (synchrony), the number of sources, and the distance of the sources to the electrode (Lindén et al., 2011; Łęski et al., 2013). Unlike for intracortical LFPs, full-scale biophysically detailed simulations have not been used to understand CSEPs (Miller et al., 2009; Hill et al., 2018). The layers of a cortical column are composed of different numbers of neurons with distinct morphologies and synchron, and are by definition at different distances from the surface electrode (Markram et al., 2015). Thus, a full-scale, biophysically detailed cortical column model is required to determine the precise laminar and cellular origins of evoked CSEPs.
Given the large number of cells in infragranular layers, we hypothesized that evoked high gamma is spatially localized to a cortical column and primarily generated by action potentials from neurons in layers V and VI. To test this hypothesis, we combined direct experimentation with biophysical modeling.
Materials and Methods
Experimental model and subject details
Data from four rats (female Sprague Dawley) were used in this study. All animal procedures were performed in accordance with established animal care protocols approved by the Lawrence Berkeley National Laboratory, Institutional Animal Care and Use Committees.
Method details
Electrophysiological recordings
All neural data were recorded with a multichannel amplifier optically connected to a digital signal processor (Tucker-Davis Technologies). Signals were acquired at 12 kHz and low-pass filtered to the Nyquist frequency (6 kHz).
Rodent Preparations
We performed experiments in four anesthetized female Sprague Dawley rats. Animals were given a 1 mg/kg subcutaneous injection of Dexamethasone the night before a procedure to reduce cerebral edema. An anesthetic state was induced with an inductive dose of ketamine (95 mg/kg, i.p.) and xylazine (10 mg/kg, i.p.). Anesthetic state was assessed using toe pinch reflex and monitoring respiration rate. Additional doses of ketamine (55 mg/kg, i.p.) and xylazine (5 mg/kg, i.p.) were administered as needed to maintain a negative reflex and a regular reduced respiration rate. Respiration was supported with a preoperative subcutaneous injection of atropine (0.2 mg/kg) and a perioperative nose cone supplying 0.8L/min of O2. A water heating bed provided thermostatic regulation. To prevent dehydration over the 10 h surgery and recording session, subcutaneous saline injections (1 ml/kg) were provided every 3 h. Once anesthetized, the rodent was affixed to a snout stereotax without earbars.
After stable anesthetic state was achieved, an incision was made along the sagittal midline. All the soft tissue on top of the skull was removed to reveal the lambda and bregma fissures. Two 1 mm burr holes were drilled over nonauditory cortical areas–one between the lambda and bregma on the left hemisphere and another anterior to bregma on the right hemisphere. These serve to reduce intracranial pressure and provide a reference for electrophysiological recordings. The right masseter muscle was then transected to uncover the portion of cranium lying over the right auditory cortex. Using a 1 mm diamond tapered round Stryker dental drill, a craniotomy was performed to expose the cortex. Given the flexible grids we used, anatomic coordinates relative to a skull landmark are not possible. Auditory cortex was identified electrophysiologically by its anteroposterior tonotopic organization (Schreiner and Winer, 2007).
Recording devices
Rodent electrophysiological recordings were made with custom designed 128-channel μECoG grids (initially fabricated in house, purchased at the time through Cortera Neurotechnologies). These arrays were placed over the primary auditory cortex as identified by anatomic and physiological properties. For the rat experiments, the dura was surgically removed, and the μECoG grid was placed directly on the pial surface. The μECoG grids were placed over the primary auditory cortex and grounded via a silver wire inserted into a nonauditory cortical area in the contralateral hemisphere, which also served as the reference. Each contact on the grid had an impedance of 30 ± 10 kΩ after electroplating with platinum black (measured in 1× PBS at 1 kHz) and had an exposed diameter of 40 μm, with a 200 μm interelectrode pitch.
Auditory stimuli
Pure tone pips (50 ms in duration) were played through a DVD player (Sony), attached to an amplifier (Tucker-Davis Technologies) and played through an electrostatic speaker located near the left ear of the animal (80 dB SPL, A-scale). The frequency and amplitude of the tone pips were parametrically varied. The frequencies spanned six octaves (500 Hz to 32 kHz) in 30 increments and eight attenuations from 0 to −70 dB (0 dB attenuation at 80 dB SPL). Each frequency-attenuation combination was presented 25 times in pseudorandom order with 250 ms of interstimulus silence.
Data analysis
All analysis was performed using code written in MATLAB (MathWorks) or Python.
Spectral analysis of CSEPs
We calculated the spectrogram of the entire recorded electrical potential time series for each electrode from 4 to 1200 Hz (54 bins) using a constant-Q wavelet transform (Schorkhuber and Klapuri, 2010). Constant-Q refers to a time-frequency decomposition in which frequency bins are geometrically spaced, and Q-factors (ratios of the center frequencies to bandwidths) are equal for all bins. The noncausal component of displayed responses (Fig. 1d) is because of the large bandwidth at lower frequencies of our constant-Q time-frequency transform.
The amplitude for each frequency bin in the neural signal was normalized relative to baseline statistics by z-scoring. The baseline statistics were computed from the periods of silence between each stimulus presentation (The 50 ms immediately following each stimulus presentation is excluded from the baseline period.). Z-scoring largely removes the canonical ∼1/fα falloff of power with frequency (characteristic of many natural signals), highlighting stimulus-evoked changes. We note that this procedure is preferred to examining the residuals from a fit and subtract methodology, given the potential issue of fitting power laws especially at the extremes of the frequency range.
Responses to stimuli were taken as the average z-scored activity ±5 ms around the peak response time after the onset of the auditory stimulus. Peak response time was calculated using the high-gamma (65–170 Hz) component of the signal, which was computed as the average of z-scores for frequency bands whose center frequency falls in that range. The shielding and grounding of our rodent experimental recording systems was sufficient to avoid any significant 60 Hz noise in the μECoG recordings. We determined z-scores for the canonical neural frequency bands by taking the mean z-scores across frequency bins in the corresponding frequency range; Beta (10–27 Hz), Gamma (30–57 Hz), High Gamma (65-170 Hz), Ultra-High Gamma (180–450 Hz), and the Multiunit Activity Range (500–1100 Hz).
Analysis of responses to pure tones
The frequency response area (FRA) depicts the modulation of a CSEP components response as a function of frequencies (F) and amplitudes (A). For each frequency-amplitude pair (f, a) in the stimulus set (see above, Auditory stimuli), we identified the response peak,
We estimated an FRA boundary to define a set of frequency-amplitude pairs that evoked a response. Intuitively, this is the stimulus-driven portion of the FRA and corresponds with the part of the FRA that resides within the canonical V shape. These boundaries were extracted using an approach identical to one used in a previous study characterizing rodent auditory responses (Guo et al., 2012). A response-frequency function,
The
To fit the FRA boundary,
Portions of the FRA boundary that fell outside of the FRA were removed, resulting in the final FRA boundary. This method yielded FRA boundaries that closely resemble those curated manually.
We determined whether a recording site was tuned using a permutation test. For a given FRA, we once again define
A tuned response exhibits a larger SD across frequencies than an untuned (flat)
We used the FRA boundary to further characterize the response properties of each recording site. We defined a best frequency (BF) as the weighted average of the FRA boundary, and calculated bandwidth (BW) as the full width of FRA boundary at half-maximum,
Spatial analysis of CSEPs
We assessed the spatial spread of CSEPs by fitting a general linear model. For a given CSEP component, the responses at an individual single electrode (Ri) was modeled as a linear combination of responses at all other active electrodes (Rj ≠ i), corrupted by independent and identically distributed Gaussian noise (ε) as follows:
The UoILasso algorithm
We used a novel statistical inference procedure (UoILasso; Bouchard et al., 2017) to fit the general linear model in Equation 7. Although a detailed description of this algorithm is outside the scope of this manuscript, here we provide the motivation, outline the innovations, and summarize the main statistical result of these innovations. The interested reader is encouraged to see Bouchard et al. (2017) for further details.
Generally speaking, in regression and classification, it is common to employ sparsity-inducing regularization to attempt to achieve simultaneously the following two related but quite different goals: to identify the features important for prediction (i.e., model selection) and to estimate the associated model parameters (i.e., model estimation). For example, the Lasso algorithm in linear regression uses L1-regularization to penalize the total magnitude of model parameters (
It is well known that this type of regularization implies a prior assumption about the distribution of the parameter (e.g., L1-regularization implicitly assumes a Laplacian prior distribution). However, strong sparsity-inducing regularization (i.e., large values of λ), which is common when there are many more potential features (p) than data samples (n; i.e., the so-called small n/p regime) can severely hinder the interpretation of model parameters. For example, although sparsity may be achieved, incorrect features may be chosen, and parameters estimates may be biased. In addition, it can impede model selection and estimation when the true model distribution deviates from the assumed distribution.
To overcome these and other issues, we have recently introduced a novel statistical-machine-learning framework called Union of Intersections (UoI; Bouchard et al., 2017). Methods based on UoI perform model selection and model estimation through intersection and union operations, respectively, leading to enhanced model selection and estimation. Focusing on linear regression, the UoILasso algorithm has three central innovations: (1) calculates model supports (Sl) using an intersection operation over bootstrap resamples for a range of regularization parameters λ (increases in λ shrink all values of β toward 0), efficiently constructing a family of potential model supports
Biophysically detailed simulation of cortical column
Using the NEURON simulation environment running on 64 nodes of Cori (Cray XC70) at the National Energy Research Scientific Computing Center, we implemented and successfully executed a biophysically detailed model of a cortical column and the CSEP produced by the activity of this neuronal network. Our simulation is a compartmental model, in which the electrical activity of one or more neurons is simulated by modeling the neuron(s) as a series of small cylindrical segments over which the membrane potential and currents can be taken as approximately constant. The evolution of the membrane potential is given by the following cable equation:
We constructed a compartmental model of the neurons in one column of rat sensory cortex based on publicly available data from the Blue Brain Project (BBP; Markram et al., 2015). These data included the spatial location and connectivity matrix of all cells in the column, as well as tuned and experimentally validated models of individual neurons from all cortical laminae including their electrical characteristics (ion channel models and associated parameters such as spatially varying ionic conductances, membrane resistance/capacitance, etc.) and their detailed morphologies based on full reconstructions of neurons observed experimentally (or algorithmically generated clones thereof), reflecting the full diversity of neurons known to be found in rodent somatosensory cortex. To save computation time, we instantiate 80% of the cells in the model, selected at random. We confirmed that there was no difference in the spectrum between 80% and 100%. The BBP dataset provides five reconstructed (or cloned) morphologies for each of 207 distinct cell types (Markram et al., 2015). Each neuron in our model is represented by one of the five morphologies for the cell type of that neuron, chosen at random and rotated by a random angle about a line passing through the soma of the cell parallel to the longitudinal axis of the column.
Neurons in our simulated column are innervated by synapses from the following three populations: 1)
Synapses from all populations produce membrane currents
Parameters of model synapses
We applied a modest amount of hand-tuning of these parameters to achieve a reasonable baseline firing rate (3–10 Hz) during time periods when the thalamocortical spike trains fire at
CSEP of the simulated column
The line source approximation (LSA) is used to simulate the extracellular potential at the cortical surface because of the transmembrane currents in each segment of each neuron in the simulation, assuming an isotropic and purely Ohmic extracellular medium (Logothetis et al., 2007). These contributions from individual neuronal segments are summed to compute the total CSEP because of the entire simulated column as follows:
Dependence of response magnitude and frequency on input amplitude
We determined the relationship between z-scored response magnitude (between 10 and 200 Hz) and the frequency of CSEPs (between 10 and 200 Hz) by varying the magnitude of inputs in both the experimental data and in the simulations. For simulation data, we varied the frequency
Laminar lesions and isolations
By summing only neuronal segments belonging to neurons in an individual cortical layer, we obtain the in silico contributions
Similarly, we perform the superposition of LSA-computed potentials excluding those segments belonging to neurons in a particular layer in the following:
Finally, by summing only the neuronal segments located within a particular range of depths below the surface, we obtain the in silico contributions to the CSEP from 200 μm slices of the column as follows:
Normalization of simulated CSEP components to baseline
To assess the contributions of sources at different depths to CSEPs (see Figs. 6, 7), we would like to normalize the CSEP contributions to baseline in a way that preserves their relative magnitudes. Using the z-score of each contribution to its own baseline does not preserve relative magnitudes between contributions. For example, two contributions that differ only by a constant multiplicative factor will have the same z-score. That is, if
As with the z-score analysis used for the full CSEP, this normalization is done independently for each neural frequency bin. In separate analyses, we confirmed that our conclusions are robust to the choice of normalization. This normalization was also done for the analysis (see Fig. 4).
Dependence of simulated CSEP contributions on number of segments/neurons, depth, and synchrony
The CSEP contribution from a given subset of the column (anatomically defined layer, or 200 μm slice) will depend on the number of neuronal segments in the subset, the distance of the neurons from the recording electrodes, and the synchrony of the activity of those neurons. To determine the relative importance of these three factors, we performed an L2-regularized regression to fit the magnitude of the high-gamma peak (the maximum of each normalized CSEP contribution across frequency bins) of each contribution as a linear function of the number of simulated neurons in each layer, or the number of neuronal segments in each 200 μm slice, the average depth of the contributing segments below the surface, and the synchrony between somatic membrane potentials, defined as the average of Pearson's correlation coefficient of the membrane potentials over all pairs of somas in the subset. Each of the three independent variables, and the dependent variable (high-gamma contribution magnitude) was normalized by dividing by the maximum across layers or slices. The L2 regularization parameter was chosen to be
Data availability
The datasets are publicly available from https://crcns.org/data-sets/ac/ac-7/about-ac-5, and the software is publicly available from https://github.com/BouchardLab/OriginsOfECoG. Further information and requests for resources should be directed to and will be fulfilled by Kristofer E. Bouchard (kebouchard{at}lbl.gov).
Results
We designed and fabricated lithographically defined, high-resolution μECoG arrays with customizable electrode geometry and spacing (Rubehn et al., 2009; Ledochowitsch et al., 2011; Viventi et al., 2011). Each contact on the μECoG arrays had a diameter of 40 μm (approximately the size of a cortical mini column), and an impedance of 30 ± 10 kΩ (measured in 1× PBS at 1 kHz; Ledochowitsch et al., 2011). To ensure that we fully captured the range of functional time scales associated with the diversity of neurobiological signals, we recorded wide band (2–12,000 Hz) electrophysiological data.
Stimulus-evoked cortical surface electrical potentials exhibit large peaks in the high-gamma range
An example of a 128-channel μECoG grid is shown on the cortical surface (subdural placement) of an anesthetized rat in the photomicrograph in Figure 1a. In this preparation, we recorded large amplitude, fast CSEPs in response to the presentation of auditory tone pips (see above, Materials and Methods). We played back stimuli consisting of short (50 ms) pure tone pips with varying frequency and intensity (amplitude) at five recording locations in four rats. An example recording from a single electrode during four consecutive stimulus presentations is shown in Figure 1b. Here, the top depicts stimulus amplitude and frequency, the middle displays the (normalized) neural spectrogram of the full response, and the bottom shows the amplitude of the high-gamma component of the neural response. For each recording electrode and CSEP frequency component, we normalized (z-scored, relative to baseline statistics taken during the interstimulus interval) the time-varying amplitude for each CSEP frequency separately, removed untuned electrodes (those that exhibit no stimulus frequency selectivity), and computed the best frequency for each tuned electrode (see above, Materials and Methods). The best frequency is the stimulus frequency, which maximally drives activity on that electrode. Figure 1, c and d, shows the average evoked potential (c) and neural spectrogram (d) derived from the recorded electrical potential in response to the best frequency of the electrode at a single amplitude (N = 25 trials). In both single trials (Fig. 1b), as well as the average (Fig. 1d), the best frequency of the electrode evoked large-amplitude, rapid CSEP deflections. The Hγ (65–170 Hz) component of the evoked response was observed to be the most robust, often exceeding 5 SDs of the baseline in response to the best frequency (e.g., Fig. 1b, second stimulation, BF). Frequencies below 10 Hz exhibited little to no response modulation (Fig. 1b,d); thus, we analyze frequencies between 10 and 1100 Hz, with a focus on high gamma.
Stimulus-evoked cortical surface electrical potentials exhibit large peaks in the high-gamma range. a, Photomicrograph of an 8 × 16 μECoG grid (pitch, 200 μm; contact diameter, 40 μm) on the surface of rat A1. b, Top, Tone stimulus played during experimental recordings. Middle, z-Scored spectral decomposition of single-trial evoked cortical surface electrical potentials from a single electrode. Bottom, High-gamma component of single-trial evoked cortical surface electrical potentials indicated by horizontal dashed lines (middle). c, Trial-averaged evoked cortical surface electrical potential on one μECoG electrode in response to presentations of the best tuned frequency of that electrode. d, Trial-averaged neural spectrogram for the electrode shown in c in response to presentations of its best tuned frequency. Dashed vertical lines in c and d represent stimulus onset and offset. Red vertical lines in c and d correspond to the time window of extracted evoked response used for subsequent analysis. e, Grand-average (mean ± SE) z-scored amplitude as a function of frequency across all tuned electrodes (N = 333).
To summarize the frequency content of evoked CSEPs, we averaged across presentations of the best stimulus at one amplitude. For each electrode, we extracted mean z-scored responses across all neural frequency components in a ±5 ms window around the time of the peak high-gamma response (Fig. 1c,d, red vertical lines). We included all electrodes with a tuned response in the high-gamma band (N = 333 electrodes from 5 μECoG placements on auditory cortex in four rats). Figure 1e plots the averaged (N = 333 electrodes, mean ± SE) z-scored response as a function of frequency. On average, we found that evoked responses were unimodally peaked around the Hγ band, with notable responses in the multiunit activity range (MuA; >500 Hz).
Robust frequency tuning and high-resolution tonotopic maps derived from μECoG
We next determined auditory receptive field properties and the spatial organization of evoked CSEPs. The plots of Figure 2, a and b, present FRA heat maps derived from Hγ activity (Fig. 2a) and multiunit activity (tMuA, after application of a threshold to the MuA band for event detection; see above, Materials and Methods; Fig. 2b) in response to these stimuli. In each FRA, pixels correspond to stimulus frequency-intensity pairings, and are colored according to the mean evoked z-scored signal (N = 25 stimuli for each), with blank spaces indicating untuned responses (see above, Materials and Methods). Data are displayed for several electrodes spanning 1.8 mm anteroposterior (AP) and 0.4 mm dorsoventral (DV) over rat A1. We found that the response profiles exhibited clear frequency tuning and the canonical V-shaped profiles expected of auditory cortical neurons (Polley et al., 2007; Guo et al., 2012; Escabí et al., 2014). Response profiles were largely similar between Hγ and tMuA. Additionally, there was a smooth gradation of best frequencies across the AP axis, with low frequencies posterior, high frequencies anterior, and similar best frequencies along the DV axis, suggestive of tonotopic organization.
Robust frequency tuning and high-resolution tonotopic maps from μECoG. a, b, FRA surfaces recorded from a μECoG array. Subplots correspond to responses of a single electrode and are organized according to electrode position on the grid/brain. In each subplot, pixels correspond to a stimulus frequency-intensity pairing and are colored according to the mean evoked z-score; Ηγ (a), tMuA (b). c, High-resolution tonotopic organization of multiple auditory cortical fields derived from Ηγ activity. Each pixel is color coded according to the best frequency of that electrode. The 8 × 16 μECoG array displayed here covered multiple auditory cortical fields (A1, PAF, and VAF) and the approximate boundaries are demarcated (black lines). d, Differential tuning at neighboring electrodes. FRAs are plotted for four electrodes (numbered as in c) and show that neighboring electrodes (1 vs 2; 3 vs 4) can have different response properties. e, f, Average normalized response surface for all electrodes with significantly tuned Ηγ (N = 333) and tMuA (N = 113) auditory responses. White line in each plot demarcates the FRA response boundaries. g, Across all tuned electrodes, the average (mean ± SE) FRA response boundaries for CSEP components (demarcated by colors) where similar. h, i, Distributions (25th–50th–75th percentiles) of best-frequencies (h) and bandwidths (i) for all tuned responses for CSEP components.
As Hγ activity had the largest number of tuned channels, we first visualized tonotopic organization by coloring each electrode (pixel) according to its best frequency extracted from the Hγ band (Fig. 2c). The 8 × 16 μECoG array displayed here covered multiple auditory cortical fields [A1, posterior auditory field (PAF), and ventral auditory field (VAF)], and the functionally defined boundaries are demarcated (Fig. 2c, black lines; Polley et al., 2007; Schreiner and Winer, 2007). Within a given auditory cortical field, there was a smooth gradation of best frequencies, with low frequencies posterior and high frequencies anterior. Tuning across the dorsoventral direction was largely similar within an auditory field. Interestingly, although frequency tuning generally varied smoothly as a function of distance between electrodes within an auditory area, we observed examples of different tuning at neighboring electrodes located in different auditory fields. For example, the FRAs for the electrodes demarcated 1 and 2 and 3 and 4 in Figure 2c are plotted in Figure 2d and show that neighboring electrodes (1 vs 2, 3 vs 4) can have different response properties, with 3 versus 4 being a particularly stark contrast. This suggests a high degree of spatial localization of the recorded signals. These results demonstrate the ability to resolve the tonotopic organization of multiple auditory cortical fields with very high resolution using high-frequency signals of CSEPs.
Figure 2a,b suggest that auditory responses of CSEP components from the same electrode are similar. The plots in Figure 2e and f, display average normalized FRAs across all channels with tuned responses in the Hγ and tMuA components, which were indeed similar. To quantify this, we first determined a boundary that separated the responsive portions of the FRA from the unresponsive portions (e.g., Fig. 2, e and f, white lines; see above, Materials and Methods). The average FRA boundaries for different frequency components across all tuned channels (β = 260, γ = 292, Hγ = 333, uHγ = 302, tMuA = 113) are displayed in Figure 2g (mean ± SE) and were highly overlapping. We found that median BFs (auditory, extracted from FRA) were ∼8.5 kHz (Fig. 2h), and there was a mild effect of CSEP component on BF (Kruskal–Wallis, df = 4, χ2 = 28, p < 0.001). The range of values in our dataset (interquartile ranges; Fig. 2h) makes the functional relevance of the marginal statistical significance for best frequency questionable. Indeed, the best frequencies extracted from different components at an electrode were highly correlated (R, median = 0.89; range = [0.8, 0.93], p < 10−5 for all). This indicates that different high-frequency CSEP components are generated by neurons with very similar tuning properties, perhaps from the same cortical column. Additionally, the median (auditory) BWs (extracted from FRA as full-width at half-maximum) were ∼1.5 octaves (Fig. 2i), and there was a robust effect of CSEP component on BW (Kruskal–Wallis, df = 4, χ2 = 116, p < 10−22). Thus, there was a systematic decrease in bandwidth with increasing CSEP component in the electrical potential (Fig. 2i), perhaps reflecting greater spatial spread of the lower frequency signals.
CSEPs are anisotropically localized to a cortical column
The spatial spread of electrical signals in the brain is of great interest, both for its importance in interpreting recorded electrical potentials and for its practical implications for sensor design. In principle, two primary factors that contribute to the spatial spread of correlations in electrophysiology recordings are diffusion because of the electrical properties of the tissue and the spatial organization of neuronal function (e.g., tonotopy). If diffusion is the main determinant, then the spatial spread is expected to be isotropic. However, given the tonotopic organization of rat auditory cortex (Fig. 2c), we hypothesized stronger correlations in the isotonic dimension (DV) versus the heterotonic dimension (AP). Furthermore, because tonotopy arises from the organization of cortical columns that contain neurons with similar tuning properties (Schreiner and Winer, 2007; Tischbirek et al., 2019), we hypothesized that the CSEP signals would be tightly localized to the approximate diameter of a column (200–500 μm; Mountcastle et al., 1997; Schreiner and Winer, 2007; Tischbirek et al., 2019).
To assess spatial spread, for each CSEP frequency component, we fit a general linear model to single-trial tone responses for each target electrode as a function of the responses at all other electrodes. Our analysis method, based on sparse general linear models, largely mitigates potential confounds because of pairwise chaining of local correlations by factoring out the covariance matrix of the regressors. We fit the linear model using the UoILasso algorithm (Bouchard et al., 2017; see above, Materials and Methods), which provides accurate estimates of parameter values (R2 > 0.9 for all fits). For each target electrode, the parameter values from the fit model (weightings on responses at other electrodes) were normalized to their maximum value across frequency to ease comparisons across frequencies and electrodes. Figure 3, a and b, displays the spatial distributions of median fit parameters as a function of location relative to the target electrode (demarcated as X) for all frequency tuned electrodes in the Hγ and tMuA components. For both of these signals, we found that model parameters were extremely localized in both AP and DV directions (values quickly go to zero), and exhibited a marked anisotropy, with larger values for DV than AP. Also, parameters were larger for Hγ than for tMuA (grayscale), indicating that less tMuA variance could be explained by surrounding electrodes and thus suggesting a more localized signal.
CSEPs are anisotropically localized to a cortical column. a, Spatial distribution of weights from a regularized linear model of Hγ responses during the tone stimuli as a function of the other electrodes on the grid. Locations are all relative to the electrode used as the dependent variable in linear regression. Values are median across all N = 333 tuned (in Hγ) electrodes. b, Spatial distribution of normalized weights for tMuA. Values are median across all N = 113 tuned (in tMuA) electrodes. c, Median ± SD of normalized linear weights across all electrodes as a function of distance in the AP (solid lines, left axis) and DV (dashed lines, right axis) dimensions along the grid. Different frequency bands are demarcated with colors. Note inverted orientation of distances along x-axis for AP (black) versus DV (gray).
We summarized the results of this analysis for the β through tMuA components (Fig. 3c, colors; N: β = 260, γ = 292, Hγ = 333, uHγ = 302, tMuA = 113) by plotting the median model parameters as a function of distance in the AP (solid lines, left y-axis, black distances on x-axis) and DV (dashed lines, right y-axis, gray distances on x-axis). Across frequency components, we found that ∼70% of the parameter magnitudes were concentrated at ±200 μm (4/143 grid locations with nonzero values, ∼3%), indicating that the vast majority of explanatory variation was localized immediately surrounding the electrode. Furthermore, the parameter values were significantly greater in the dorsoventral than the anteroposterior direction for all CSEP components (Wilcoxon signed-rank test, p < 10−4-10−28). Within both the AP and DV directions, a significant effect of CSEP frequency component on parameter magnitude at 200 μm was observed, with lower frequencies having larger values (Kruskal–Wallis; df = 4; DV, χ2 = 35, p < 10−6; AP, χ2 = 51, p < 10−12). This indicates that lower frequencies have a greater spatial spread and is in line with lower frequencies having broader tuning (Fig. 2). A cortical column is ∼350 μm in diameter (Mountcastle et al., 1997). Thus, these results indicate that cortical surface electrical potentials are localized to single cortical columns and that the degree of localization increases with increasing CSEP frequency. Furthermore, the localization is anisotropically distributed and aligned with tonotopic organization, indicating that differentiation of function across the cortical surface is the primary determinant of spatial correlations of evoked μECoG signals.
Biophysical in silico cortical column reproduces in vivo observed μECoG response
The experimental results presented above demonstrate that evoked CSEPs were tightly localized along the cortical surface, approximately to a single cortical column in the rat. This suggests that the sources of the μECoG signal are mostly located within the column directly underneath the electrode. To investigate the laminar and cellular origin of sources within a column that generate CSEPs, we simulated a full-scale, biophysically detailed model of a cortical column where the full morphology of each neuron is represented by hundreds to thousands of connected cylindrical neuronal segments (Markram et al., 2015). The column model and μECoG electrode is depicted in Figure 4a, where circles indicate the locations of (a subset of) neuronal somas (black, excitatory neurons; red, inhibitory neurons). Stimulus-evoked input to the column is provided by activating (with Poisson spike trains) thalamocortical synapses located throughout the column according to the distribution shown in Figure 4b, which also displays the cortical layers.
Biophysical in silico cortical column reproduces in vivo observed μECoG response. a, Rendering of a random subselection of 626 neurons in the simulated column (∼2% of the total). Black, excitatory neurons; red, inhibitory. Circles represent somas, lines represent dendritic structures. The position of the simulated μECoG electrode relative to the column is shown above. b, Distribution of synapses from the thalamus along the depth axis of the simulated cortical column. c, Data from one simulated stimulation and prestimulus/poststimulus silence. i, Population spike rate of thalamic and background cortical spike trains activating synapses in the column. ii, Spike raster of all neurons in the column versus soma depth (y-axis). Note that differences in raster density in part reflect differences in neuron density across cortical layers. iii, Population spiking (fraction of neurons spiking in 1 ms) of biophysically detailed cortical neurons. iv, Cell-averaged spike rate of biophysically detailed neurons in each layer. Darker shades indicate deeper layers. d, CSEP computed by the Line Source Approximation from all neurons in the column during a 150 ms window centered around the 50 ms tone pip stimulation. e, Spectrogram of the CSEP in d. Top, Mean normalization. Bottom, z-Score normalization. f, Frequency content of CSEP during 10 ms centered at the response peak (indicated with dotted red lines in d and e). Top, Mean normalization. Bottom, z-Score normalization. Individual electrode averages from experimental results are in gray, black is grand average. Individual stimulus presentations from simulations are in pink, red is grand average. All traces are normalized to their respective maxima. g, Whisker plots (median, IQR, 95% CI) of correlation coefficients comparing the frequency content of experimental results and average simulation results for z-score and mean normalizations.
An example of the activity of the column is displayed in Figure 4c. The biophysical neurons in the column received thalamic input in the form of Poisson spike trains that were modulated in time to emulate our tone stimulus (Fig. 4ci, black) and background Poisson spike trains (Fig. 4ci, gray) that were not modulated by the stimulus. In Figure 4cii we show the spike times (black, excitatory neurons; red, inhibitory neurons) in response to one presentation of the input stimulus (Fig. 4ci). Neurons are arranged by depth below the surface, which allows us to visualize the laminar boundaries as sharp changes in the density of firing reflecting different cell densities across layers. The fraction of neurons in the column firing action potentials is displayed in Figure 4ciii as a function of time. For most cortical layers, the time to peak of 15–20 ms (Fig. 4civ), as well as the following period of slightly elevated activity until stimulus offset, are both consistent with the in vivo recordings (Atencio and Schreiner, 2010b, 2012). Note that each curve in Figure 4civ is individually scaled, and the early timing of Layer I reflects the spiking of a small number of inhibitory neurons receiving convergent input from pyramidal cells in other layers. Layer I had essentially no contribution to the CSEP (see Fig. 6).
The biophysical model produces CSEPs (Fig. 4d–g) consistent with the high-frequency transient onset response observed in vivo. We computed the electrical potential at the cortical surface of the simulated column using the line source approximation (Holt and Koch, 1999) and processed the simulated data identically to the experimental data. Average (mean ± SE, N = 60 stimulus presentations) raw evoked cortical surface electrical potential from the model is plotted in Figure 4d. The average evoked spectrogram for simulated CSEPs is shown in Figure 4e (dashed black lines, stimulus; dashed red line, 10 ms window around the peak high-gamma response), and the extracted response as a function of frequency for the simulated CSEP is shown in Figure 4f (red), as well as the experimental data (black). We first examined the degree to which the simulation recreated the first moments (i.e., means) of the frequency content of the experimental data by plotting the mean normalized neural spectrogram (Fig. 4e,f, top). There is striking agreement in the mean frequency content of CSEPs collected experimentally and CSEPs generated by the simulations (Fig. 4f, top). Both experimental and simulated CSEPs exhibit a peak frequency of ∼100 Hz, and the spread of the signal around the peaks are highly overlapping. We next examined the degree to which the simulation recreated the second-order moment (i.e., variances) of the experimental data by examining the z-scored frequency content (Fig. 4e,f, bottom). This revealed that although the high-frequency structure was well preserved, the simulation had reduced z-scored responses in the lower frequency ranges compared with the experimental data. We quantified the similarity between experimental and simulated frequency content as the correlation coefficient between the average simulated response (Fig. 4f, dark red lines) and the average response at individual electrodes (Fig. 4f, light gray traces). The distributions of correlation coefficients for both mean and z-scored normalized responses are plotted in Figure 4g (mean normalization, median R = 0.91, (IQR) = [0.85, 0.93]; z-score normalization, median R = 0.70, IQR = [0.60, 0.83]). Thus, biophysical simulations accurately recreate key aspects of experimentally acquired data, indicating that they are a good forward model of CSEP generation.
In silico cortical column predicts experimentally observed relationship between response magnitude and frequency
The model makes testable predictions regarding the relationship between the magnitude and frequency content of CSEP responses. In the simulations, we varied the magnitude of the excitatory input to the network by increasing the mean firing rate of the thalamic spike trains during the stimulus. Analogously, in the experiments, we monitored the evoked CSEP in response to varying sound amplitudes at the best frequency of each electrode.
For the simulations, Figure 5a displays the average z-scored evoked response to stimulation of different amplitudes as functions of frequency (input magnitude given by color saturation, indicated by color bar, inset). We observed that the magnitude of CSEP response depended monotonically on input magnitude. More interestingly, we found that as the magnitude of the response increased, so did the frequency content of that response. This can be seen as a sweep toward the upper right of the individual traces as input magnitude increases. We quantified the relationship between response magnitude (z-scored response between 10 and 200 Hz) and the frequency content (frequency at peak response between 10 and 200 Hz). The pink-to-red squares in Figure 5c display the normalized maximum response magnitude versus normalized frequency at maximum response for varying input amplitudes (color saturation demarcates magnitude of input, circled square demarcates input used in Figs. 4, 6, 7). Intuitively, these effects were mediated by an increase in the population mean firing rate and spike synchrony resulting from increased input spike rate.
In silico cortical column predicts experimentally observed relationship between response magnitude and frequency. a, Average z-score as a function of frequency in eight simulations with variable input amplitude. b, Average z-score as a function of frequency in the experimental data for six different stimulus amplitudes. c, Normalized response magnitude versus normalized response frequency for experimental data (black, mean ± SD) and for simulations (red). Each data point corresponds to the response frequency and magnitude associated with a distinct input magnitude (response magnitude increases monotonically with input magnitude). Circled point indicates the input magnitude used in Figures 4, 6, 7. Orange dashed line is unity.
Next, we sought to determine whether this relationship between magnitude and frequency existed in the experimental data. Figure 5b displays the z-scored CSEP at an example electrode as a function of frequency in response to the BF stimulus presented at different amplitudes (inset, color bar). Similar to Figure 5a, we observed a sweep toward the upper right of the individual traces with increasing input amplitude (Fig. 5b). For frequency tuned electrodes, we calculated the same quantities (maximum response magnitude and frequency at maximum response) as a function of the amplitude of auditory input at the best frequency of the electrode in the tone stimuli (Fig. 5c, gray-to-black circles, mean ± SD, N ∈ [206, 299]). As in the simulations, we observed that increasing the input magnitude resulted in an increase in both the magnitude of the peak response and frequency at the peak response. Further, there is a striking correspondence in the curvature of response frequency versus response magnitude plots derived from experimental and simulation data (Fig. 5c). These results demonstrate a prediction made by the model that was confirmed by a novel experimental finding.
Evoked μECoG responses originate in infragranular layers
We next used the model to understand the spatial distribution of the generating sources of the CSEP. A key feature of the biophysical model is that the CSEP calculation is separate from the numerical simulation of the neurons in the column, enabling us to calculate CSEPs from arbitrary samples of neuronal segments in the column without perturbing the activity at all. We first examined the contributions to the CSEP from cortical layers by computing contribution of each layer to the CSEP individually (see above, Materials and Methods).
Figure 6a plots the raw evoked CSEP (inset, scale bar) as a function of time for each layer and indicates the average depth of neuronal somas for the layers. Surprisingly, we found that layers V and VI produce the largest evoked potentials, despite being the farthest away, whereas neurons in superficial layers contribute very little to the total CSEP. Figure 6b shows the frequency content of each contribution during a 10 ms window surrounding the response peak. The inset shows the relative magnitudes of the layer contributions in the band centered at 94 Hz, the apex of the high-gamma peak, shown as a dotted vertical line. As with the raw evoked potential, we found that infragranular layers also contribute most to the high-gamma component of ECoG responses, 51% from layer V, 35% from layer VI, and the remaining 14% coming from layers I and IV. To further probe the dependence of the high-gamma component on individual layers, we also performed lesion studies examining the CSEP produced when activity in one cortical layer is excluded, which confirmed that layer V lesions caused the largest reduction in high-gamma (data not shown).
Evoked μECoG responses originate in infragranular layers. a, Contributions to the simulated CSEP from anatomic layers. Top to bottom, Cortical layers I through VI. The sum of these contributions is the total CSEP. b, Frequency content of the laminar contributions during stimulus peak. Layer V and VI contributions dominate the high-gamma peak. c, Magnitude at peak frequency of the CSEP contribution of each cortical layer versus number of neurons in the layer. d, Magnitude at peak frequency of the CSEP contribution of each cortical layer versus average distance of cell bodies in the layer from the recording electrode. e, Magnitude at peak frequency of the CSEP contribution of each cortical layer versus synchrony of somatic membrane potentials averaged over all pairs of neurons in the layer. f, Pie chart showing the relative importance of these three factors in a linear model of the high-gamma peak contribution magnitudes of anatomic layers.
The results above appear counterintuitive when the contribution of sources is viewed only as a function of distance. However, in addition to distance, the number of sources and their correlations are additional biophysical factors that dictate the contribution of neuronal populations to a distally recorded signal (Lindén et al., 2011). A priori, the relative importance of these factors to determining laminar contributions to evoked CSEPs in a full-scale cortical column model is not clear. Thus, we plotted the peak high-gamma responses of each layer as a function of the number of simulated neurons in (Fig. 6c), the average distance of somas in each layer from the recording electrode (Fig. 6d), and the synchrony between somatic membrane potentials in each layer during the stimulus (Fig. 6e). We note that in Figure 6d the peak high-gamma response versus depth shows a positive slope, contrary to the physical principle that individual neurons farther from the electrode will contribute less to the signal. However, as is evident from these plots, there are correlations between depth and the other variables. For example, deeper layers tend to contain more neurons. Thus, we fit a regularized linear model to predict peak high-gamma magnitude across layers as a function of depth, number of segments, and synchrony of neuronal somas simultaneously, which fit the data well (R2 = 0.98). The relative magnitudes of the fit coefficients are plotted in Figure 6f, which shows that the number of segments and between-cell synchrony are the dominant factors that determine source contributions to CSEPs, whereas depth was a minor factor. To determine how robust these results were to baseline normalization, we performed the same analysis with a different normalization procedure and found very similar results (data not shown). Thus, infragranular layers contribute ∼86% of evoked CSEP responses because of their increased number of neurons and increased synchrony.
Evoked μECoG responses originate in sources 800–1400 μm below the surface
The previous results indicate that layers V and VI are the dominant sources to evoked CSEPs. However, because of the large, extended morphology of some neurons relative to the column depth, knowledge of the largest contributing anatomic layers does not necessarily imply precise knowledge of the spatial distribution of segments generating CSEPs. For example, the apical tufts of many layer V pyramidal neurons reach into layer I. Thus, we next isolate contributions to the CSEP from 200 μm axial slices of the column.
Most slices contain segments from neurons in more than one layer, and a given neuron can contribute to more than one slice. The breakdown of segments in each slice by anatomic layer is shown in Figure 7a, where each color represents one slice. For each slice, five bars are shown displaying the number of segments in that slice belonging to neurons in the five cortical layers. For example, the top slice is dominated by segments from layer V neurons (Fig. 7a, fourth column). The total number of neuronal segments in each slice is shown in Figure 7b, which makes clear that the slices between 800 and 1200 μm have the most segments. Figure 7c shows CSEPs calculated only from segments in the slices as a function of depth (inset, CSEP scale bar). The largest contributors to the evoked responses are the slices located from 800 to 1400 μm below the surface, that is, somas in layer V. We extracted a 10 ms window around the peak of the CSEP response and analyzed the frequency content of the contribution of each slice within that window. The results are shown in Figure 7d. The inset shows the relative magnitudes of the contributions of the slices at 94 Hz, the apex of the high-gamma peak, shown as a dotted vertical line. Here, we see that the slices spanning 800–1400 μm are also the ones contributing most to the high-gamma peak (56% total), which is where layer V somas are located. Thus, this analysis demonstrates that layer V somas are the major generating source of evoked ECoG signals.
Evoked μECoG responses originate in sources 800–1400 μm below the surface. a, Proportional breakdown of segments by anatomic layer. Most slices contain segments from neurons in multiple cortical layers. Bars represent proportion of total segments in the slice, different slices not to scale. b, Total number of simulated neuronal segments in each 200 μm axial slice of the column. c, Contributions to the CSEP from 200 μm slices, organized by depth. Top, Cortical surface. The sum of these contributions is the total CSEP shown in Figure 4b. d, Frequency content of the slice contributions during stimulus peak, colored by slice depth. Slices containing somas of layer V neurons dominate the high-gamma peak. e, Magnitude at peak frequency of the CSEP contribution of each slice versus number of neuronal segments in the slice. f, Magnitude at peak frequency of the CSEP contribution of each slice versus average distance of segments in the slice from the recording electrode. g, Magnitude at peak frequency of the CSEP contribution of each slice versus average synchrony in the slice. h, Pie chart showing the relative importance of the three factors in our linear model of the high-gamma peak contribution magnitudes of the slices.
As with the layer contributions, we sought to ascertain the relative importance of the number of segments in the slice (Fig. 7e), the depth of the slice below the surface (Fig. 7f), and the synchrony of membrane potentials of segments within the slice (Fig. 7g) in determining the high-gamma peak contribution magnitude. The results of a regularized linear regression predicting high-gamma peak from those parameters (R2 = 0.91) are shown in Figure 7h. As with the layer contributions, very similar results were observed using a different normalization procedure (data not shown). Similar to the layer contributions, we find that the number of segments and the synchrony are the most important factors determining the magnitude of the contribution of a slice to the CSEP.
Discussion
We found that evoked CSEPs had strongly nonmonotonic frequency structure, with a large peak at the 65–170 Hz range (Hγ). In experimental data, we quantitively demonstrated that evoked CSEPs were localized to a cortical column. Full-scale biophysical simulations of a cortical column reproduce experimentally observed evoked spectrum and indicate that evoked μECoG high-gamma is generated by infragranular neurons.
The spatial spread of electrical potentials is of long-standing debate and great interest for basic neuroscience with applications to sensor design and brain-machine interfaces. We demonstrate tight, anisotropic surface localization of CSEPs to less than or equal to ±200 μm, with higher frequencies exhibiting greater spatial localization. We observed examples of neighboring electrodes (in different auditory areas) with very different tuning. A column is ∼300 μm in diameter (Mountcastle et al., 1997). Thus, our results demonstrate that CSEPs are localizable to individual cortical columns.
Our results indicate a much more localized signal than has been directly quantified in previous studies (Escabí et al., 2014; Muller et al., 2016; Hermiz et al., 2018; Dubey and Ray, 2019), although there have been qualitative descriptions of tight localization (Tchoe et al., 2022). Most studies that directly quantify evoked signal spread calculate the pairwise correlation between electrodes and examine the decay with distance (Escabí et al., 2014; Muller et al., 2016). However, such analyses potentially confound long-range correlations because of direct interactions with correlations because of chaining of pairwise interactions. The sparse general linear models used here largely mitigate such confounds. Our results corroborate and extend a study of the spatial spread of low-frequency electrical potentials inside the cortex (Katzner et al., 2009). However, we demonstrate that higher frequencies are more spatially localized and have narrower stimulus tuning. Our estimate of spatial spread should be considered an upper bound as it is at the limit of the interelectrode spacing (which is still half that of Utah arrays).
We observed large anisotropy in the spatial spread of evoked CSEPs, with greater values dorsoventral than anteroposterior. This anisotropy thus reflects the functional organization (tonotopy) of the cortical tissue (Polley et al., 2007; Schreiner and Winer, 2007; Guo et al., 2012). Therefore, the major determinant of the spatial scale of correlations in stimulus-evoked CSEPS across the surface is not the passive diffusion of signals (which is expected to be isotropic (Łęski et al., 2013), but instead is correlations in function. Together, our results imply that maximization of ECoG grids should be matched to the spatial resolution of functional differentiation of the underlying cortex and subsampled relative to this to ensure robustness to electrode loss over time.
ECoG records electrical potentials reflecting weighted linear superpositions of all electrical sources in the brain (Lindén et al., 2011; Buzsáki et al., 2012; Łęski et al., 2013). Several studies have examined the relationship between intracortical spiking activity and the frequency content of recorded intracortical LFPs but report divergent results. Steinschneider et al. (2008), characterized the spectral content of A1-evoked LFPs and the relationship to MuA in nonhuman primates. Several studies have suggested that high-frequency content of intracortical LFPs (broadband activity) is directly modulated by neuronal spiking activity (Manning et al., 2009; Scheffer-Teixeira et al., 2013). Ray and Maunsell (2011) showed that evoked intracortical gamma band activity can be dissociated from high-gamma band activity, whereas a report from Leszczynski et al. (2020) suggests that broadband activity can be dissociated from MuA in different layers. However, these studies have primarily focused on intracortical recordings. Thus, little is known about the relative contributions of the superposed neuronal sources at the cortical surface.
Our study addresses this gap with a full-scale, biophysically detailed model (Markram et al., 2015) of cortical sources and a forward model of their superposition at the surface. The model produces stimulus-evoked mean responses with frequency content strikingly similar to that observed experimentally. The simulation predicted that as the magnitude of the input increases, there should be concomitant increases in both the frequency and amplitude of the evoked response with a concave shape. This prediction was validated with a novel observation in the experimental data. Previous efforts to understand CSEPs and intracortical LFPs produced by large populations of neurons use simplified neuronal models (Miller et al., 2009; Mazzoni et al., 2015) or simplified network models that omit several cell types or even entire cortical layers (Hagen et al., 2018). The Blue Brain Project model (Markram et al., 2015) used here, represents the state of the art in biophysically detailed simulation of neurons and cortical columns. However, the simulation did not perfectly reproduce the sharp onset in the raw CSEP or the variance at low frequencies. Nonetheless, together, our results indicate that for the first time a computational model accurately captures the biophysical processes giving rise to the evoked ECoG response observed in the data.
The biophysical model demonstrates that the intuition that CSEP signals must be generated primarily in superficial layers is incorrect. Instead, our results implicate neurons in cortical layer V as the primary source of the signal recorded at the surface, with layer VI also contributing substantially. Similarly, analysis of contributions to the CSEP by depth showed that slices of the column containing layer V somas produce most of the signal observable at the surface, not dendritic compartments. Thus, despite the fact that our model has detailed dendritic morphologies and voltage-gated Ca2+ channels, we did not find evidence of layer V dendritic compartments as major contributors to evoked high-frequency CSEPs, in contrast to previous results (Suzuki and Larkum, 2017; Leszczynski et al., 2020). Note that in line with the slow time course of dendritic Ca2+, the signal reported in Suzuki and Larkum (2017) is much slower than the high-gamma signal we focus on. The number and synchrony of neurons were found to be more important than depth in determining the contribution of a population to the surface signal. Although layers II–IV are closer to the recording electrode than layers V and VI, they have fewer neurons (Shepherd, 2004; Markram et al., 2015) and reduced synchrony (Atencio and Schreiner, 2013; Adesnik and Naka, 2018). Layers V and VI are composed of predominately pyramidal cells (Markram et al., 2015; Adesnik and Naka, 2018) whose spikes contribute most to high frequencies (Buzsáki et al., 2012). Thus, we conclude that evoked CSEP high gamma is a biomarker of layers V and VI pyramidal neuron firing rates. Our results also indicate that differences in laminar architecture across cortical areas and species may impact the precise origins of ECoG signals.
The simulation results suggest a word of caution for the interpretation of multiunit activity (and LFPs more broadly) recorded both at the surface and intracortically (Khodagholy et al., 2015; Trautmann et al., 2019; Leszczynski et al., 2020; Paulk et al., 2021). Previous reports of single-unit recording from the cortical surface as well as intracortically recorded multiunit activity may contain contributions from distal, but numerous and synchronous, neurons. Further, experimental analysis of contributions of spiking activity to CSEPs must use well-isolated single units, not multiunit activity, as has been done previously (Dougherty et al., 2019; Leszczynski et al., 2020). To provide mechanistic precision, we conducted experiments in anesthetized animals and likewise simulated a single cortical column. Although the stimulus-evoked responses of cortical neurons are well preserved from the anesthetized to awake states (Niell and Stryker, 2010), this simplification likely reduces the impact of, for example, top-down inputs that are present in the awake animal.
The finding that evoked ECoG high-gamma is primarily generated by synchronous neurons in layer V provides a potential explanation of the robust tuning to exogenous variables found here and previously (e.g., auditory stimuli, Bockhorst et al., 2018; vocal tract articulators, Bouchard et al., 2013, etc.). In particular, neurons in layer V have previously been found to have sharper tuning curves than neurons in layers II and III (Atencio and Schreiner, 2010b, 2012; Harrison et al., 2012). Thus, although cortical surface electrical stimulation may activate broadly connected neurons in layers II and III (Graziano et al., 2002; Harrison et al., 2012), recordings from the surface can reflect the finely tuned responses and precise projections of neurons in layer V (Atencio and Schreiner, 2010b; Harrison et al., 2012). The homogeneity of columns within an area (Mountcastle et al., 1997) and the linear properties of cortical tissue suggest that results here derived from µECoG will extrapolate to larger electrodes used clinically.
In summary, our results indicate that evoked ECoG high-gamma responses are primarily generated by the population spike rate of pyramidal neurons in layers V and VI of single cortical columns. Together, these results highlight the possibility of understanding how microscopic sources (specific neuronal populations) produce mesoscale signals (i.e., ECoG). For example, in some cases we observed a pronounced secondary peak at ∼375 Hz in the experimental data (Fig. 5b); this novel spectral component has not been previously characterized and was not reproduced by the simulation. We conjecture that this novel spectral component may reflect the activity of neurons in layers II and III in response to input from adjacent cortical columns (which were not explicitly modeled here). More broadly, we propose that different high-frequency components of ECoG signals reflect spiking activity of neurons in different cortical layers. As neurons in different layers perform distinct computations (Adesnik and Naka, 2018), this proposition implies that different components of ECoG signals could be biomarkers for these computations.
Footnotes
This work was supported by Lawerence Berkeley National Laboratory LDRD for the Neural Systems and Data Science Lab, National Institute of Neurological Disorders and Stroke Grant RNS118648A (K.E.B.), Marco Microelectronics Advanced Research Corporation Grant 2009-BT-2052, and National Science Foundation Grants EFRI 1240380 and DBI-1152658 (M.M.M.).
The authors declare no competing financial interests.
- Correspondence should be addressed to Kristofer E. Bouchard at kebouchard{at}lbl.gov
