Theta–Alpha Connectivity in the Hippocampal–Entorhinal Circuit Predicts Working Memory Load

Participants

We used data from 13 adult human patients [mean ± SD (range): 35 ± 13 (18–56); six females] in this study. All patients were implanted with depth electrodes (1.3 mm diameter, eight contacts of 1.6 mm length, and 5 mm spacing between contact centers; Ad-Tech, Racine, WI, www.adtechmedical.com) in the medial temporal lobe for evaluation of the surgical treatment of epilepsy. Electrode placement was exclusively guided by clinical needs. There were no seizures recorded during any of the recording sessions, and any trials with interictal epileptiform activity were excluded from analysis. All patients provided written informed consent before participating. This study has been approved by the relevant institutional ethics review board (Kantonale Ethikkommission Zürich, PB 2016-02055) and is in agreement with the Declaration of Helsinki.

Experimental design

We used a modified Sternberg task in which the encoding of memory contents, their maintenance, and their recall were temporally separated (Fig. 1A). Each trial started with a fixation period (1 s) followed by the stimulus for 2 s. Participants were instructed to memorize a set of four, six, or eight letters that were presented at the center of the screen. The number of letters indicated the memory load. After the disappearance of the stimulus, the maintenance interval followed (3 s). After the presentation of the probe letter, the participants responded with a button press (“IN” or “OUT”) to indicate whether the probe was part of the stimulus letter set. After the response, the participants were encouraged to relax before they initiated the next trial. The participants performed 50 trials per session, which lasted approximately 10 min. During the recording period of several days, several participants performed more than one session of the task up to a total of eight sessions. Table 1 contains detailed information about the number of trials and sessions for each participant. To eliminate any potential confusion arising from mixing categories with varying proportions of incorrect trials, subsequent analyses exclusively utilized trials with correct responses.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Working memory task, recording sites, and schematic of multivariate pattern analysis. A, An example trial of the task. Each trial consisted of a set of consonants, including four, six, or eight letters for 2 s (encoding). The number of letters was thus specific for the memory load (load 4, load 6, and load 8). After a delay for 3 s (maintenance, red), a probe letter was shown, and the participants indicated whether the probe was or was not contained in the stimulus letter set. B, Response accuracy decreased with load 4 (light red), 6 (light green), and 8 (light blue) across the group of participants. **p < 0.01. C, Channel locations of all participants in MNI152 space. Recording regions are indicated by different colors (red, hippocampus; blue, EC). D, Schematic of the multivariate pattern analysis. Phase locking value (PLV), Granger causality (GC) index, and phase–amplitude coupling (PAC) were calculated separately between the Hipp (red) and the EC (blue) for trials with WM load 4, load 6, and load 8. The patterns of PLV (e.g., theta/alpha PLV), GC (e.g., theta/alpha GC), and PAC (e.g., theta/alpha–gamma PAC) were separately used to train the SVM classifier to classify the WM load (load 4, load 6, or load 8). Specifically, we converted the features to a feature vector (M), fed them into a linear SVM classifier, trained the classifier on the data from N−1 participants (N−1 × M), and tested it on the remaining one participant (1 × M). We used inter-participant cross validation by leaving one participant out for validation and replicated the classification. The accuracy was used as the performance metric.

Channel localization and selection

Channel localization was performed using postimplantation CT scans and structural T1-weighted MRI scans. For each patient, the CT scan was co-registered to the postimplantation scan, as implemented in the FieldTrip toolbox (Oostenveld et al., 2011). The channels were visually marked on the co-registered CT-MR images. The channel positions were then normalized to the MNI 152 space and assigned to specific brain regions using the Brainnetome Atlas (Fan et al., 2016). Channel positions were verified by the neurosurgeon (L.S.) after merging pre-operative MRI with postimplantation CT images of each individual patient in the plane along the electrode (iPlan Stereotaxy 3.0, Brainlab).

