Control of working memory by phase–amplitude coupling of human hippocampal neurons

Patients

A total of 36 patients (44 sessions; 21 female individuals; 15 male individuals; age: 40.47 ± 13.76 years; Supplementary Table 5; no statistical methods were used to predetermine sample sizes) participated in the study. All of the patients had Behnke–Fried hybrid electrodes (AdTech) implanted for intracranial seizure monitoring and evaluation for surgical treatment of drug-resistant epilepsy. Their participation was voluntary, and all of the patients gave their informed consent. This study was part of an NIH Brain consortium between three institutions (Cedars-Sinai Medical Center, Toronto Western Hospital and Johns Hopkins Hospital) and was approved by the Institutional Review Board of the institution at which the patient was enrolled. A pre-operative magnetic resonance imaging (MRI) image together with MRI or computed tomography post-operative images were used to localize the electrodes using Freesurfer as previously described³⁴. Electrode positions are plotted on the CITI168 Atlas Brain⁵⁷ in MNI152 coordinates for the sole purpose of visualization (Fig. 1b). The 3D plot was generated using FieldTrip (v.20200409) and the Brainnetome Atlas⁵⁸. Coordinates appearing in white matter or outside of the target area is due to usage of a template brain. Electrodes that were localized outside of the target area in native space were excluded from analysis (8 out of a total of 265 recording sites). Data used in this study are available in the DANDI archive⁵⁹.

Task

The task consists of 140 trials and 280 novel pictures. Each trial started with a fixation cross presented for 0.9 to 1.2 s (Fig. 1a). Depending on the load condition, the fixation cross was followed by either one (load 1; 70 trials) or three (load 3; 70 trials) consecutively presented pictures, each remaining on the screen for 2 s. In load 3 trials, pictures were separated by a blank screen randomly shown for 17 to 200 ms. Picture presentation was followed by a 2.55 to 2.85 s long maintenance period in which only the word “HOLD” was presented on the screen. The maintenance period was terminated by the presentation of a probe picture, which was either one of the pictures shown earlier in the trial (match) or a picture already presented in one of the previous trials (non-match; see below). The task was to indicate whether the probe picture matched on of the pictures shown earlier in the same trial or not. The probe picture was shown until the patients provided their response through a button press. The response mapping switched after half the trials, which was communicated to patients during a short half-time break. Responses were provided using a Cedrus response pad (RB-844; Cedrus). All of the pictures were novel (that is, the patient had never seen this particular image) and were drawn from five different visual categories: faces, animals, cars (or tools depending on the version), fruits and landscapes. Images (width × height: 10.5 × 7 visual degrees) were presented in the centre of the screen and never more than twice (that is, when serving as the probe picture). Pictures were only repeated when presented as the probe stimulus. To make sure that also the non-match probe pictures were never completely new to patients (as were the matching probe pictures), which could have been used as a strategy to solve the task without using WM, we always used a picture that the patients had seen already in one of the earlier trials, randomly drawn from one of the categories not used during encoding. If a patient participated in more than one session, we used a completely new set of pictures in each session to ensure that all pictures were novel in all of the sessions. The overall longer duration of load 3 as compared to load 1 trials ensured increased cognitive control demands in trials with higher load. The maintenance period was the same length regardless of load, and all analysis of neural activity during the maintenance period was performed within this time window (0–2.5 s after the maintenance period onset). Note that, in load 3 trials, the three encoded items were from three different categories, assuring that the participants always had to maintain pictures from three different categories. Thus, when comparing trials between load 1 and 3 for preferred trials, each load condition always contained exactly one item from the preferred category.

Spike sorting

