Research ArticleNEUROSCIENCE

Persistent hippocampal neural firing and hippocampal-cortical coupling predict verbal working memory load

See allHide authors and affiliations

Science Advances  27 Mar 2019:
Vol. 5, no. 3, eaav3687
DOI: 10.1126/sciadv.aav3687


The maintenance of items in working memory relies on persistent neural activity in a widespread network of brain areas. To investigate the influence of load on working memory, we asked human subjects to maintain sets of letters in memory while we recorded single neurons and intracranial encephalography (EEG) in the medial temporal lobe and scalp EEG. Along the periods of a trial, hippocampal neural firing differentiated between success and error trials during stimulus encoding, predicted workload during memory maintenance, and predicted the subjects’ behavior during retrieval. During maintenance, neuronal firing was synchronized with intracranial hippocampal EEG. On the network level, synchronization between hippocampal and scalp EEG in the theta-alpha frequency range showed workload dependent oscillatory coupling between hippocampus and cortex. Thus, we found that persistent neural activity in the hippocampus participated in working memory processing that is specific to memory maintenance, load sensitive and synchronized to the cortex.


Working memory (WM) describes our capacity to prospectively store sensory input and translate it into an appropriate behavioral response. This capacity is necessary for several higher-order cognitive functions, such as reading and dialing a telephone number (1). If the input is represented in WM for prospective use, then how are these representations maintained in the brain?

The basic process that underlies the temporary storage of information in WM is persistent neuronal firing in a widespread neural network. Sustained neuronal firing—while information is maintained during the delay before a response—is a hallmark of WM (2). To date, single-neuron recordings in humans found persistently active neurons during WM maintenance for neurons with high stimulus selectivity (concept neurons) that showed, however, no workload dependence (3, 4).

At a systemic level, WM processing is known to correlate with sustained neuronal oscillations in the theta-alpha range (3 to 12 Hz) (511). The anatomical basis of WM involves a widespread network of brain areas, as shown noninvasively with electroencephalography (EEG) (5, 6, 1012) and functional magnetic resonance imaging (fMRI) (13), and invasively with intracranial EEG (iEEG) (79, 14, 15). The hippocampus has been proposed to be a subcortical node associated with frontal theta oscillations (16). Evidence for hippocampal involvement in WM originated from iEEG studies that revealed changes in theta and gamma power (15) and cross-frequency coupling that was load dependent during periods of WM activation (17). However, other studies using simultaneous EEG-fMRI (18) and fMRI (19) have failed to identify evidence supporting hippocampal involvement in the maintenance of WM; thus, the involvement of the hippocampus in WM is far from clear.

Synchronized oscillations have been proposed as a mechanism for functional interactions within the neuronal networks that underlie cognitive tasks (20, 21). These oscillations support phase-locked high-frequency firing within local networks and temporal coupling of the low-frequency phase for long-range communication between cortical areas (9, 10, 12, 22). This synchronization suggests an active maintenance process through reverberating signals between brain regions.

To provide direct evidence of hippocampal involvement in the WM neural network, we recorded single-neuron activity from microelectrodes and iEEG in the medial temporal lobe (MTL, which includes the hippocampus, entorhinal cortex, and amygdala) simultaneously with scalp EEG during the presurgical evaluations in patients with epilepsy. We used a WM task known to elicit frontal theta and parietal alpha oscillations during the maintenance period (5, 6) to clarify both the involvement of the hippocampus in WM and the mechanism of the functional hippocampal-cortical interaction. We hypothesized that hippocampal neurons will show load-dependent firing in WM. Since previous studies demonstrated that highly stimulus-selective cells in the hippocampus did not exhibit load dependence (3, 4), we tested load dependence without the constraint of high stimulus selectivity using letters in a modified Sternberg WM paradigm. We also expected that neuronal firing in the MTL would distinguish periods of the trial in WM and that oscillatory coupling between hippocampal iEEG and scalp EEG would be associated with WM load during maintenance.


Task and behavior

Subjects performed a modified Sternberg WM task (36 total sessions from nine subjects). The task differs from the original Sternberg paradigm in that the items were presented all at once rather than sequentially, thus separating the encoding period from the maintenance period. In each trial, subjects were instructed to memorize a set of four, six, or eight letters presented for 2 s (encoding). The number of letters was thus specific for the memory load. After a delay (maintenance) period of 3 s, a probe letter prompted the subjects to retrieve their memory (retrieval) and to indicate by button press (“IN” or “OUT”) whether or not the probe letter was a member of the letter set held in memory (Fig. 1A). Throughout this paper, the term “workload” refers to the number of letters held in memory. In their effort to solve the task, the subjects reported that they relied on audio-verbal transformation of the letters. Although stimuli were visually presented in this task, letter strings are thought to activate verbal WM through the phonological loop (1).

Fig. 1 Task, behavioral results, and recording sites.

(A) In the task, sets of consonants were presented and had to be memorized. The set size (4, 6, or 8 letters) determined WM workload. In each trial, presentation of a letter string (encoding period, 2 s) was followed by a delay (maintenance period, 3 s). After the delay, a probe letter was shown, and subjects indicated whether the probe was presented during the encoding period. (B and C) Behavioral results. (B) Accuracy of all sessions. The jitter reflects the rank order of the accuracy in trials with set size 8. Each dotted line connects an individual session. The nine subjects performed a total of 36 sessions. (C) Median reaction time (relative to onset of the probe letter) as a function of workload. The jitter reflects the rank order of the reaction time in trials with set size 8. Each dotted line connects an individual session. In (B) and (C), thick and light blue lines represent the means and SEM across all sessions, respectively. (D) Location of the microelectrodes at the tip of the depth electrodes in MNI’s MNI152 space (see Methods). Recording locations are projected on the parasagittal plane x = −25.2 mm and are color-coded (cyan, hippocampus; magenta, entorhinal cortex; yellow, amygdala).