For each participant we analyzed the iEEG from a maximum of three electrodes per hemisphere targeting the hippocampal head (anterior), the hippocampal body (posterior), and the entorhinal cortex. Targeted regions and hemispheres varied across participants for clinical reasons and included the hippocampus in the left (n = 12) and right (n = 13) hemisphere and the entorhinal cortex in the left (n = 12) and right (n = 11) hemisphere. We selected the two most medial channels on each electrode targeting the hippocampus or the entorhinal cortex as in previous studies (Oehrn et al., 2014; Pacheco Estefan et al., 2019). The final number of selected channels in each region for each participant is listed in Table 1. We included only ipsilateral channel pairs in the analysis. The final data set contained 87 channels in the hippocampus and 46 channels in the entorhinal cortex across all patients. There were 6.7 ± 1.5 (range 4–8) channels per patient in the hippocampus and 3.5 ± 0.9 (range 2–4) channels per participant in the entorhinal cortex. Figure 1C presents all recording locations with the BrainNet Viewer toolbox (Xia et al., 2013) in MATLAB (MathWorks).

Data acquisition and preprocessing

Intracranial data were acquired using a Neuralynx ATLAS recording system, sampled at 4 kHz and analog-filtered above 0.5 Hz. The data were recorded against a common reference where the outermost electrode contact in the temporal cortex served as an electrical reference. After data acquisition, neural recordings were downsampled to 1 kHz and band-pass filtered from 1 to 200 Hz using the zero-phase delay finite impulse response (FIR) filter with a Hamming window. Line noise harmonics were removed using a discrete Fourier transform. The filtered data were manually inspected to mark any contacts or epochs containing epileptiform activity or artifacts for exclusion and were then re-referenced. The continuous data were segmented into 8 s trials with a 1 s fixation period as the baseline, 2 s encoding period, 3 s maintenance period, and 2 s retrieval period. Here, we focused on the maintenance period. The trials with residual artifacts were rejected after a visual inspection. In total, we rejected 65 trials with load 4 (5.5%), 36 trials with load 6 (3.9%), and 39 trials with load 8 (5.1%) across all participants. Preprocessing routines were performed using the FieldTrip toolbox (Oostenveld et al., 2011) and customized scripts in MATLAB.

Statistical analyses

To assess the significance of a value, we created a null distribution estimated from 1,000 permutations on data with scrambled labels using a non-parametric permutation test. The significance was defined as exceeding the threshold that obtained from the 95th percentile of the empirically estimated null distribution.

To compare the metrics of phase synchronization, GC, PAC, and PSI between two loads (load 4 vs load 6, load 4 vs load 8, load 6 vs load 8), the statistical significance was then estimated using a permutation test, in which a null distribution was created by randomly assigning trials into two loads, computing the differences between loads, and repeating this procedure 1,000 times. We also applied paired t tests to directly compared measurements from the two directions within each load condition. Multiple comparisons were corrected by false discovery rate (FDR). p < 0.05 was considered significant.

To rule out potential confounding effects of aperiodic activity on WM load, we first performed separate repeated-measures analyses of variance (ANOVAs) for the hippocampus and the EC, with a within-group effect of load across the frequency spectrum. The aperiodic activity was chosen as it captures non-oscillatory effects at a specific frequency, allowing us to examine the precise frequency ranges sensitive to load effect on non-oscillatory components. Additionally, we applied a linear-mixed effect model to explore whether the aperiodic activity from the hippocampus and the EC contributed to the hippocampus–EC interaction, with the load, the aperiodic activity as fixed factors, the participants as a random factor, and the electrophysiological indexes as dependent variables.

For all the decoding analyses, we used a non-parametric permutation approach to test the significance of the accuracy values. We created a null distribution of the decoding accuracy by shuffling the data labels 1,000 times. For each decoding analysis, the null distribution was generated for each test, and we took the maximum value of the null distributions across tests as a final null distribution for multiple corrections, as a previous study did (Mamashli et al., 2021). The averaged decoding accuracy exceeding the 95th percentile of such null distribution (p < 0.05) was considered significant.

Time–frequency analysis

Time–frequency power was separately computed within the hippocampus and the EC for trials with different WM loads. For each trial and each channel, we convolved the signal with complex-valued Morlet wavelets (six cycles) to obtain power information at each frequency from 1 to 100 Hz (in steps of 1 Hz) with a time resolution of 1 ms. The task-induced power was analyzed per trial using a statistical bootstrapping procedure [methods have been described in more detail in our previous publication (Li et al., 2022)]. Then, the raw power for each time point during the task was z-scored by comparing it to the null distribution to generate the z-scored power. For each participant, the z-scored spectral power in the theta/alpha band was averaged across the maintenance period within the hippocampus and the EC separately for each WM load.