Each hybrid depth electrode contained eight microwires from which we recorded the broadband LFP signal between 0.1 and 8,000 Hz at a sampling rate of 32 kHz (ATLAS system, Neuralynx; Cedars-Sinai Medical Center and Toronto Western Hospital) or 30 kHz (Blackrock Neurotech; Johns Hopkins Hospital) depending on the institution. All recordings were locally referenced within each recording site by using either one of the eight available micro channels or a dedicated reference channel with lower impedance provided in the bundle, especially when all channels contained recordings of neuronal spiking. To detect and sort spikes from putative single neurons in each wire, we used the semiautomated template-matching algorithm OSort (v.4.1)⁶⁰. Spikes were detected after band-pass filtering the raw signal in the 300–3,000 Hz band (single-cell quality metrics are shown in Extended Data Fig. 1). All analysis in this paper (including the LFP) is based on signals recorded from micro wires. We isolated 360 neurons in the hippocampus, 496 in the amygdala, 204 in the pre-SMA, 188 in the dACC and 206 in the vmPFC. Of the LFP channels, 586 channels were in the hippocampus, 421 in the amygdala, 283 in the pre-SMA, 325 in the dACC and 307 in the vmPFC.

LFP preprocessing

Before analysing the LFPs, we removed spike waveforms (action potentials) and excluded trials with interictal discharges and high-amplitude noise. First, to avoid leakage of spiking activity into lower frequency ranges^61,62, we removed the waveforms of detected spikes from the raw signal by linear interpolating the raw data from −1 to 2 ms around each spike onset in the raw recording before downsampling. As the same spike can, in rare instances, be recorded on more than one wire, we not only interpolated the data for the wire on which the neuron was detected but also for all other wires in the same wire bundle. We then low-pass filtered the raw signal using a zero phase-lag filter at 175 Hz and downsampled to 400 Hz. Line noise was then removed between 59.5 and 60.5 Hz as well as between 119.5 and 120.5 Hz using zero phase-lag band-stop filters. Extended Data Fig. 1f,g shows the raw LFP as well as the log–log power spectrum for an example channel from the hippocampus. The slope of log–log power spectra did not differ significantly between load 1 and load 3 trials in hippocampal channels (n = 586 channels; mean slope −1.7526 ± 0.3902 versus −1.7517 ± 0.3928, t₅₈₅ = −0.86, P = 0.39).

Artefact and inter-ictal discharge detection was performed on a per trial and wire basis using a semiautomated algorithm together with subsequent visual inspection of the data. To detect high-amplitude noise events as well as inter-ictal discharges, we z-scored the amplitude in each channel across all trials. To avoid artefactual amplitude biasing, we first capped the data at 6 s.d. from the mean and re-performed the z-scoring on the capped data^63,64. If a single time sample in each trial and wire exceeded a threshold of 4 s.d., the trial was removed from the analysis for that wire. Jumps in the signal were detected by z-scoring the difference between every fourth sample of the capped signal. Trials in which any jump exceeded a z-score of 10 s.d. were removed. The result of this cleaning process was visually inspected in every recording and any remaining artefactual trials were removed manually. If a wire or brain region showed excessive noise or epileptic activity, it was entirely removed from the analysis. On average, 20.4 ± 13.9 trials (14.6% of the data) were removed per wire.

PAC

We measured the strength of PAC for a wide set of frequency combinations in all of the recorded micro channels (except those excluded, see above) using the modulation index (MI) as introduced previously⁶⁵. As the cleaning process described above produced a different set of available trials for each channel, we first randomly subsampled from all correct trials in each channel such that the number of trials were the same for both load 1 and load 3. We then extracted the LFP starting at −500 until 3,000 ms following the maintenance period onset in each selected trial. We then filtered (using pop_eegfiltnew.m from EEGLAB, v.2019.1)⁶⁶ each trial separately within the respective frequency bands of interest (see below for more details). We then extracted the instantaneous phase from the lower-frequency signal and the amplitude from the higher-frequency signal using the Hilbert transform. Lastly, we cut each trial to the final time window of interest of 0–2,500 ms relative to maintenance period onset. This last step ensures that filter artefacts that arise at the edges of the signal are removed. All analysis of neural activity during the maintenance period was performed in this 2.5-s-long time window that started at the onset of the maintenance period. The length of the analysis window was the same in both load conditions. Next, we concatenated the phase and the amplitude signal across trials and computed the MI as described previously⁶⁵ (18 phase bins). We computed MIs separately for load 1 and load 3 trials. All subsampled trials from both load conditions were used to select for significant PAC channels in an unbiased fashion (see below). Example code to reproduce parts of the results in this study in published at Zenodo⁶⁷.