The average correct response rates were 90.2 ± 7.8% for IN and 93.1 ± 4.3% for OUT. The rate of correct responses decreased with set size from a set size of 4 (98.5% correct responses) to set sizes of 6 (90.5%) and 8 (84.7%) [permuted repeated-measures analysis of variance (ANOVA), F2,70 = 40.25, P = 0.0002; Fig. 1B]. Across all sessions, the capacity averaged 6.8 [Cowan’s K, (correct IN rate + correct OUT rate − 1) × set size], which indicates that the subjects were able to maintain at least six letters in memory. The mean response time (RT) for the correct trials (1630 trials) increased with workload (48 ms per item; permuted repeated-measures ANOVA, F2,2007 = 21.10, P = 0.0002; Fig. 1C). Correct IN/OUT decisions were made more rapidly than incorrect decisions (1.36 ± 0.64 s versus 1.87 ± 1.21 s, respectively; permutation t test, t = 5.06, P = 5.0 × 10−5). Together, the subjects performed the task within normal limits (7 ± 2 items) and easily performed the low-workload trials (set size 4) (5, 23). In summary, these data show that our subjects were able to perform the task and that the difficulty of the task increased with set size.

Single-unit electrophysiology

To find out whether hippocampal neurons participated in this task, we performed multielectrode recordings of broadband neuronal activity from microelectrodes separated into single- or multiunit activity and local field potentials (LFPs). We refer here to a putative unit by the term “neuron.” We recorded the activity of 1526 (51 ± 35 per session) neurons from the hippocampus (n = 676), entorhinal cortex (n = 449), and amygdala (n = 401) across all microelectrodes. Neurons active during multiple sessions entered the analysis independently for each session. We quantitatively assessed the spike sorting quality and identified the mean firing rate of all neurons as 2.55 ± 3.30 Hz (fig. S1). Figure 1D shows the microelectrode recording sites projected on a parasagittal plane in the Montreal Neurological Institute (MNI) space. In our first analysis of individual neuron types, we focused on the maintenance period, when stimuli were absent, and the retrieval period after the presentation of the probe letter. To reduce noise in the classification, unless stated otherwise, subsequent analyses included only trials with correct responses.

First, we identified neurons that fired persistently in the maintenance period. These maintenance neurons were defined by their higher firing rates during the maintenance period than during the fixation period (permutation t test; Fig. 2A; further examples in fig. S2). We identified maintenance neurons in all recorded areas (Fig. 2B); however, they occurred most frequently within the hippocampus (χ21 = 47.30, P = 6.1 × 10−12), where 20.6% of all recorded neurons were maintenance neurons (permutation test against scrambled data, P = 0.002). Were the firing rates of maintenance neurons related to workload? During maintenance, a subset of maintenance neurons in the hippocampus (15%, 21 neurons; permutation test against scrambled data, P = 0.002) increased their firing rates systematically with workload (set size 4 versus set sizes 6 and 8; Fig. 2C). Across MTL, the hippocampus stood out with a significant percentage of maintenance neurons with persistent activity that correlated with the WM load.

Fig. 2 Persistent activity of maintenance neurons.

(A) Example of a maintenance neuron recorded from hippocampus. Top: Peristimulus time histogram (bin size, 500 ms; step size, 20 ms). Shaded areas represent ± SEM across trials of all spikes associated with the neuron. Inset: Mean extracellular waveform ± SEM. Middle: Periods of significance (black) between low-workload trials (set size 4, 45 trials) and high-workload trials (set sizes 6 and 8, 46 trials; P < 0.05, cluster-based nonparametric permutation test). Bottom: Raster plot of trials reordered to set size and RT for plotting purposes only. Compared to set size 4 (blue), the neuron fires more for set size 6 (green) and set size 8 (red). Figure S1 provides further examples of single neurons. (B) Percentage of all recorded neurons identified as maintenance neurons in the hippocampus (Hipp), entorhinal cortex (Ent), and amygdala (Amg). The number of maintenance neurons is provided below the area label. The hippocampus contained significantly higher proportions of maintenance neurons than the two other areas. (C) Percentage of maintenance neurons for which the firing rate during maintenance differed as a function of load. The hippocampus stood out as containing a significant percentage of neurons that increased their firing rate under high workloads. For (B) and (C), significance was assessed by comparing with a null distribution with permuted labels. ***P = 0.002.

Next, we identified neurons, which fired specifically after the presentation of the probe letter that initiated the retrieval period. These probe neurons were defined based on increased firing only during the presentation of the probe stimulus relative to encoding and maintenance (permutation t test; Fig. 3). We identified probe neurons in all areas and, most prominently, in the entorhinal cortex (χ21 = 23.16, P = 1.5 × 10−6; Fig. 3E), in which 9% of neurons were probe neurons (41 neurons; permutation test against scrambled data, P = 0.001).

Fig. 3 The activity of probe neurons is related to WM retrieval.

(A) Example probe neuron recorded from the entorhinal cortex. Firing rates are shown separately for trials in which the probe was held (IN, cyan) or not held (OUT, magenta) in memory. Top: Peristimulus time histogram (bin size, 200 ms; step size, 20 ms; shaded areas represent ± SEM). Bar: Time with significant differences between IN and OUT trials (P < 0.05, cluster-based nonparametric permutation test). Bottom: Raster plot of reordered trials, black marker at button press. Button press was at median RT 1.1390 s after probe letter presentation. (B) Same neuron as in (A), with trials aligned to button press. The peak response was reduced (27 Hz versus 18 Hz; permuted t test, P = 0.005). (C) Firing rates of an example probe neuron recorded from the entorhinal cortex, shown separately for low-workload trials (set size 4, blue) or high-workload trials (set size 6, green, and set size 8, red). Bar: Time with significant differences between trials of low and high workloads (P < 0.05, cluster-based nonparametric permutation test). Bottom: Raster plot of reordered trials, black marker at button press. (D) Same neuron as in (C), with trials aligned to button press (button press was at median RT 1.1390 s after probe letter presentation). The peak response was reduced (15 Hz versus 10 Hz; permuted t test, P = 0.005). (E) Percentage of probe neurons in each area. Probe neurons were most frequent in the entorhinal cortex. Significance was assessed by comparing with a null distribution with permuted labels. ***P = 0.001.