A previous study reported that periodic properties of electrophysiological data are highly variable and also coexist with variable and dynamic aperiodic activity (Donoghue et al., 2020a,b; Donoghue and Watrous, 2023). To exclude possible confounding effect of aperiodic activity on neural oscillations of the hippocampus and EC, we have adopted a distinct approach by separately characterizing the aperiodic properties of power spectra originating from both the hippocampus and the EC, with the Fitting-Oscillations-and-One-Over-F (FOOOF) toolbox (Donoghue et al., 2020b). We isolated the aperiodic component of power spectra across the entire frequency spectrum for each load and region and compared these components among WM load for each region with the repeated-measures ANOVA.

Phase locking value

To explore the potential interaction between the hippocampus and the EC during WM maintenance, we employed a phase locking value (PLV) to assess the degree of consistency for each channel pair phase relationship independent of their absolute phases and amplitudes among repeated trials, with PLV = 1 referring to strong phase synchrony where all trials are synchronized without any variations between the two channels. Using the same parameters of time–frequency analysis, we computed the PLV in the time–frequency domain from 1 to 100 Hz during maintenance for each channel pair within the same hemisphere (one channel from the hippocampus and one from the EC) for the trials with each WM load.

To evaluate the statistical significance of PLV, a null distribution was created by randomly shuffling the trials with load 8 for each channel pair and computing the corresponding PLV spectrogram and repeating the same procedure for 1,000 times, as our previous study did (Boran et al., 2019). Then, the null distribution of all channel pairs was averaged and only the time–frequency PLVs above the threshold (95% of the null distribution) were kept as significant PLVs for further analyses. In addition, the workload-dependent increases in the PLV were subsequently assessed by subtracting the PLV for trials with one load from the PLV for another load, and the statistical significance between two loads was then estimated using a permutation test, as mentioned in the section of statistical analyses.

Granger causality analysis

After establishing the phase synchrony, which measures undirected connectivity between the hippocampus and EC, we proceeded to investigate the directionality of their interaction using two complementary measures: a frequency-domain GC and PSI. GC measures the degree to which the signal from a region (i.e., the hippocampus) can be better predicted by incorporating information from another signal (i.e., the EC) in a specific frequency band, and vice versa (Zheng et al., 2019). For each channel pair, the trial-wise mean was subtracted from each trial before fitting to an autoregressive model and computing the spectral GC. We then used the Multivariate Granger Causality MATLAB Toolbox (Barnett and Seth, 2014) based on the Akaike information criterion to define the model order for each pair. Based on the observations from the PLV above, we computed the GC index across 2–12 Hz (in steps of 0.25 Hz) for both directions (from the hippocampus to the EC and the reverse direction) with the trials of load 4, load 6, and load 8, separately. To test the statistical significance of GC, we created a null distribution by randomly shuffling the signal between the channel pairs 1,000 times and averaged the null distribution of all channel pairs. Only the value above the threshold (95% of the null distribution) was kept as significant GC for further analyses. For the GC index across 2–12 Hz from the hippocampus to the EC as well as the opposite direction, we also applied the permutation test for comparisons between two loads, as mentioned in statistical analyses. To rule out the bias of the aperiodic activity, with the linear mixed-effect model, we considered WM load, the aperiodic activity of the hippocampus, and the EC within 2–12 Hz as fixed factors and the participants as a random factor. The hippocampal driven GC as well as the EC driven GC across 2–12 Hz was set as the dependent variable.

Phase slope index analysis

On the other hand, PSI examines whether the slope of the phase differences between the channel pairs remains consistent across several adjacent frequency bins (Nolte et al., 2008). A positive PSI signifies that the channel in the first structure (e.g., hippocampus) leads the channel in the second structure (e.g., EC), whereas a negative PSI indicates the reverse. For the trials with WM load, using the FieldTrip toolbox (Oostenveld et al., 2011), the data segments during maintenance were zero padded and multiplied with a Hann taper from 2 to 12 Hz with 1 Hz step, from which we computed the theta/alpha PSI at each channel pair within the same hemisphere in each participant (i.e., one from the hippocampus and the other from the EC) and pooled all possible channel pairs between the hippocampus and EC for each participant. To correct for any spurious results, we randomly shuffled the trials and recomputed the PSI at each channel pair. This step was repeated 1,000 times to create normal distributions of channel pair-resolved null PSI data.