To standardize the MI in each channel, frequency and condition, we computed 200 surrogate MIs by randomly combining the phase and amplitude signals from different trials (trial-shuffling), again separately for load 1, load 3 and for all trials. We fit a normal distribution to these surrogate data (normfit.m) to obtain the mean and s.d. of each distribution. These values were then used to z-transform the raw MI values. Standardizing MI values eliminates potential systematic differences that might arise due to load-related power or phase differences, which could drive observed differences in PAC. Moreover, low frequencies are more vulnerable to non-specific correlations to high-frequency power due to non-stationarities in the LFP signal caused by factors such phase resets. Comparing raw modulation indices to trial-shuffled surrogates within the same condition will reduce PAC that is caused by such non-specific interactions (discussed in detail previously⁶⁸). In addition to providing a measure of significance, normalizing the MI values therefore allows for comparisons across conditions, frequencies and channels⁶⁹. A channel was indicated as having significant PAC present if the normalized MI computed across all subsampled trials (both loads) exceeded a z-score of 1.64 (P < 0.05, right-sided).

We repeated the above procedure for all frequency combinations. The phase signals were extracted for centre frequencies between 2 and 14 Hz in steps of 2 Hz (2 Hz fixed bandwidth). The amplitude signals were extracted for frequencies between 30 and 150 Hz in steps of 5 Hz. The bandwidth of the amplitude signals was variable and depended on the centre frequency of the low-frequency signal. It was chosen such that it constituted twice the centre frequency of the phase signal (for example, if combined with an 8 Hz centre frequency for the phase signal, the bandwidth of the amplitude signal was chosen to be 16 Hz). This procedure ensures that the side peaks that arise if the amplitude signal is modulated by a lower-frequency phase signal are included⁶⁸.

To determine the influence of theta waveform shape on PAC, we tested for differences in theta waveform peak-to-trough as well as rise-to-decay asymmetries between the two load conditions, which could potentially cause differences in TG-PAC^70,71. To extract and characterize each theta cycle during the delay period in all significant hippocampal PAC channels, we used the bycycle toolbox⁷² in Python. We averaged estimates for peak-to-trough as well as rise-to-decay asymmetries across cycles during the maintenance period from the same trials used for our PAC analysis within each load and tested the estimates between the conditions. Results of this analysis are presented in Extended Data Fig. 3c.

Moreover, we determined the number of significant PAC channels that showed theta–high-gamma nesting using the method described previously⁶⁹. To do so, in each PAC channel we determined the theta phase bin for which gamma amplitude was maximal, that is, the preferred theta phase of gamma amplitude. In the band-pass-filtered and phase-binned theta (3–7 Hz) signal, we then determined all instances during the delay periods of all of the correct trials in which this phase bin occurred, and extracted the precise timepoint at which the concurrent instantaneous gamma amplitude (70–140 Hz) was maximal within each bin. To obtain the average waveform, we selected a window of 500 ms centred on each timepoint in the raw (unfiltered) LFP recording and averaged the signals across all windows in each channel. Example average waveforms from two channels are shown in Extended Data Fig. 3e. In accordance with ref. ⁶⁹, we characterized a waveform as being nested if at least three local maxima fell within a window of 45 ms (that is, 3 cycles at 70 Hz) around the preferred phase. Results are presented in Extended Data Fig. 3e.

Relationship between single-trial PAC and FR or RT