Was the response of these neurons related to either probe letter presentation or movement initiation? We compared the maximal firing rate after alignment to the probe onset with the rate after alignment to the button press. Aligning the response with the probe stimulus onset resulted in significantly higher peak firing rates (mean for all probe neurons, 5.47 ± 4.97 Hz versus 4.28 ± 4.37 Hz; permutation t test, t75 = 4.40, P = 0.0005; Fig. 3, A and B). Thus, firing of probe neurons reflects the processing of the probe stimulus more than direct preparation for movement. However, only a small number of probe neurons (eight neurons, 10.5% of probe neurons) reflected the subject’s decision (IN versus OUT; Fig. 3, A and B), and only one neuron showed workload dependence (Fig. 3, C and D), while the substantial majority of the probe neurons did not reflect specific trial demands. Thus, the responding neurons may represent the perception of the probe letter and the subsequent change in behavioral demand.

Attractors in neural population dynamics

We subsequently focused on the neural population firing rate during the periods of the trial. We used demixed principal component analysis (dPCA) (24) to project the firing rates from all neurons onto a low-dimensional component space. To demix the effect of stimulus identity, dPCA was informed by the workload of the trials. We formed a three-dimensional space from the first three demixed principal components (dPCs 1, 2, and 3, explaining 44.7% of the variance; Fig. 4A, fig. S3, and movie S1). dPCA distinguished between workloads 4, 6, and 8. This speaks for a population activity that distinguished between workloads 6 and 8 very early during the trial, since the projection on the dPC2-dPC3 plane shows three angles of 120° each, i.e., the optimal balanced distinction between the three workload conditions.

Fig. 4 Population firing predicts behavior.

(A) Mean trajectories in the neuronal state space constituted by the three largest dPCs during fixation (starting at the origin), encoding, maintenance, and retrieval. Set sizes are color-coded (4, blue; 6, green; 8, red) (see also fig. S3 and movie S1). (B) Multidimensional speed of the population in the four periods of a trial. The speed during maintenance was slowest. (C) Multidimensional pairwise distance between all possible pairs of attractors during the different periods of the task. The distance during maintenance was largest. The combination of low speed and large mutual distance provides evidence for attractors in state space during maintenance. The analysis in (B) and (C) included the first 15 dPCs that explain 79% of the signal variance. Significance was assessed by permutation t test. ***P = 0.0005. (D) During the encoding period of trials with high workloads (set sizes 6 and 8), the activity of neurons in the hippocampus (but not entorhinal cortex or amygdala) carried information about whether the subject would later respond correctly or not. (E) During maintenance, workload (set size 4 versus set sizes 6 and 8) could be decoded from hippocampal neurons. Sufficient decoding was possible even when using only those hippocampal neurons that we had identified as maintenance neurons. (F) During retrieval (−500 to 0 ms relative to button press), the subpopulations in all recorded areas were predictive of the response IN or OUT. In (B) to (F), boxplots represent quartiles (25%, 75%); horizontal lines represent medians; whiskers show ranges up to 1.5 times the interquartile range; dots outside whiskers show outliers. Markers below bars indicate significance versus chance performance. Significance was assessed on the basis of a null distribution with scrambled labels (gray boxes). ***P = 0.002.

In this neural state space, we first analyzed the rate of change over time (speed; Fig. 4B). The speed was highest during encoding, when maintenance neurons ramped up their firing. By contrast, speed was lowest during maintenance (permutation t test against the other trial periods, P = 0.0005). When comparing the two periods without visual stimulation, the speed during maintenance was much lower than during fixation. Analyzing the pairwise distance between trajectories that corresponded to different set sizes, we found that the distances during maintenance were larger than during encoding (permutation t test, P = 0.0005; Fig. 4C). Thus, during maintenance, the observed pattern of the firing rates of all recorded neurons resembled an attractor as a location in state space. Over the trial, the neural trajectories distinguished the periods within trials, among which the maintenance period stood out in clustering neural activity around attractors. Hippocampal population activity distinguished between trials of set sizes 4, 6, and 8 and was thus indicative of workload.

Decoding of information during the trial

So far, our findings did not answer the question of whether the neuronal activity observed was related to the subject’s behavior. Therefore, we analyzed neuronal population firing on a trial-by-trial basis. A population analysis of single trials seems adequate in this cognitively demanding task, which, during maintenance, involves internal processing in the absence of external stimuli. For this analysis, we created pseudopopulations of neurons that were recorded separately in different subjects; however, they were treated as if they were recorded simultaneously. Furthermore, we selected subpopulations that consisted of, e.g., only hippocampal maintenance neurons. We subsequently trained a decoder on a subset of trials and tested its performance on an independent set of test trials (25). Thus, we assessed whether neuronal firing was sufficiently salient to be representative of task demands or performance in single trials. Given the substantial variability in neuronal dynamics during a trial (Fig. 4A), we performed separate analyses for the different periods within the trials.

For the encoding period (−5 to −3 s), we found that the activity of hippocampal neurons, but not neurons in other areas, could be used to decode whether the subject would later reply correctly or incorrectly in the trial. The median decoding accuracy was 80% (permutation test against scrambled data, P = 0.002; Fig. 4D). This finding suggests that hippocampal neurons carry information on the memory content of the trial.

For the maintenance period, we found that the workload (set size 4 versus set sizes 6 and 8) of each trial could be decoded from hippocampal neurons but not from neurons in the other areas (permutation test against scrambled data, P = 0.002; Fig. 4E). The median decoding accuracy was 80%. The decoder did not distinguish between trials of set sizes 6 and 8, which points to a ceiling effect as illustrated by firing rates of the neuron in Fig. 2A. To ensure that some subjects did not drive the results, we performed a leave-one-out analysis. The exclusion of any subject changed the decoding accuracy to a value in the range [71 to 81%]. We also performed the same analysis using only maintenance neurons; there was little difference in the decoding accuracy between the use of only maintenance neurons and the use of the whole population of hippocampal neurons. Conversely, the exclusion of maintenance neurons reduced the prediction accuracy to 71%. This evidence indicates that hippocampal neurons carry information on the memory workload of the trial and that maintenance neurons carry the principal part of this information.

For the retrieval period, we found that neurons from all areas encoded whether the subject would later reply by pressing the IN or OUT button in the trial (permutation test against scrambled data, P = 0.002; Fig. 4F). The median decoding accuracy reached 90% if firing from all neurons was included, thus indicating that all regions contributed to the response selection during retrieval from WM.

