Author response: Spike-phase coupling patterns reveal laminar identity in primate cortex

Abstract

Full text Figures and data Side by side Abstract Editor's evaluation Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract The cortical column is one of the fundamental computational circuits in the brain. In order to understand the role neurons in different layers of this circuit play in cortical function it is necessary to identify the boundaries that separate the laminar compartments. While histological approaches can reveal ground truth they are not a practical means of identifying cortical layers in vivo. The gold standard for identifying laminar compartments in electrophysiological recordings is current-source density (CSD) analysis. However, laminar CSD analysis requires averaging across reliably evoked responses that target the input layer in cortex, which may be difficult to generate in less well-studied cortical regions. Further, the analysis can be susceptible to noise on individual channels resulting in errors in assigning laminar boundaries. Here, we have analyzed linear array recordings in multiple cortical areas in both the common marmoset and the rhesus macaque. We describe a pattern of laminar spike–field phase relationships that reliably identifies the transition between input and deep layers in cortical recordings from multiple cortical areas in two different non-human primate species. This measure corresponds well to estimates of the location of the input layer using CSDs, but does not require averaging or specific evoked activity. Laminar identity can be estimated rapidly with as little as a minute of ongoing data and is invariant to many experimental parameters. This method may serve to validate CSD measurements that might otherwise be unreliable or to estimate laminar boundaries when other methods are not practical. Editor's evaluation Authors demonstrate powerful methods that can be applied across species to find reliable markers that characterize activity in different cortical layers. Authors provide compelling evidence for these methods that enable systematic comparisons between slow extracellular voltage fluctuations and spiking across cortical columns. The results are timely since linear multielectrode array recordings have become a state-of-the-art technique. https://doi.org/10.7554/eLife.84512.sa0 Decision letter Reviews on Sciety eLife's review process Introduction Linear array electrodes have become a ubiquitous electrophysiological tool for understanding the functional roles of neural populations across the layers of the cortex, their interactions, and the computations they perform. This understanding requires reliable assignment of neurons to their respective laminar compartments. Precise localization of individual neurons can be obtained by electrolytic lesion to mark the position of an electrode; however, this procedure is not practical in experiments where the same animal is used in multiple experimental sessions as is nearly always the case in experiments involving non-human primates. The gold standard for identifying laminar compartments from functional recordings is current-source density (CSD) analysis (Mitzdorf, 1985; Mitzdorf and Singer, 1978; Schroeder et al., 1998). The CSD represents the second spatial derivative of local field potential (LFP) activity averaged to repeated events. These events are often sensory-evoked stimulation as in flashed visual stimuli in visual cortex (Maier et al., 2010; Schroeder et al., 1991; Wang et al., 2020) or short duration sounds in auditory cortex (Happel et al., 2010; Lakatos et al., 2007; Szymanski et al., 2011) which produce feedforward activation of the canonical cortical columnar circuit via activation of the input layer (Gratiy et al., 2011; Mitzdorf, 1985; Mitzdorf and Singer, 1978). Laminar compartments can then be assigned from the timing of subsequent patterns of current sources and sinks across electrode contacts on the linear array (Mitzdorf, 1985). The earliest current sink reflects the feedforward activation of the input layers (Mitzdorf, 1987; Mitzdorf, 1985), and the current sources above and below are used to estimate the boundaries with respect to superficial (layers 1–3) and deep (layers 5 and 6) cortical layers (Self et al., 2013). CSD analysis has been validated histologically (Mitzdorf and Singer, 1979; Schroeder et al., 1991; Takeuchi et al., 2011) and reliably captures known laminar differences in the functional properties of cortical layers, such as differences in evoked latencies across layers in response to driving input (Bode-Greuel et al., 1987; Einevoll et al., 2013; Plomp et al., 2017). There are some limitations on the use of CSD for identifying cortical layers in electrophysiology recordings. CSD analysis has largely been limited to primary sensory areas where the types of sensory stimuli that can robustly and reliably generate the evoked responses necessary for averaging CSDs are well established. Although CSD analysis has been used in higher-order cortical areas where sensory-evoked responses are apparent, such as visual responses in frontal cortex (Bastos et al., 2018; Godlove et al., 2014), it is less clear what to use as a triggering event in non-sensory cortical areas to reveal similar laminar patterns. Another limitation of CSD analysis is that depth estimates can only be taken as the average across the period of sensory response collection, making it difficult to track electrode drift during a recording. Further, noise in recordings due to bad channels, variability in filtering parameters, or ambiguity in the identifying source–sink pattern could potentially lead to the elimination of otherwise useful data or inaccurate laminar assignment. Without an alternative means of estimating the depth of electrode contacts, analyses of cortical circuit function are at risk of using biased definitions of laminar identity resulting in spurious conclusions about layer-dependent neuronal features. Here, we present a new method for determining laminar boundaries in cortical recordings based on spike-phase coupling patterns across linear array electrodes. Previous work has shown ongoing spiking activity is strongly coupled to the phase of LFP fluctuations in the cortex (Davis et al., 2022; Destexhe et al., 1999; Dotson et al., 2014; Eeckman and Freeman, 1990; Esghaei et al., 2018; Haegens et al., 2011) and these spike–field relationships can be diagnostic of circuit interactions in cortical systems (Safavi et al., 2023). We hypothesized that, as phase shifts occur across the layers of the cortex, the preferred phase angle of this spike–LFP relationship should be influenced by or reflect changes in laminar circuitry that contribute to differences in sources and sinks. We find that there is a phase reversal in spike–field coupling that reliably corresponds to the laminar boundary between input and deep layers as estimated by CSD. We therefore propose a novel methodology, called laminar-phase coupling (LPC), for identifying the laminar boundary between input and deep cortical layers. This pattern reliably reversed in awake recordings from both common marmosets and rhesus macaque in multiple cortical areas including extrastriate visual and prefrontal cortex. This method can be applied to the same data used for CSD analysis, but can also be applied to any cortical linear array recording data without the need for specific sensory stimuli as with CSD analysis. Data recorded under a variety of conditions, such as during behavioral tasks, spontaneous activity during fixation, or continuous data recorded in the dark all reliably reveal the same pattern of laminar-phase separation. LPC is largely invariant to filtering or referencing, and can be done on single- or multi-unit activity. The analysis does not require averaging and can be calculated online to estimate recording depth throughout the duration of an experiment. The preferred phase angle in the spike–field relationship is a simple yet effective measure for estimating contact depth in linear array recordings with respect to laminar boundaries and may help inform the study of cortical circuits in situations where CSD measurements are ambiguous or where CSD methods are impractical or ineffective. Results We recorded electrophysiology data from 32-channel linear probes (Atlas Neuroengineering) inserted perpendicular to the cortical surface into cortical regions V4 in two macaque monkeys and middle temporal (MT) or pre-frontal cortex (PFC) in two marmoset monkeys performing head-fixed visual tasks under various experimental conditions (Figure 1a). CSD analyses were performed using stimulus-locked LFP epochs (Figure 1b). These epochs were locked to stimulus flashes in the case of V4 and MT recordings or full field flashes in the case of PFC recordings. As standard in CSD analysis (Franken and Reynolds, 2021; Nandy et al., 2017; Schroeder et al., 1998; Self et al., 2013; Westerberg et al., 2022), we took from the resulting evoked source–sink patterns the earliest sink (red) as the location of the input layers and assigned our reference depth relative to the input layer as the bottom of the current sink, with the boundaries to the sources (blue) above and below identifying the boundaries of the superficial and deep cortical layers (Figure 1c, d). Figure 1 with 2 supplements see all Download asset Open asset Current-source density (CSD) analysis reveals laminar boundaries in cortical recordings. (a) Schematic of recording locations in two marmosets (PFC or MT) and two macaques (V4). Recordings were made with a 32-channel linear silicon probe in an acute recording chamber penetrating through a silicone artificial dura (AD). (b) Average evoked local field potential (LFP) responses to a stimulus flash (N = 141 trials) in an example macaque recording session in area V4. Each line is the response measured on a single contact. The depth is in absolute distance from the most proximal contact. (c) CSD measurement from the example recording session in b. The input layer is defined as the bottom and top of the earliest current sink (red), with the superficial and deep layers defined as above and below the input layer, respectively. Depth is measured relative to the bottom of the input layer. (d) Same as c, but in an example marmoset MT recording session (N = 225 trials). In order to test the relationship between spiking and LFP phase across the layers of the cortex, we first identified multi-unit spike times and measured the generalized phase (GP) of the LFP filtered from 5 to 50 Hz (Davis et al., 2022; Davis et al., 2020a). The use of a wider filter than traditionally used (i.e. 4–8 or 8–15 Hz) helps reduce phase distortions that occur when applying narrowband filters to signals that have broad spectral content (Figure 1—figure supplement 1). GP provides a measure of wideband phase that corrects for errors that may arise from applying the Hilbert Transform to signals with broader spectral content (Figure 1—figure supplement 2), although our results did not depend on our specific use of these techniques. We binned electrodes into superficial, input, and deep layers depending on their location relative to the source–sink pattern defined in the CSD analysis. We first calculated the degree of spike coupling on each channel to the LFP phase measured on the same channel grouped by cortical layer (Figure 2a). This was done by taking the length of the circular resultant of the spike-phase distribution, which we call the spike-phase index (SPI). This index value ranged from 0 (perfectly uniform spike-phase distribution) to 1 (all spikes occur at a single phase). The superficial channels had an average SPI value of 0.162 ± 0.012 (mean ± standard error of the mean [SEM]; N = 34 sessions; Figure 2b). This was significantly greater than the input layer (SPI = 0.122 ± 0.010 mean ± SEM; p = 0.0003, Wilcoxon signed-rank test; Figure 2c), and both were significantly greater than the deep layers (SPI = 0.073 ± 0.009 SEM; p = 0.000003 S. v. D. and p = 0.000009 I. v. D.; Figure 2d, e). The preferred LFP phase angle (i.e. the circular mean of the phase distribution) for spiking in the superficial and input layers was 2.86 and −2.68 rad, respectively, which corresponds toward the trough of ongoing LFP fluctuations. This was significantly different from the preferred phase angle for deep layer spiking (−1.44 rad; F = 38.52, p < 1 × 10−7 and F = 21.15, p = 0.00002, respectively; Watson–Williams test; Figure 2f) suggesting that indeed the deeper layers may be distinguished by a change in the preferred spike–LFP phase angle relative to the superficial electrodes and this could be read out from the spike–LFP relationship across the depth of the cortex. Figure 2 Download asset Open asset Within channel spike–local field potential (LFP) phase coupling strength and preferred phase angle varies across layers. (a) Spike–LFP phase coupling was measured by taking the phase of the LFP on a contact at the times when multi-unit spikes were detected on the same contact. (b) Spike-phase distribution for superficial contacts averaged across all recordings sessions (N = 34 sessions across 4 monkeys; error bars denote standard error of the mean [SEM]). The spike-phase distribution was strongly peaked toward ±π rad. (c, d) Same as b, but for contacts in the input and deep layers. (e) The average spike-phase index (SPI) was significantly weaker for input and deep layers relative to the superficial layer (0.16 vs. 0.12 and 0.07; p = 0.017 S. v. I. and p < 1 × 10−7 S. v D.; Wilcoxon signed-rank test; *p < 0.05, ***p < 0.0001). (f) There was no difference in the preferred LFP phase angle between superficial and input layers, but a significant difference between these layers and deep layers (−2.86 and −2.68 vs. −1.44 rad; F = 38.52, p < 1 × 10−7 and F = 21.15, p = 0.00002, respectively; Watson–Williams test). While the previous analysis examined the spike-phase relationship on each channel as a function of cortical depth on average, because of variability in the preferred mean phase on each channel (due to factors such as the strength of spike–LFP coupling and the number of spikes recorded), there was no significant correlation between the preferred phase angle and the depth of the electrode as measured with CSD (circular–linear correlation r = 0.28; p = 0.24). A more sensitive measure is to relate the preferred phase angle in the spike–LFP relationship across the depths of all channels in the cortex. This was achieved by calculating the SPI from the multi-unit spiking activity on a given channel relative to the LFP phase angle measured on each contact across the depth of the recording (Figure 3a). We call the resulting matrix the cross-channel LPC. As an example, the first matrix row is derived from the spikes recorded on the first channel compared against the phase measured from the LFP on each of the 32 channels. The second row is derived from the spikes on the second channel compared against the phase measured from the LFP on each channel, and so on. The result is a 32-by-32 matrix where each cell represents the relation between the spikes on one channel and the LFP phase on another channel. However, not all of these channels are in the cortex, so we realigned the matrix based on estimates of cortical depth from the bottom of the input layer as identified from CSD analysis. Figure 3 with 6 supplements see all Download asset Open asset Cross-channel spike–local field potential (LFP) phase coupling reversal correlates with putative input layer. (a) Schematic displaying cross-channel phase coupling analysis. The spike times on each channel are compared against the phase of the LFP on each other channel, yielding a matrix of spike-phase coupling. (b) Example recording session cross-channel spike-phase index (SPI) values from marmoset MT. Red dashed lines indicate the estimated depth of the bottom of the input layer from current-source density (CSD) analysis. (c) Preferred phase angles for the cross-channel spike-phase relationships in b. Spikes across all channels preferred ±π rad phase angles from LFP measure on superficial and input electrodes, but preferred 0 rad phase angles from LFP measured on deep phase angles. The LFP phase reversal aligns well to the estimate of the input layer from CSD analysis (red dashed line). (d) Grand average cross-channel SPI across all recording sessions in MT, V4, and PFC aligned to the putative input layer from CSD analysis (N = 34 sessions). (e) Grand average of the preferred phase angle for the data in d. (f) Scatter plot shows a significant correlation between the depth of the bottom of the input layer estimated from CSD analysis (x-axis) and the depth of the phase reversal (y-axis; Pearson’s r = 0.64, p = 0.00005). We first looked at the SPI values across the cross-channel LPC matrix. Figure 3b shows an example from a single recording (note that the diagonal in this matrix, from top left to bottom right, represents the approach in Figure 2). In this recording, we found that superficial spiking activity (defined based on CSD analysis) was strongly coupled to the LFP phase recorded on other superficial channels, which dropped off sharply below the input layer, and recovered when computed relative to the LFP phase on deeper channels. We next looked at the preferred phase angle for spiking activity on each channel relative to the LFP phase measured for different laminar compartments (Figure 3c). We found that spiking activity on all channels, regardless of cortical depth, tended to spike during ±π radian LFP phase angles measured on superficial channels. This phase relationship flipped such that spiking on all channels, regardless of cortical depth, tended to spike during 0 radian phase angles for LFP measured on deep channels. This phase reversal occurred about the channel we identified from CSD analysis as the boundary between the input and deep layers in the cortical column. In order to see if this pattern held across our recordings, we then aligned all of our recordings (N = 34) relative to the boundary that separated input and deep layers (from each recording’s CSD) and computed grand averages for SPI (Figure 3d) and preferred phase angle (Figure 3e). The average pattern across recordings showed the laminar-phase reversal was well aligned to the estimate of the input layer from the CSDs across recordings. This was also true when we broke out the recordings to only average across sessions in each cortical area (MT, V4, PFC; Figure 3—figure supplement 1). To quantify how well the cross-channel LPC and CSD aligned on each individual recording session, we compared the depth of the laminar boundary estimated from CSD analysis to the depth at which the preferred spiking phase angle reversed (Figure 3f). There was a significant correlation between the depth estimated from the CSD and the phase reversal (r = 0.64, p = 0.00005; Pearson’s correlation) and no significant difference between the estimated depth values (mean difference = 271 ± 49 μm SEM; p = 0.98, Wilcoxon signed-rank test). We next asked whether the observed phase relationship was parameter dependent. We first asked whether or not the specific use of a 5–50 Hz filter or the calculation of GP was necessary to see the laminar-phase reversal. To test this we filtered the data in a variety of commonly used frequency bands (4–8, 8–14, 15–30, and 30–50 Hz) and calculated phase using the Hilbert Transform (Figure 3—figure supplement 2). We found that the observed phase reversal was apparent in each case, indicating the results were not dependent on the filter or the method used to compute phase. We also tested whether the results depended on the use of multi-unit spiking activity. We performed the same analysis aligning cross-channel LPC from single units across recording sessions based on CSD depth (Figure 3—figure supplement 3). We observed the same relationship between laminar depth and preferred phase angle as when using multi-unit spiking activity. We also separated out different conditions during the same recording session where the electrode placement was unchanged (Figure 3—figure supplement 4). These included fixation during a visual task, freely viewing natural scene images, freely viewing in total darkness, and fixation during receptive field mapping. The cross-channel LPC revealed qualitatively similar depths for the laminar-phase reversal across experimental conditions within the same penetration. Finally, we explored the impact of referencing on the presence of the phase reversal. In our recordings, the electrode probe had a reference contact at the base of the shaft. In order to test whether our results depended on the location of the reference contact we re-referenced the LFP data to the shallowest contact, the deepest contact, or a common average reference (Figure 3—figure supplement 5). The phase reversal persisted under all referencing conditions. These results indicate the phase relationship is a robust feature of columnar recordings and not dependent on a particular set of parameters. One of the advantages of the LPC method is it can be done without the need for a triggering event, like a stimulus onset. Further, as little as a minute of continuously recorded data is sufficient to recover the phase reversal across channels. We find estimates can be derived from as few as hundreds of multi-unit spikes on each channel, although estimates are more reliable when spikes number in the thousands. As a result, depth estimates can be sampled at arbitrary points during recording sessions making this technique sensitive to tracking putative electrode drift over the course of a recording. To demonstrate this, we calculated LPC on sequential 3-min epochs during the first 15 min of an example recording session (Figure 3—figure supplement 6). We found the depth of the phase reversal moved over the course of 15 min. When comparing the depth after the first or last 2 min across all recordings, we found a significant difference in the estimated depth of the input layer from the spike-phase reversal (p = 0.005, Wilcoxon signed-rank test) such that the estimate of input layer depth consistently drifted deeper across recording sessions by an average of 132 ± 26 μm (SEM). This change in the depth estimate across the recording is consistent with electrode drift from movement of the probe relative to the cortex as the tissue settles at the start of the recording. The cross-channel LPC relies on both spikes and LFP, but if a similar phase pattern can be observed in the LFP alone, it may not require spiking activity to detect the laminar boundary. To test this we measured the circular–circular cross-correlation between LFP phase on each channel and LFP phase on each other channel. On many recording sessions we found a characteristic phase correlation pattern that seemed to align well with both CSD and LPC depth estimates (Figure 4a–c), and across many recordings the phase correlations within layers could dissociate laminar compartments. However, this was not always the case. Some sessions yielded phase correlation patterns that were noisy and difficult to interpret, although knowing where the boundary was made it easier to identify the characteristic pattern (Figure 4d, e). While, on average, LFP phase correlations can be informative (Figure 4f), the spike–LFP phase relationship appears to provide a more robust measure of contact depth. Figure 4 Download asset Open asset Local field potential (LFP)–LFP phase correlations can also reveal laminar boundaries, but are less consistent. (a) Preferred phase angle across all 32 channels from an example marmoset recording. The red dashed line indicates the current-source density (CSD) estimate of the bottom of the input layer. (b) Circular LFP phase correlation across channels. There is a strong within compartment correlation that aligns with the boundary between channels above and below the input/deep boundary (red dashed lines). (c, d) Same as panels a and b, but for a macaque recording. (e, f) Same as above, but an example of poor phase correlation within compartments despite strong phase reversal alignment with CSD. (g, h) Grand averages for preferred phase angle and LFP phase correlations across all recording sessions (N = 34). Plots were aligned to the putative input/deep layer boundary identified by CSD analysis. Across our recordings there were some instances (14/34) where there was 200 μm or more disagreement between CSD and LPC estimates for the location of the input layer relative to the probe contacts. Figure 5a shows an example CSD that has the characteristic source–sink pattern used to identify the input layer. This estimate was well aligned to the cross-channel LPC pattern previously described (Figure 5b, c). However, Figure 5d shows an example of a CSD profile with a similar source–sink pattern that is not well aligned to the cross-channel LPC pattern, which shows the characteristic phase reversal 600 μm below the CSD estimate (Figure 5d, e). In order to determine which measurement was more accurate in identifying cortical depth in these cases of disagreement, we sought to use each measure to replicate two findings from the literature that varied with cortical depth. Figure 5 Download asset Open asset Current-source density (CSD) analysis is not always consistent with laminar-phase reversal. (a) Example CSD from macaque V4 with estimate of input layer boundaries indicated by red dashed lines (N = 81 trials). (b) Cross-channel spike-phase index (SPI) for the same example recording as the CSD in a. Red line is the CSD estimate of the bottom of the input layer. Green line is the estimate from the phase reversal. (c) Cross-channel spike-phase pattern for same recording as in a and b. The phase reversal is well aligned to the CSD estimate of the input layer. (d) CSD from a different example recording session in marmoset MT (N = 661 trials). Laminar depth estimated the same as in a. (e, f) Cross-channel SPI and phase angle as in b and c. The CSD and phase reversal estimates do not align. First, previous work in multiple cortical areas has shown that the superficial layers have lower firing rates as compared to the input and deep layers (de Kock and Sakmann, 2009; Haegens et al., 2015; Lakatos et al., 2005; Leszczyński et al., 2020). We identified recording sessions where the CSD and LPC estimate of input depth differed by more than 200 μm (N = 14 sessions) and we measured the average firing rate as a function of depth relative to either the CSD or LPC estimate of the input layer (Figure 6a). We found that on these sessions, the CSD aligned firing rates failed to recapitulate the finding of higher firing rates in the input and deeper layers (p = 0.09, Wilcoxon signed-rank test). Conversely, we did find significantly stronger input firing rates when aligned to the LPC estimate (p = 0.0003, Wilcoxon signed-rank test), suggesting the LPC estimate on these sessions was better aligned than the CSD to the ground truth. Figure 6 Download asset Open asset Laminar-phase coupling (LPC) pattern better replicates spike rate and local field potential (LFP) power dissociations than current-source density (CSD) when the two measures disagree. (a) Normalized spontaneous firing rates across sessions when CSD and LPC disagreed as a function of contact depth by more than 200 μm (N = 14; shaded regions denote standard error of the mean [SEM]). LPC captures the expected higher spontaneous firing rates around the input layer. (b) Mean z-scored LFP power across sessions in lower frequencies (8–30 Hz, blue line) and higher frequencies (65–100 Hz, red line) as a function of electrode depth referenced to the CSD estimate for sessions. (c) Same as in a, but referenced to the depth estimate from the LPC phase reversal across sessions. The laminar power relationship is more pronounced when using the phase reversal instead of CSD estimate. (d) Scatter plot showing the disagreement in the depth of the input layer estimated from CSD (x-axis) and the depth from LPC (y-axis) in this subset of recording sessions. (e) Scatter plot showing the alignment of the CSD input depth (x-axis) crossover in LFP power (y-axis). (f) Same as e, but for LPC. The LPC depth measure was more correlated with the LFP power reversal than the CSD measure (Pearson’s r = 0.65 vs. 0.87 CSD vs. LPC). We next examined a second previously reported finding that varied with laminar depth. Previous groups have reported an inversion in the spectral content of the LFP across layers, with the superficial layers exhibiting more high-frequency power (e.g. 80–200 Hz) and the deep layers exhibiting more low-frequency power (8–20 Hz) (Maier et al., 2011; Maier et al., 2010). The crossover between the power in higher and lower frequencies was reported to occur around the input layers defined from CSD measurements (Bastos et al., 2018), and validated histologically in multiple cortical areas (Mendoza-Halliday et al., 2022). To test whether we could identify a similar relationship in these ambiguous recordings, and whether that relationship was stronger when aligned to CSD or LPC estimates of depth, we calculated low- and high-frequency power on each channel in each session and aligned them to each depth estimate (Figure 6b, c). The depth estimate from the LPC measure was significantly correlated with the high/low power inversion (Pearson’s r = 0.87, p = 0.00005) whereas the CSD estimate was more weakly, although still significantly correlated (Pearson’s r = 0.65, p = 0.012), indicating that on the recording sessions where there was some disparity between the CSD