We calculated single-trial estimates of TG-PAC for all significant PAC channels of both MTL regions and the vmPFC. We used mixed-effect GLMs to assess whether RT is related to PAC in a trial-by-trial manner (using only correct trials). We included load as a confounder and modelled random intercepts for each significant PAC channel nested into patientID. To examine whether there was a correlation between FR of category neurons (see below) and single-trial estimates of PAC, we used a mixed-effects GLM with load as a confounder and modelled random intercepts for each neuron to significant PAC channel combination. Only correct trials were used.

Category cell selection

We selected for neurons of which the response after stimulus onset during encoding differed systematically between the picture categories of the stimuli shown. To do so, for each trial, we counted the number of spikes a neuron fired in a window between 200 to 1,000 ms after stimulus onset (all encoding periods and the probe period). We then grouped spike counts based on the category of the picture shown in that trial. For each neuron, we computed a 1 × 5 permutation-based analysis of variance (ANOVA) with visual category as the grouping variable, followed by a post hoc one-sided permutation-based t-test between the category with maximum spike count and all other categories. We classified a given neuron as a category neuron if both tests were significant (P < 0.05, 2,000 permutations (see below)). We refer to the category with the maximum FR as the preferred category of a cell. To test whether the observed number of category cells was significantly larger than that expected by chance in each area, we repeated the above selection for 500 times after shuffling the category labels for each stimulus across all picture presentations. If the observed number of category cells in the unshuffled data was higher than the 99th percentile (P < 0.01) of the resulting shuffled distribution (which, across all five brain areas, corresponds to a Bonferroni-corrected alpha level of 0.05), we considered the number of category cells observed in a given area as significant. Note that category cells are selected only using spiking activity from the encoding period, leaving the FRs during the maintenance period independent for later analyses.

SFC

To measure how strongly the spiking activity of a neuron followed the phase of an LFP in a certain frequency, we computed the SFC. We measured SFC as the mean vector length (MVL) across spike-phases for all neuron-to-channel combinations available within a bundle or across regions (within the same hemisphere) in correct trials⁷³. To estimate the instantaneous phase from LFPs in different frequency ranges in each trial, we applied continuous wavelet transforms using 40 complex Morlet wavelets⁷⁴ spanning from 2 to 150 Hz in logarithmic steps. The number of cycles for each wavelet changed as a function of frequency from 3 to 10 cycles, also in 40 logarithmic steps⁷⁵. This ensured a higher temporal precision for longer wavelets at low frequencies and higher frequency precision for faster wavelets at high frequencies as compared to using a constant number of cycles across all frequencies⁷⁶. Extended Data Figure 1 shows the temporal and spectral characteristics across all wavelets used in this analysis. To assess the quality of our wavelet transform, we tested how well we were able to reconstruct the original signal after applying the wavelets to our data. To reconstruct the signal, we extracted the real-valued (bandpass-filtered) signal after applying each wavelet to the data and then summed up these signals across all frequencies. This resulted in a signal that closely followed the original recording in each trial (an example trial is shown in Extended Data Fig. 1j). We assessed how well the reconstructed signal predicted the original signal by computing R² values extracted from linear models using the reconstructed signal as predictors and the original signal as response variables in each trial and channel. As the quality of the reconstruction could change as a function of frequency or time, we performed this analysis for several time and frequency bins. First, we band-pass filtered both signals within the spectral bandwidth of each wavelet and then applied the linear model in sliding windows of 500 ms with a step size of 25 ms. The results of this analysis are presented in Extended Data Fig. 1.