To construct a directional effect of the hippocampus–EC on a population level, we averaged the raw PSI across channel pairs and participants. Correspondingly, the null distributions were also averaged across channel pairs and participants. Consequently, the raw PSI outputs can be compared to the distribution of null PSI to derive a z-score in the theta/alpha band (for a similar approach, Solomon et al., 2019). To examine if the null distribution of PSI by randomization is a normal distribution, we assessed the normality of the null distribution for different WM loads using the Jarque–Bera test, a widely used statistical test that examines the skewness and kurtosis of a sample to determine its normality. The null distribution of PSI for load 4, load 6, and load 8 is normally distributed (Jarque–Bera test: load 4, p = 0.13; load 6, p = 0.11; load 8, p = 0.50). As a result, raw PSI outputs were z-scored and significant PSI was thresholded at |z| > 1.96, in which the hippocampus leads were defined as z > 1.96 and the EC leads as z < −1.96 as in our previous study (Li et al., 2022, 2023).

To assess the statistical significance of PSI differentiation for WM load, the permutation test described previously was also used here to create a null distribution of PSI differences between two loads. Then the real PSI differences were obtained between loads and were then compared with the corresponding null distribution to estimate a z-score with the positive value indicating PSI increase in the high load condition versus the low load condition.

Cross-regional PAC

Next, we investigated the cross-regional PAC between the hippocampus and the EC using the modulation index (MI) (Tort et al., 2010; Vaz et al., 2017), which reflects the coordinated activity between brain regions. We first calculated the time series of phase and amplitude envelope. This was achieved by applying the standard Hilbert transform to the low-frequency (2–30 Hz) phase, extracted from the hippocampus/EC channel, using a 2 Hz step. Similarly, the high-frequency (30–150 Hz) amplitude was obtained from the channel in the other structure (EC/hippocampus) with a 5 Hz step. Subsequently, we partitioned the continuous phase values of the modulating frequency into 20 evenly spaced phase bins. For each phase of the low-frequency modulating signal, we determined the corresponding amplitude of the high-frequency modulated signal and assigned it to the respective phase bin. To assess this coupling, we employed the MI that quantifies the disparity in entropy between the computed phase–amplitude distribution and a uniform distribution using a normalized Kullback–Leibler distance between each pair of low modulating frequencies and high modulated frequencies.

To assess the statistical significance of PAC, we generated a null distribution with a trial shuffling procedure. Specifically, we created shuffled versions by associating the phase series of trial k with the amplitude series of trial l, with k and l randomly chosen among the trial numbers. We then generated 1,000 surrogate MI values, from which we could infer the MI chance distribution. To construct a directional effect of PAC on a population level, we averaged the raw PAC across channel pairs and participants. Correspondingly, the null distributions were also averaged across channel pairs and participants. Consequently, the raw PAC outputs can be compared to the distribution of null PAC to derive a z-score in the phase and amplitude frequency band. Significant zPAC was thresholded at |z| > 1.96.

We also examined the relationship between cross-regional theta/alpha–gamma zPAC and WM load with the permutation tests for comparisons between two loads. Thresholding was performed at the 95th percentile level, as stated in the statistical analyses section. Similarly, as done in GC analysis, the linear mixed-effect model was applied to examine the effect of the aperiodic activity on zPAC. Specifically, we treated WM load, the aperiodic activity of the hippocampus (2–12 Hz) and of the EC (30–100 Hz) as fixed factors, and the participants as a random factor. The hippocampal theta/alpha phase–EC gamma amplitude zPAC and the opposite direction were set as the dependent variables. We also tested the impact of phase synchrony to the functional effect of zPAC within the same model, as a previous study conducted (Wang et al., 2021).

Machine learning analyses

In addition to conventional univariate analysis, multivariate analysis detects subtle load-related distribution pattern changes missed by univariate methods (Grootswagers et al., 2017), and enhances findings’ generalizability and reliability through inter-subject validation (Wang et al., 2020). Therefore, we investigated next whether the neural activity and inter-regional communication within the hippocampal–EC circuit were modulated by WM load. We used the patterns from PLV, GC, zPAC, and z-scored power from the trials with WM load as our features. Here we used support vector machine (SVM) (Chang and Lin, 2011) as a classifier to classify the WM load (load 4/6/8). SVM is widely used in decoding analyses in neuroimaging studies (Mamashli et al., 2021) because of its suitability for analyses with a relatively small number of samples. It is provided by the COSMOMVPA package (Oosterhof et al., 2016) in MATLAB. And the approach of leave-one-out cross-validation (LOOCV) was applied to validate the decoding accuracy. Considering the inherent difficulty of generalizing across different subjects (Poldrack et al., 2009; Poldrack, 2011), LOOCV is shown to be a suitable method for obtaining dependable accuracy estimates, especially when working with data sets that have a restricted number of samples (Wong, 2015). Details of our multivariate pattern analysis (MVPA) decoding analyses were as follows:

1. PLV patterns: We considered the theta/alpha (2–12 Hz) PLV patterns between the hippocampus and the EC during maintenance to investigate whether the phase synchrony pattern could decode the WM load. For each participant and each load, there were M features (11 frequency bins × 3,000 time bins) converted to a feature vector. We trained a linear SVM classifier and applied LOOCV at subject level by splitting the data set of all subjects (N = 14) into a training set from N−1 subjects and a testing set from the remaining one subject. This process was repeated N times to ensure comprehensive validation. For each iteration, we used the feature vectors labeled as load 4, load 6, and load 8 from the training data set, which encompassed the data of N−1 participants. This resulted in a training data set including (N−1) participants × 3 sets × M features. Subsequently, we calculated the average classification accuracy by averaging the results across the N repetitions of the cross-validation procedure. Meanwhile, to reduce the feature dimensionality, principal component analysis (PCA) was applied to the training data set to keep several principal components (K components) that explained 99% of the variance in the data. We also applied the K components matrix on the remaining data set from one participant and tested the SVM classifier. This procedure was replicated N times for cross-validation. The schematic of the MVPA using the feature patterns is shown in Figure 1D. The accuracy of the classifier was averaged across N cross-validations as a measure of performance. To test the significance of the accuracy, we created a null distribution by shuffling the training labels 1,000 times. And the averaged decoding accuracy exceeding the 95th percentile of the null distribution (p < 0.05) was considered significant.

2. GC patterns: We considered the GC patterns from two directions, from the hippocampus to the EC and the reverse direction, to allow us to investigate whether there was a specific information flow pattern that could decode the WM load. The GC patterns were calculated in the theta/alpha band separately for trials with different WM load. For each participant and each load, the GC pattern included M values (M = 41), and these values were converted into a feature vector. As described above, we used the feature vectors labeled as load 4, load 6, and load 8 from N−1 participants as a training data set (N−1 × 3 sets × M features) and tested these on the remaining one participant data set. Similar to the LOOCV performed in the previous analysis, we left one participant out for validation and replicated this procedure N times. The accuracy of the classifier was averaged across all replications. In total, we separately performed this classification process for the two directions: for the hippocampus modulating the EC and for the EC modulating the hippocampus. Similar to the PLV patterns, we generated a null distribution with each direction of GC patterns and took the maximum value of the two null distributions as the final null distribution for multiple corrections, as a previous study did (Mamashli et al., 2021). The averaged decoding accuracy exceeding the 95th percentile of the null distribution (p < 0.05) was considered significant.

3. Z-scored PAC patterns: We calculated the theta/alpha phase-gamma band zPAC between the hippocampus phase modulating the EC amplitude and the opposite direction differences for the trials of load 4, load 6, and load 8 separately. For each participant and each load, there were M features (11 phase bins × 36 amplitude bins) converted to a feature vector. Similar as described for the PLV patterns, we combined N−1 participants’ data set from trials with WM load as training data set, applied zPAC to the training data set to K components that explained 99% of the variance, fed the features (N−1 × 3 sets × K components) into a linear SVM classifier and trained the classifier, and tested it on the remaining one participant data set that already applied K components matrix to the testing data set. The accuracy of the classifier was averaged across N replications. To test the significance of the accuracy, we created a null distribution by shuffling the training labels 1,000 times. And the averaged decoding accuracy exceeding the 95th percentile of the null distribution (p < 0.05) was considered significant.

4. Z-scored power patterns: To address whether local activity in the hippocampus and the EC contributed to WM load, we used a frequency specific z-scored power pattern at the theta/alpha band (2–12 Hz, 11 frequency bins) during maintenance from the hippocampus and EC to decode the WM load. The training data set for the linear SVM classifier from N−1 participants (N−1 × 3 sets × M features) and the classifier was tested on the remaining one participant data set. The accuracy of the classifier was averaged across N replications by LOOCV. We performed two classifications (two regions) in this decoding analysis, and the statistical analysis that was performed aligns with the above analyses.