Electrophysiology of cortico-hippocampal functional coupling

Since it was thus established that hippocampal neurons process memory workload and contribute to retrieval from WM, we next examined how the hippocampal neuronal activity was involved with the network of cortical areas activated during the task. To this end, we recorded oscillatory activity in the EEG from scalp electrodes and the iEEG from macroelectrodes implanted in the hippocampus, entorhinal cortex, and amygdala.

In our iEEG and EEG data, we here first show the results for selected electrodes and electrode pairs that illustrate patterns present in the data in subject 1 (Fig. 5). The depth electrode positions for iEEG are summarized in Fig. 1D, and the positions of the scalp electrodes were based on the 10-20 system. We then extend these analyses to the connectivity pattern across the entire dataset in Fig. 6.

Fig. 5 Neuronal responses in the WM network in subject 1.

(A) The task induced high alpha power in the scalp EEG at electrode site Pz, predominantly during maintenance (only trials with set size 8 in this analysis; n = 55). (B and C) During the last 2 s of the maintenance period, EEG power was higher for high-workload trials (set size 8, red line) than for low-workload trials (set size 4, blue line; n = 61 trials) in (B) the band (8.5 to 16 Hz, black bar) at electrode site Pz and (C) the band (8 to 12.5 Hz, black bar) at electrode site Fz. (D) The task increased the PLV between the hippocampal iEEG and the scalp EEG at electrode site Pz specifically in the alpha band and during the last 2 s of the maintenance period (−2 to 0 s, only trials with set size 8). (E) During this period, the PLV between the hippocampus and Pz was higher for trials with high workloads (set size 4 versus set size 8) in alpha (9.0 to 12.0 Hz, black bar), and (F) the alpha PLV was maximal to occipitoparietal scalp electrodes. (G) Theta-alpha power (7.5 to 10.0 Hz) in the left hippocampal iEEG increased with workload during maintenance. (H) During maintenance, the hippocampal iEEG power for set size 8 (red) exceeded that for other set sizes in the 7.5- to 10.0-Hz band (normalized to fixation, corrected for multiple comparisons). (I) The average firing rate of maintenance neurons in this subject was higher for a high workload during encoding and maintenance (n = 5 hippocampal neurons with more than 100 spikes during both encoding and maintenance). (J) Spike-field coherence between the LFP and the neurons from (I) appears in the alpha frequency range during maintenance (thin line, individual neurons; thick cyan line, average across neurons) but not during encoding (thick gray line, average across neurons). (B, C, E, and G to I) Black bar: Frequency range of increase with workload (set size 8 versus set size 4); magenta bar: frequency range of significant PLV; golden bar: frequency range of gating (maintenance versus encoding). P < 0.05, cluster-based nonparametric permutation test.

Fig. 6 Maintenance-induced hippocampal-cortical PLV.

(A to H) PLV spectra for the electrode pair of maximal PLV for subjects 2 to 9; for subject 1, refer to Fig. 5E [set sizes: 4 (blue), 6 (green), and 8 (red)]. High PLV between hippocampus iEEG and scalp EEG appeared for set size 8 in all subjects during the last 2 s of the maintenance period (magenta bar: frequency range of significant PLV with P < 0.05, randomization test against a null distribution with scrambled trials). PLV was higher for high-workload trials (set size 8) than for low-workload trials (set size 4) in eight of nine subjects (black bar: frequency range of elevated PLV with P < 0.05, randomization test against a null distribution with scrambled labels). The PLV was higher during the last 2 s of the maintenance period than during encoding (golden bar: frequency range of elevated PLV with P < 0.05, randomization test against a null distribution with scrambled labels). (I) Maximal PLV for each subject during retention plotted against the location of the electrode site along the anterior-posterior hippocampal axis. Electrodes in left (right) hippocampus are marked by circles (squares). Marker face colors denote subjects. The linear fit shows a gradual increase in the PLV toward the posterior end of the hippocampal axis (34 electrode pairs; R2 = 0.118, P = 0.0467, permutation t test).

As illustrated for subject 1 in Fig. 5, at the cortical level, the scalp EEG power increased during maintenance for a high workload (set size 8), with a peak within the 8- to 12-Hz alpha range over parietal cortex (Fig. 5, A and B). This workload-dependent EEG power increase extended to frontal electrode Fz (Fig. 5C). While this power increase is known from healthy subjects for the same or similar task (5, 6), it did not appear in the other subjects of our study group (fig. S5 and table S1). At the hippocampal level, the iEEG power increased in the alpha range with workload during the maintenance period (Fig. 5, G and H). At the neuronal level, the mean firing rate of the maintenance neurons in this subject distinguished between low and high workloads (set size 4 versus set size 8) during the maintenance period (Fig. 5I). The firing rate became strongest before the onset of the test stimulus, similar to the time course of power in scalp EEG and hippocampal iEEG (Fig. 5, A and G). Last, the neuronal action potentials showed an enhanced spike-field coherence to the LFP in the alpha range during maintenance (Fig. 5J).

Thus, we found workload-dependent increases of scalp EEG power in frontoparietal areas, hippocampal alpha power, and the firing rate of hippocampal maintenance neurons. But how are hippocampal and cortical activities related? To help in answering this crucial question, we examined the interregional synchronization with the phase-locking value (PLV) as a measure of functional connectivity between hippocampal iEEG and scalp EEG. The WM task enhanced the hippocampal-cortical PLV specifically in the 8–10 Hz alpha range and during the maintenance period (Fig. 5D and fig. S4). The PLV increased with workload (Fig. 5E) and was highest at occipitoparietal electrode sites (Fig. 5F). The hippocampus thus appears to be a node in the brain network activated during the WM task in this subject.

Cortico-hippocampal low-frequency coupling in the group of subjects