For our SFC analysis, we first extracted data between −500 and 3,000 ms around the maintenance period onset from all clean trials in each channel and then computed a complex Morlet wavelet convolution to extract the instantaneous phase of the LFP as described above. The trials were then cut to the final time window of interest of 0 to 2,500 ms after the maintenance period onset to remove filter artefacts at the edges of each trial. To further avoid a bias of the MVL based on differences in spike count, we subsampled spikes such that an equal number of spikes was available in each condition. We included neurons that had at least 50 spikes available in each condition (we used a minimum of ten spikes for the preferred versus non-preferred analysis in category neurons due to a potentially low spike count in the non-preferred condition⁷⁷). Next, we extracted the phase in the LFP closest to the timestamp of each spike, averaged across all spike-phases in polar space, and computed the MVL for each of the 40 frequencies. We repeated this subsampling 500 times and averaged the resulting MVLs across all repetitions within conditions. To avoid potential bias of load within the preferred versus non-preferred (category neurons; Fig. 3) or fast versus slow RT SFC comparison (cross-regional analysis; Fig. 5), we computed the SFC estimates within each load condition and then averaged across the loads.

The resulting MVL in each neuron-to-channel combination was further normalized using a surrogate distribution, which was computed after adding random noise to the timestamps of all spikes within a condition 500 times. Potential biases of the MVL based on systematic differences between the conditions (such as power differences between conditions within a given frequency band) were thereby reduced. Like for the measure of PAC (see above), we fit a normal distribution to the surrogate data and used the mean and the s.d. of that distribution to z-score the raw MVL within each condition.

To compare SFC between preferred and non-preferred trials, we computed SFC for all category neuron-to-channel combinations within the same region in frequencies between 2 and 150 Hz during the maintenance period and compared trials in which preferred or non-preferred stimuli were correctly maintained. We used cluster-based permutation statistics to identify ranges of frequencies with significant differences (Fig. 3f). To determine whether the observed gamma SFC difference between preferred and non-preferred trials was dependent on gamma amplitude, we tested whether gamma SFC (averaged across 70–140 Hz) for category neurons in the hippocampus differed between preferred and non-preferred trials for high and low gamma amplitudes separately (median split).

Whether theta or gamma-band SFC was related to the preference of the cell and/or load was tested by averaging SFC within the theta (3–7 Hz) or gamma band (70–140 Hz) and computing a 2 × 2 permutation-based ANOVA with the factors load and preference for all category neuron to PAC channel combinations in the hippocampus.

To examine whether cross-regional SFC differed between the two load conditions, we computed SFC for all neuron-to-channel combinations between the respective areas in each load condition. We then used cluster-based permutation statistics to identify ranges of frequencies with significant differences (Fig. 5b; alpha level Bonferroni-corrected for all tests across two MTL areas, three frontal areas and two cell populations). To further determine a relationship to RT (Fig. 5f), we performed a median split of RTs for all correct trials within each load condition and compared cross-regional SFC between hippocampal PAC neurons and the vmPFC, averaged within and the significant theta range, between fast and slow RTs (averaged across both load conditions).

Selection of PAC neurons

We selected for neurons whose FR was correlated with both theta phase and gamma amplitude during the maintenance period of the task. For all neuron-to-channel combinations within a bundle of microwires, we extracted the data from correct trials between −500 and 3,000 ms relative to the maintenance period onset and estimated the phase of theta signals by filtering between 3 and 7 Hz and computing a Hilbert transform in each trial. Gamma amplitude was determined by computing wavelet convolutions for frequencies between 70 and 140 Hz in frequency steps of 5 Hz (each wavelet using 7 cycles). Trials were cut to 0 to 2,500 ms after maintenance period onset to remove edge artefacts, and were then concatenated. The extracted amplitudes in each gamma frequency were z-scored across all trials and averaged across all frequencies. Computing wavelet convolutions in 5 Hz steps and z-scoring the data before averaging avoided biasing power estimates to lower frequencies due to the power law. Next, for each neuron–channel pair, we performed a median split of gamma amplitudes and binned all amplitudes into low and high gamma, respectively. In each of the two gamma groups, we further binned the corresponding theta phases into 10 groups (36° bins), resulting in a total of 20 bins (Fig. 4a). In each of those bins, we then counted the number of spikes that occurred in each theta–gamma bin.

We fit three Poisson GLMs for each neuron-to-channel combination. In model 1, spike count (SC) was a function of theta phase (10 levels), gamma amplitude (2 levels), and the interaction between theta phase and gamma amplitude. We included theta separately as cosine and sine due to the circularity of phase values⁷⁸, which enabled us to treat theta phase as a linear variable. Model 2 included the theta phase and gamma amplitude as main effects but not the interaction term. Model 3 included a main effect for theta phase and an interaction term but no main effect for gamma amplitude:

We next compared pairs of models using a likelihood-ratio test between model 1 and the two other models (using compare.m). A neuron qualified as a PAC neuron if model 1 explained variance in spike counts significantly better than both of the other two models (P < 0.01, FDR corrected for all possible channel combinations). The rationale behind each model comparison was as follows. First, we were specifically interested in neurons that followed the interaction of theta phase and gamma amplitude, that is, PAC, and not just theta phase or gamma amplitude alone. Selecting neurons for which model 1, including the interaction term, explained spike count variance of a given neuron significantly better than model 2, lacking the interaction term, ensured extracting those neurons. Second, we also compared model 1 to model 3, lacking the gamma term, for the following reason. Assume that a given neuron–channel combination has an LFP with strong PAC at the field potential level, that is, strong interactions between theta phase and gamma amplitude, and a neuron of which the FR is not related to neither theta phase nor gamma amplitude. Nevertheless, this situation would result in a significant interaction term in model 1 because the spikes that fall into the low and high gamma amplitude groups will have different theta phases (due to PAC). This is only the case if the underlying PAC in the LFP is very strong (an illustration and further discussion is provided in Extended Data Fig. 6). However, in this scenario, the gamma amplitude term (or the theta phase term) would not be significant. Comparing model 1 to model 2 and model 3 therefore ensures that cells were selected only at PAC neurons in which the interaction term explained variance above and beyond the main effects and interactions alone.

As we did not observe strong PAC nor persistently active category neurons in frontal regions, we restricted this analysis to channels from the MTL regions and performed it separately in each load condition. If spike count variance was significantly better explained by model 1 than the two other models in either of the load conditions for at least one neuron-to-channel combination, we included this neuron as a PAC neuron. If a neuron was selected in more than one neuron-to-channel combination, we selected the combination with the highest R² in the full model (model 1). This combined channel was later used for within-region SFC as well as FR correlation analyses. Lastly, to determine whether the number of selected PAC neurons per area was significantly higher than chance, we repeated the entire selection process 200 times after pairing spikes and LFPs from different, randomly selected trials, therefore destroying their relationship with theta phase and gamma amplitude. The P values indicate the proportion of repetitions that resulted in a higher number of selected neurons using the shuffled data than the original number of PAC neurons determined using the unshuffled data.

Properties of PAC neurons

We used mixed-effects GLMs with load as a confounder and modelling a random intercept for each PAC neuron-to-channel combination nested into patientID (using only correct trials and the LFP channel selected for each neuron; see above) to examine the relationship between FR of PAC neurons and single-trial estimates of PAC. Note that trial-by-trial correlations are independent from the selection procedure as PAC neurons were selected on the basis of trial-averaged theta–gamma interactions, irrespective of their trial-by-trial variance.

Noise correlations and population category decoding

We investigated the effect of noise correlations among groups of simultaneously recorded neurons on population decoding accuracies for the image category currently held in mind and on WM behaviour during the maintenance period. To estimate noise correlations among pairs of category and PAC neurons, for each neuron, we counted spikes in bins of 200 ms that slid across the maintenance period (0–2.5 s after the last picture offset) in steps of 25 ms. We then computed the correlation coefficient across all 101 time bins in each single trial for each pair of neurons and averaged across all considered trials within each condition. We used only correct trials for this analysis, and paired only neurons that were recorded in the same session and within the same brain region. Pairs of neurons recorded on the same channel were not considered as a precaution against spurious correlations caused by spike sorting inaccuracies. To assess the significance of noise correlations among pairs of neurons, we shuffled trial labels within conditions, that is, within the preferred and non-preferred category as well as within each load, 1,000 times in each pair and recomputed the average correlation coefficient across all pairs. The original average correlation coefficient was then compared against the distribution of all average correlation coefficients obtained from the 1,000 trial shuffles (Fig. 6a (right)).