We then extended the PLV analysis to all subjects. We included the electrode contacts that targeted the hippocampus and all scalp electrodes. For the maintenance period, we calculated the PLV between the scalp EEG and the iEEG and then identified the electrode pair with the largest significant PLV for each subject. All subjects showed a significant PLV peak for at least one pair of signals recorded from a scalp electrode and a hippocampal channel (table S1). The frequency band of the maximal significant PLV indicated individual differences; we used the frequency bins of significant PLV to define an individual frequency band for each subject (permutation test, P < 0.05), which had its lower bound in the theta-alpha band across subjects and extended up to 17.5 Hz (median band limits, 6.0 to 11.5 Hz; Fig. 6 and table S1). In the majority of the subjects, there were no differences in the iEEG or EEG power in the individual frequency band during maintenance, which suggests that the differences in PLV were primarily a result of the differences in phase consistency.

Was the workload dependence of PLV (Fig. 5E) present across the group of subjects? During maintenance, the PLV was increased for high workload in eight of nine subjects in the individual frequency band (set size 4 versus set size 8; permutation test, P = 0.05; range, 5.5 to 16.5 Hz; median band limits, 7.8 to 10 Hz; Fig. 6 and table S1). In the evolution of the PLV over time (as illustrated for subject 1 in Fig. 5D), the PLV during maintenance exceeded the PLV during encoding in eight of nine subjects (Fig. 6 and table S1). Thus, the functional connectivity between the hippocampus and cortex in the theta-alpha frequency range increased specifically during maintenance and with the cognitive demand of the task.

To look for even more localized contributions of subregions within both hippocampi, we selected the scalp electrode from the pairs with the maximal PLV and plotted the maximal PLV in the individual frequency band for all hippocampal contacts (Fig. 6I). There was no evidence that the PLV differed between the left and right hippocampus. For the subject group, there was a gradual increase in the PLV from the anterior end to the posterior end of the hippocampal axis during maintenance (R2 = 0.118, P = 0.0467) and the PLV to posterior sites (y < −20 mm) was higher than to anterior sites (permutation t test, t = 2.09, P = 0.0405). This posterior preference might be associated with the functional organization proposed for the hippocampal longitudinal axis (26).

Recordings in the epileptogenic zone

In our analysis, we also included electrodes in brain regions that were ultimately resected. Among biomarkers that are thought to indicate the epileptogenic zone, we analyzed high-frequency oscillations (HFOs) in all patients (27). The rates of clinically relevant HFOs were not affected by the performance of the task. This points to a clear dissociation between the ripples observed during cognitive tasks in rodents and the clinically relevant HFOs that are used to predict seizure outcome in the individual patient (27).

When comparing hippocampal channels within and outside the seizure onset zone (SOZ) (table S1), the increase in the PLV was not significantly different between the non-SOZ electrodes and SOZ electrodes. Although epilepsy may have reduced functional connectivity to SOZ channels (28), we observed task-dependent connectivity, suggesting residual function in the epileptogenic hippocampus.

Spike-field coherence

Last, we tested whether the action potentials of maintenance neurons were temporally aligned with the LFP. For this analysis, we included all maintenance neurons that had at least 100 spikes during both maintenance and encoding. We found significant spike-field coherence to the LFP but not to the iEEG or the scalp EEG. For all frequency bins, the percentage of neurons with spike-field coherence was higher during maintenance than during encoding (fig. S6). This finding suggests that during maintenance, the functional coupling in the WM network was enhanced.


We recorded scalp EEG, temporo-medial iEEG, and single-unit activity from nine epilepsy patients as they participated in a WM task. Our main findings refer to neural firing in the MTL and coupling between scalp EEG and iEEG. The trial-by-trial variability in neural firing in the hippocampus (i) differentiated between success and error trials already during the encoding period and (ii) predicted workload during maintenance, while (iii) firing in the entorhinal cortex, hippocampus, and amygdala predicted the behavioral response during the retrieval period. In addition, workload enhanced the coupling between EEG and hippocampal iEEG, particularly (iv) during WM maintenance and (v) in theta-alpha frequencies, and (vi) coupling increased with the size of the letter set to be memorized. We interpret this activity as reflecting enhanced functional connectivity between the cortex and hippocampus during the maintenance of information about WM. We thereby showed the interaction between the hippocampus and previously identified frontoparietal sites engaged in WM, which was only feasible using the combination of modalities used here.

Neuronal firing patterns during the trial

Within the MTL, we identified different roles for neuronal firing in the hippocampus, entorhinal cortex, and amygdala during the different task periods.

Encoding. First, during encoding, hippocampal population firing predicted whether subjects later responded correctly or incorrectly (Fig. 4D). Several studies (3, 4) have identified individual neurons that represent a stimulus with a specific concept (e.g., a familiar item on a picture). In our study, we could not identify concept neurons. The stimulus material used in our task included letters, which are highly familiar but represent a far more abstract concept than a picture and convey minimal emotional content. This aspect might explain why the encoding of a set of letters is not assigned to neurons that form attractors for the individual items in memory (3). Rather, hippocampal neurons here seem to participate in generic processes of WM, possibly by establishing a contextual index in the network with distributed cortical sites involved in letter processing (2).

Maintenance. Our main question was whether neuronal activity during maintenance predicts workload. We found maintenance neurons that fired in the absence of stimulation (Fig. 2). The maintenance neurons in hippocampus even predicted the number of items held in memory in each trial (Fig. 4E). This is remarkable because neurons that code for the numerosity of memory content were known only in frontal cortex (3). We now found workload-dependent neurons also in the MTL.

Our finding of workload dependence seems at variance with a previous finding that hippocampal neurons were not predictive of the number of items in memory (3). The different types of stimuli may explain this difference between the two tasks. While the maintenance of concepts carried by visual representations in the hippocampus has been nicely demonstrated (3, 4), there is also evidence that the hippocampus is capable of ordering sequences of a number of simple items (29). Simple items like letters may be retained using rapid sequence coding (30) and, after audio-verbal transformation, rehearsal in the phonological loop (1, 7). This is less likely to be the case for complex pictures that evoke elaborate semantic representations and of which fewer items can be maintained. The results of these studies point to a distinction between simple load-dependent and complex concept-dependent maintenance in WM. Recent studies suggest that the longitudinal axis of the hippocampus reflects this distinction where anterior parts might support processing of items with lower spatial resolution (26). Our data thus provide the first direct evidence that hippocampal firing predicts workload, which matches the trial difficulty as derived from the subjects’ performance (Fig. 1, B and C).

In previous research on WM maintenance, frontoparietal cortical regions provided the most important nodes of the circuit (31). The MTL was long considered unimportant for WM because patients with MTL resection had apparently normal short-term memory. This finding implied that episodic memory, but not WM, is hippocampus dependent (32). However, MTL activity during WM tasks has been demonstrated in neuroimaging research (29) and neurophysiology research (3, 15). We extend these findings by reporting a workload-dependent increase in neural firing across our pool of neurons and the subsequent decoding of subject behavior.

Retrieval. Last, during retrieval, the highest proportion of neurons reacting to the probe letter was identified in the entorhinal cortex, whereas neurons in all three of our MTL areas predicted the IN or OUT responses of subjects (Fig. 4F). This finding supports a contribution of the MTL in stimulus-response mapping. The differential firing (Fig. 4F) for letters in memory, compared to letters that were not in memory, suggests a sensitivity of memory content contingent on the task demands, as all letters had been presented multiple times in the task, whereas their relevance for memory varied from trial to trial.

Together, hippocampal neuron population firing was content specific for letter numerosity, showed persistent stimulus-selective delay activity, and predicted a memory-dependent behavioral response (right/wrong; in/out). The hippocampus is thus an active node of the verbal WM network.

Functional connectivity by gated PLV. How does the hippocampus communicate within the WM network (21)? Maintenance consistently elicited enhanced PLV in a specific band in the theta-alpha frequency range in all but one subject (Fig. 6). Similarly, the spike-field coherence between neuronal action potentials and LFP was higher during maintenance than during encoding (Fig. 5J and fig. S6). The PLV had a well-defined onset in the last 2 s of the maintenance period and rapidly diminished following the probe letter. The PLV increased with workload. To explain a similar finding in iEEG power at selected cortical electrode sites, the term “gating” has been introduced (8). In our task design, this would indicate that neural communication between brain regions is channeled through theta-alpha coupling within a limited time window between the onset and offset of the maintenance phase during trials with a high workload.

The PLV synchronization occurred strongly with occipitoparietal scalp electrodes, a common locus for alpha waves. This finding is consistent with scalp EEG findings in several healthy subjects recorded while performing the same task (5, 11) and in similar studies (6). However, the hypothesis inspired by rodent research—that hippocampal theta would drive human frontal midline theta in the scalp EEG (5, 16)—is not conclusively supported by our results. Our finding of significant PLV modulation in the theta-alpha band extends the results of a recent study (22) that focused uniquely on the 3- to 8-Hz theta band for widespread synchronization. The enhanced hippocampal coupling between neuronal firing and LFP during maintenance suggests that neuronal firing may be integrated into the WM network via coupling to LFP and long-range connectivity in the alpha-theta band. Thus, the load dependence and the temporal alignment with neural activity at the population level and the single-neuron level further suggest that the PLV directly gates neural activity within the hippocampal circuits that participate in WM.

Communication through coherence

Given that we have demonstrated the plausibility of our results, we now suggest how they integrate into the current knowledge about WM mechanisms. The integration of distributed components of information across brain regions has been proposed to rely on long-range recurrent connections that support oscillatory signals (20), for example, during WM tasks (10, 12). Regarding the frequency of oscillation, a relationship between the rhythm frequency and the spatial extent of engaged brain networks has been proposed, with low frequencies binding large-scale networks and high frequencies coordinating smaller networks (22). There are several examples in which neural oscillations play a role in coordinating functional neuronal assemblies thought to be responsible for computation and communication in large-scale brain networks (10, 33). Thus, the interpretation that action potentials generated in the hippocampus are embedded within the WM network through coupling with LFP and long-range recurrent connections seems consistent with the current view on neural communication within and between brain regions.

Hippocampus in the WM network

In the task, the visual perception of the Roman characters was transformed to an abstract representation with reduced complexity. This verbal representation was maintained over the delay period in the phonological loop, with the aim of producing an appropriate behavioral response (1). The preferred coupling from hippocampus to frontal and parietal scalp electrodes agrees with the cortical sites where letter stimuli elicited persistent activity (2).

Neural activity patterns reflected differences between low and high load. They leveled off at high load and did not incrementally increase with set size. This finding may reflect the processing capacity limits (29), which are more closely related to the concept of workload. It further suggests that the hippocampal activity pattern decoded the quantity of letters that could actively be maintained in memory.

The separation of WM trial periods was expected, as single-neuron data also showed period-dependent differences in WM activity (3, 4) during object maintenance and as fMRI studies showed that specific parts of the hippocampus were involved in maintenance (34), encoding, and recognition. It seems plausible that population neuron activity during encoding and recognition shares processes with long-term memory; however, maintenance neurons are specific for WM. Furthermore, the time course of activity (Fig. 5) suggests that the WM network’s function goes beyond memorization toward a task-related preparation in expectance of the probe.


We provide direct evidence for hippocampal involvement in processing workload during the maintenance of information in WM. We observed enhanced neural activity at high workload that appeared as both increased persistent firing of hippocampal neurons and increased functional interaction between hippocampus and cortex through synchronized theta-alpha EEG oscillations.



We used a modified Sternberg task in which the encoding of memory items, maintenance, and recall were temporally separated (Fig. 1A) (5, 6). Each trial started with a fixation period ([−6, −5] s), followed by the stimulus ([−5, −3] s). The stimulus consisted of a set of eight consonants at the center of the screen. The middle four, six, or eight letters were the memory items, which determined the set size for the trial (4, 6, or 8, respectively). The outer positions were filled with “X,” which was never a memory item. After the stimulus, the letters disappeared from the screen, and the maintenance interval started ([−3, 0] s). A fixation square was shown throughout fixation, encoding, and maintenance. After maintenance, a probe was presented. The subjects responded with a button press to indicate whether the probe was part of the stimulus. The subjects were instructed to respond as rapidly as possible without making errors. The hand used for the response was counterbalanced across subjects within the clinical constraints. After the response, the probe was turned off, and the subjects received acoustic feedback regarding whether their response was correct or incorrect. Before initiating the next trial, the subjects were encouraged to blink and relax. The subjects performed 50 trials in one session, which lasted approximately 10 min. Trials with different set sizes were presented in a random order, with the single exception that a trial with an incorrect response was always followed by a trial with a set size of 4. The task is freely available at During the recording period of several days, several subjects decided to perform more than one session of the task up to a total of seven sessions.