To investigate the contribution of PAC neurons to the population category decoding accuracy when noise correlations among neurons were intact or removed, we used an approach introduced previously³⁶ (Fig. 6b). To measure how much a single neuron affects the decoding accuracy of an ensemble of neurons, this approach finds optimized neuron ensembles that have maximal decoding accuracy by adding each single neuron to the ensemble in a stepwise manner. Each neuron’s contribution to the ensemble can thereby be determined. In more detail, using a linear decoder, first the decoding performance for each single neuron in each region is determined from all correct trials. The neuron with the best decoding performance is then paired with each remaining neuron to determine which pair yields the best decoding accuracy. This most informative pair of neurons is then again combined with each remaining neuron to determine the most informative triplet of neurons, and so on. These steps were repeated until all neurons were part of the decoding ensemble.

As we were most interested in decoding picture category from FRs in the maintenance period, we used trials from load 1 only. This is because the maintenance period in load 3 trials contains intermixed information about the three different categories maintained in WM. We trained a linear support vector machine (SVM) decoder (fitcecoc.m; one-versus-one) on 80% of trials and tested it on the remaining 20% using z-scored FRs. To ensure an equal amount of data for all five categories, we subsampled trials to match the lowest number of trials available in each stimulus category. Noise correlations among neurons were left intact by using the same trials for each neuron or removed by shuffling trials per neuron within each category. Shuffling trials within each category ensured that the original category label remained correct but correlations among neurons were removed. Any decoding benefit that is purely based on category-selective firing activity is therefore not affected (Fig. 6b (red)). Note that, if PAC neurons enhanced the decodability by ‘residual coding’ of category information, they should have done so also when noise correlations were removed through shuffling. We repeated each decoding analysis 500 times and averaged the results to generalize across trial selections.

To test the influence of PAC neurons as well as their noise correlations on decoding performances, we first tested contributions between intact and removed noise correlations for PAC neurons that were added to the ensemble before maximal decoding performance was reached in each session and area³⁶. This approach therefore tested the effect that single PAC neurons had on information encoding within a neural population (Fig. 6c). To determine a functional specificity of PAC neurons as a group, we further compared the maximal decoding performance before and after all PAC neurons were removed from the neuronal ensemble in each session. We did this for all sessions that had at least one PAC neuron, and at least two neurons left after removing all PAC neurons. We then compared this effect with removing the same number of non-PAC neurons from the ensembles (averaged across 500 iterations of random selections).

We quantified the effect of noise correlations on the geometry of the population response. The effect that noise correlations have on the encoded information in a population of neurons depends on the angle between the signal and the noise axis³⁷. To illustrate how the angle between the noise and the signal axes changes with noise correlations, we first simulated neural responses for a population of neurons that were partially tuned to two different categories (see Fig. 6f for a simulation of two neurons for which one was tuned and the other was untuned to category). Firing rates for each neuron were drawn from a normal distribution with variable variance. We simulated 200 trials for each category. For tuned neurons, a variable offset was added to the mean of one of the categories. To add noise correlations to the population of neurons, in each trial, we added a random number drawn from a normal distribution to the FRs of all neurons. To compare our simulation to a condition with removed noise correlations, we shuffled trials within each category to destroy noise correlations within conditions but leave signal correlations among neurons intact. We then determined the signal axis by training a linear SVM classifier on the FRs from all neurons and extracting the hyperplane (decision boundary) obtained from the model. The signal axis is defined as a vector orthogonal to that plane. The noise axis was determined by extracting the first principal component of the data across both categories.