Nine subjects participated in the study (table S1). All subjects were implanted with depth electrodes in the MTL for the potential surgical treatment of epilepsy. All subjects provided written informed consent for the study, which was approved by the institutional ethics review board (PB 2016-02055). All subjects had normal or corrected-to-normal vision and were right handed as confirmed by neuropsychological testing. The depth electrodes (1.3 mm diameter, 8 contacts of 1.6 mm length, and spacing between contact centers 5 mm; Ad-Tech, Racine, WI, were stereotactically implanted into the amygdala, hippocampus, and entorhinal cortex. Each macroelectrode had nine microelectrodes that protruded approximately 4 mm from its tip. In addition, scalp EEG electrodes were placed at the sites of the 10-20 system with minor adaptations to avoid surgical scalp lesions.

Depth electrode localization

Electrodes were localized using postimplantation computed tomography (CT) scans and postimplantation structural T1-weighted MRI scans. For each subject, the CT scan was registered to the postimplantation scan as implemented in FieldTrip (35, 36). In the coregistered CT-MR images, the electrode contacts were visually marked. The contact positions were normalized to the MNI space and assigned to a brain region using the Brainnetome Atlas (37). In addition, depth electrode positions were verified by the neurosurgeon (L.S.) after merging preoperative MRI with postimplantation CT images of each individual subject in the plane along the electrode (iPlan Stereotaxy 3.0, Brainlab, München, Germany). We grouped electrodes according to whether they were recorded from the SOZ or another area. As an illustration, we projected the three-dimensional positions of all electrode tips, i.e., the positions of the microelectrodes, on a parasagittal plane (MNI space x = −25.2 mm; Fig. 1D).

Recording setup

Intracranial data were acquired at sampling frequencies of 4 kHz for the contacts of the macroelectrodes and 32 kHz for the microelectrodes via the ATLAS recording system (0.5- to 5000-Hz passband, Neuralynx, Bozeman, MT, USA; iEEG was recorded (Neuralynx ATLAS; sampling rate, 4000 Hz; passband, 0.5 to 1000 Hz) against a common intracranial reference and subsequently transformed to a bipolar montage. Only one bipolar signal from a pair involving the deepest contacts of each macroelectrode was used for further analysis. Scalp EEG was recorded (NicOne; sampling rate, 256 Hz; passband, 0.3 to 100 Hz). Scalp EEG was re-referenced to the averaged mastoid channels.

Spike detection and single-unit identification

The Combinato package 17 ( was used for spike sorting (38). Combinato follows the same procedure of other freely available software packages: peak detection in the high-pass (>500 Hz) signal, computation of wavelet coefficients for detected peaks, and superparamagnetic clustering in the feature space of wavelet coefficients. As an advantage over other clustering procedures, Combinato provides better automated artifact rejection and is more sensitive in the detection of clusters of small size (few action potentials). We visually inspected each identified cluster based on the shape and amplitude of the action potentials and the interspike interval (ISI) distributions. We removed clusters that exhibited a low firing rate (<0.1 Hz), noisy waveforms, or nonuniform amplitude or shape of the action potentials in the recorded time interval. Moreover, to avoid overclustering, we merged highly similar clusters identified on the same microelectrode. Last, we computed several metrics of spike sorting quality (fig. S2): The percentage of ISI below 3 ms was 1.58 ± 4.18%, and the ratio of the peak amplitude of the mean waveform to the SD of the noise was 3.28 ± 1.56.


For statistical significance analysis, we used nonparametric permutation tests. To assess the significance of a value, we created a null distribution estimated from n > 200 permutations on data with scrambled labels. The P value is limited by the number of permutations as P = 1/(number of permutations + 1). Reported P values were based on the empirically estimated null distribution. If no value of the null distribution exceeded the observed value, we reported the exact P value.

To compare two distributions, we used a permutation t test wherein we randomly mixed values from the two distributions 2000 times to create a distribution of t values and computed an empirical P value from this distribution using the EEGLAB toolbox (39). For comparisons between more than two groups, we used F statistics. When comparing firing rates or power spectra, we used cluster-based nonparametric permutation tests (36). For each test, matrix sizes, number of samples, and features are provided in table S2.

dPCA of firing rates

To illustrate the firing rates of all neurons over the periods of a trial, we used dPCA as a dimensionality reduction method (24). In contrast to PCA, which explains the variance of components regardless of task conditions, dPCA incorporates the task conditions, i.e., set size in our task, when explaining the variance. We used nonoverlapping bins of 2 ms for neuronal firing and smoothed the rates in a window of 1000 ms with a Gaussian kernel. We z-scored the resulting rates based on the mean of the firing rate during fixation. We computed time-dependent and set size–dependent dPCA components and ordered them by explained variance. To test the significance of the explained variance, we created a null distribution of explained variance by scrambling labels and computing dPCA 200 times. For each component, we compared the explained variance to the null distribution.

For the component trajectories in the dPC space, we defined multidimensional distanced(p(t),q(t))=i=1n(pi(t)qi(t))2(1)and multidimensional speedV(t)=1n×i=1n(pi(t)pi(tΔt))2Δt(2)where p(t) and q(t) are neural activity vectors, Δt = 20 ms is the step size, and n is the number of dimensions in the component space. To create a distribution for multidimensional speed in different trial phases, we ran dPCA analysis 200 times, subsampling (20%, with replacement) neural firing rate trajectories for each run.

Neural decoding of firing rates

We pooled the firing rates of all neurons recorded across sessions and subjects to form a pseudopopulation of neurons. In this neural state space, we trained a decoder to perform population decoding for the prediction of correct/incorrect trials, set size of trials, and IN/OUT trials. The input was the matrix of neural firing rates of each trial of all neurons, binned in 1-ms windows, and smoothed in a 1500-ms window. The neurons were grouped according to anatomical location (hippocampus, entorhinal cortex, and amygdala). As the classifier, the decoder was trained and tested at 250-ms steps using a support vector machine as implemented in the NDT toolbox (25). Independently for each time point, classification was performed by assigning firing rates from randomly selected trials for each neuron to training and test splits in each of the 500 resampling runs. We first classified the trials based on workload (set size 4 versus set sizes 6 and 8) using 10 cross-validation splits. We selected units from sessions with at least five incorrect responses and then classified correct/incorrect trials using five cross-validation splits. To classify the IN/OUT trials after the probe, we used a narrow window size of 500 ms and 10 cross-validation splits. We tested significance by creating a null distribution of accuracies with 500 runs of the decoder with scrambled labels and comparing the true accuracy to the null distribution. This approach resulted in a minimum value of P = 0.002 for the statistical significance. We corrected the P value for multiple comparisons using Bonferroni’s method. For the 12 tests performed (Fig. 4, D to F), the corrected value was P = 0.002 × 12 = 0.024 and remained statistically significant.

EEG preprocessing

Data were visually inspected, and trials that contained eye movements or artifacts were removed using the toolbox EEGLAB (39).

Frequency spectra

To calculate spectra during maintenance (Fig. 5, B, C, E, F, and H), the data were segmented into 2-s epochs and grouped according to set size. We focused on the last 2 s of the maintenance period to avoid interference from visual evoked responses and eye blink artifacts. We used the multitaper method (±2 Hz smoothing in the frequency domain with seven tapers) as implemented in FieldTrip (36) to reduce spectral leakage and control the frequency smoothing (36).

Phase-locking value

To probe the functional connectivity between the cortex and hippocampus, we calculated the PLV between the scalp EEG and iEEG channels (21, 36).PLVi,j(f)=1N|n=1NXi(f)(Xj(f))*|Xi(f)||Xj(f)||(3)where PLVi,j is the PLV between channels i and j, N is the number of trials, X(f) is the Fourier transform of x(t), and (∙)* represents the complex conjugate.

Using the spectra of the 2-s epochs, phase differences were calculated for each electrode pair (i, j) to quantify the interelectrode low-frequency phase coupling. The phase difference between the two signals indexes the coherence between each electrode pair and is expressed as the PLV. The PLV ranges between 0 and 1, with values approaching 1 if the two signals show a constant phase relationship over all trials (21).

To test the statistical significance of the PLV, we first permuted trials 200 times and computed the 95th percentile threshold within trials of set size 8. We used this one-tailed test for increase in the PLV at each frequency bin independently. Only the electrode pairs with PLV above threshold were selected for further processing. From the PLV spectra, we selected the electrode pair and frequency range of the highest PLV for each subject (table S1). Further analysis focused on these electrode pairs.

The workload-dependent increases in the PLV for the selected electrode pairs were subsequently assessed by subtracting the PLV for trials with set size 4 from the PLV for set size 8. The statistical significance was then estimated using a permutation test, in which a null distribution was created by randomly assigning trials (i.e., 2-s epochs) into two set size conditions, computing the PLV differences between conditions, and repeating this procedure 200 times. From the frequency spectra of the PLV differences for each subject, we considered only the frequency band of the PLV difference that was above the 95th percentile to be statistically significant (table S1). The edges of this individual frequency band (black bar in Fig. 5E) are presented in table S1, and further analysis focused on this individual frequency band.

Next, we determined the PLV difference between the two behavioral conditions encoding and maintenance by subtracting the PLV during encoding from the PLV during maintenance. Similar to the results above, the statistical significance was estimated using a permutation test, in which a null distribution was created by randomly assigning 2-s epochs into two behavioral conditions, computing the PLV differences between conditions, and repeating this procedure 200 times. Only the frequency points with PLV above the 95th percentile threshold were considered statistically significant. Thus, we could relate the PLV to the behavioral condition during task performance.

Time-frequency analysis

To calculate spectral power in the time-frequency domain (Fig. 5, A and G), each trial was transformed to a time-frequency map. We used multitapers with a window width of 10 cycles per frequency point, smoothed with 0.2 × frequency, and used three tapers (36). We computed power in the frequency range of 4 to 30 Hz with a time resolution of 0.1 s. The power during fixation ([−6.0, −5.0] s) served as a baseline for the baseline correction (power − powerbaseline)/powerbaseline for each time-frequency point. Using the same parameters, we computed the PLV in the time-frequency domain (Fig. 5D).

Spike-field coherence

We next quantified the coupling between the neuronal action potentials and the surrounding field. We calculated the spike-field coherence between a neuronal spike and the LFP recorded from the same microelectrode. For this step, we used the pairwise phase consistency (PPC; P^2) (21, 36, 40). Compared to PLV, this measure is more balanced toward differences in spike counts. To test the statistical significance of the PPC, we created a null distribution by permuting trials for spikes 200 times and computed the 95th percentile threshold. Using this threshold at each frequency bin, we computed the percentage of neurons with significant PPC (fig. S6).


Supplementary material for this article is available at

Fig. S1. Metrics for spike sorting for the identification of putative neuronal units.

Fig. S2. Examples of maintenance neurons.

Fig. S3. Details of the dimensionality reduction technique (dPCA) used for state space analysis.

Fig. S4. Frontal scalp EEG and PLV in subject 1.

Fig. S5. Scalp EEG power spectra for all subjects.

Fig. S6. Details of the spike-field coherence comparison between maintenance and encoding.

Table S1. Subject characteristics.

Table S2. Statistics reporting checklist.

Movie S1. Attractor dynamics in neuronal dPCA state space throughout the task.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank the physicians and the staff at Schweizerische Epilepsie-Klinik for their assistance and the patients for their participation. Funding: We acknowledge grants awarded by the Swiss National Science Foundation (SNSF 320030_176222 to J.S.), Mach-Gaensslen Stiftung (to J.S.), Stiftung für wissenschaftliche Forschung an der Universität Zürich (to J.S.), and Forschungskredit der Universität Zürich (to T.F.). The funders had no role in the design or analysis of the study. Ethics statement: The study was approved by the IRB, and all subjects gave written informed consent. Author contributions: T.F. and J.S. designed the experiments. P.H. set up the recordings. T.F. and J.S. performed the experiments. E.B., T.F., and J.S. analyzed the data. T.G. provided patient care. L.S. performed surgery. E.B, T.F., P.K., and J.S. wrote the manuscript. All authors discussed the results and commented on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article