We then quantified and compare the angle between the signal and the noise axis in the recorded data (Fig. 6g). For each recording session, we extracted the signal and the noise axis for the neuronal ensemble at which the difference in category decoding was maximal between removed and intact noise correlations. For this analysis, we included all sessions that had at least two hippocampal neurons available. To obtain the direction of the signal axis, we extracted the hyperplane from each of the ten trained binary SVM classifiers (trained on 80% of the data; one-versus-one decoding, see above) and derived a vector orthogonal to that plane using a QR decomposition. The noise axis was determined by extracting the first principal component of the data across categories. The resulting angle between the two vectors was determined and then averaged across all 10 binary learners and all 500 decoding repetitions, resulting in one angle per session separately for intact and removed noise correlations.

To further determine the functional specificity of PAC neurons, we projected the population responses onto the signal axes and determined the variance of the projection values before and after PAC neurons were removed from the ensembles at which the difference between intact and removed noise correlations was maximal (Fig. 6h). This analysis was performed for all sessions that had at least one hippocampal PAC neuron, and at least two neurons left after removing all PAC neurons from those ensembles. The rationale of the analysis was based on the idea that the variance of the projected values should be small when the angle between the noise and signal axis is large and vice versa³⁸. For each binary classifier, we projected the population responses for each category onto the signal axis and determined their s.d. We then averaged the obtained variances across both categories, all 10 binary classifiers and all 500 iterations, and tested the variances between intact and removed noise correlations before and after all PAC neurons were removed from the ensembles.

To compare noise correlations between fast and slow RT trials, we examined all possible PAC–category cell pairs in a given session (Fig. 6i). We analysed the trials in which the preferred or non-preferred categories of the category cell were held in WM separately. Fast and slow RT trials were defined by median split, computed separately in each load condition, and then averaged to avoid a bias of load in RTs. To assess the specificity of the fast versus slow RT trial difference to PAC neurons, we randomly paired category neurons with any other non-PAC neuron and compared noise correlations between fast and slow RT trials (for n = 162 randomly selected pairs; same n as for PAC-to-category neuron pairs).

The significance of population decoding (Extended Data Fig. 7a) was assessed by comparing the original decoding accuracy to a distribution of 1,000 decoding accuracies after randomly shuffling category labels.

Statistics

Throughout this Article, we use (cluster-based) nonparametric permutation tests (statcond.m as implemented in EEGLab, using option ‘perm’, or ft_freqstatistics.m in FieldTrip), that is, tests that do not make assumptions about the underlying distributions, or mixed-effects GLMs (fitglme.m in MATLAB) to assess statistical differences between conditions. In these tests, random permutations of condition labels were performed to estimate an underlying null distribution, which was then used to assess the statistical significance of the effect. The paired permutation t-tests that we performed are equivalent to computing pair-wise condition differences and testing the differences against zero. All permutations statistics used 10,000 permutations, and t-tests were tested two-sided unless stated otherwise. The corresponding t and F estimates, which are computed based on a normal distribution, are provided as a reference only. Bayes factors were computed using the BayesFactor package⁷⁹. BF₀₁ indicates the evidence of H0 (null hypothesis; no evidence between conditions) over H1. A value of 1 indicates equal evidence for H0 and H1, and values larger than 1 indicate more evidence for H0 over H1 and vice versa. SFC estimates tested across several frequencies were corrected for multiple comparisons using cluster-based permutation statistics as implemented in FieldTrip⁸⁰ with 10,000 permutations and an alpha level of 0.025 for each one-sided cluster, which was also Bonferroni corrected for the number of tests involved. Depending on whether we used z-scored FRs or spike counts, we used mixed-effects GLMs based on a normal or Poisson distribution, respectively. Finally, error bars shown in figures show the s.e.m. unless otherwise stated.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.