Ripples reflect a spectrum of synchronous spiking activity in human anterior temporal lobe

We examined intracranial EEG (iEEG) recordings in 21 participants with intracranial electrodes placed for seizure monitoring as they performed a verbal episodic memory task (Figure 1A; see Materials and methods). In an example electrode in the medial temporal lobe, we observed transient increases in high-frequency activity (70–200 Hz; HFA) in individual trials immediately before participants vocalize their response during the retrieval period (Figure 1C). When averaging across all trials, the increases in HFA prior to vocalization appear sustained, consistent with previous studies of episodic memory retrieval (Burke et al., 2014; Yaffe et al., 2014).

Figure 1 with 6 supplements see all

Download asset Open asset

Figure 1—source code 1

Matlab code of ripple events in the iEEG signal.

https://cdn.elifesciences.org/articles/68401/elife-68401-fig1-code1-v2.zip

Download elife-68401-fig1-code1-v2.zip

We hypothesized that the transient increases in HFA observed in individual trials may be related to discrete 80–120 Hz ripples that have been previously associated with human memory retrieval (Vaz et al., 2019). We therefore identified ripples in each iEEG electrode in each participant (ripple rate .35 ± .04 Hz (mean ± SEM) across electrodes across all participants; Figure 1B; Figure 1—figure supplement 1 to Figure 1—figure supplement 5; see Materials and methods). In the same example electrode, the transient increases in HFA observed in each trial correspond to the detection of individual ripples (Figure 1C). Across all trials, ripple rates demonstrate a clear increase that coincides with the sustained increase in HFA. We examined whether the changes in HFA and ripple rate were similarly modulated by successful memory retrieval in this examplar electrode (Figure 1D). Both HFA and ripple rates increase prior to vocalization during successful memory retrieval trials compared to trials in which the participant failed to successfully retrieve the correct word.

Since ripple characteristics can vary significantly across brain regions, we examined whether the differences in HFA and ripple rate between correct and incorrect retrieval trials exhibit a similar spatial pattern across the brain. Across participants, HFA increases during the one second prior to vocalization are greater during correct memory retrieval compared to incorrect memory retrieval in several anatomic regions (Figure 1E; see Materials and methods). When examining ripple rates, we also found significant increases during correct compared to incorrect trials in similar anatomic regions during this same time period (Figure 1E; Figure 1—figure supplement 6B). Across ROIs from the entire cortex, the participant average difference in HFA between correct and incorrect trials is positively correlated with the difference in ripple rate (r = .13, p = 6.1 × 10^-4; Figure 1F). Within two specific brain regions, the medial temporal lobe and the anterior temporal lobe, this positive correlation between the participant average difference in HFA between correct and incorrect trials and the average difference in the ripple rate is greater (medial temporal lobe, MTL: r = .60, p = .00072; anterior temporal lobe, ATL: r = .23, p = .026; Figure 1G; Figure 1—figure supplement 6C,D). Together, these data suggest that the changes in HFA and ripple rate observed with successful memory retrieval obey a similar anatomic distribution.

We then examined this relation between HFA and 80–120 Hz ripples within retrieval trials in all individual electrodes in all participants. In each electrode, we computed the Pearson correlation between the average HFA and ripple rate across all retrieval trials, performed a Fisher’s $z$-transform to normalize the correlation coefficients across participants, and then computed an average across all electrodes in each participant. We similarly computed this correlation across random three second epochs throughout the experimental session, which we designated as baseline. In both cases, the distribution of correlations across participants was consistent and significantly greater than zero (baseline r = .085 ± .018; retrieval r = .093 ± .020; t(20) = 4.6, p < 1 × 10^-6, one-tailed t-test; Figure 1H). The relation between HFA and ripple rate was similar during memory retrieval and baseline (retrieval-baseline .0079 ± .0016, t(20) = 1.44, p = 0.083, paired t-test). However, the relation was stronger during correct retrieval as compared to incorrect retrieval (correct r = .316 ± .0339, incorrect r = .259 ± .0296; t(20) = 8.7, p < 2 × 10^-8; correct-incorrect .057 ± .004, t(20) = 2.40, p = .013, paired t-test; Figure 1I).

To explicitly examine the extent to which 80–120 Hz ripples contribute to 70–200 Hz power, we conducted a control analysis by removing the time indices in which ripples were detected from the iEEG trace and recomputed HFA power. While the distribution of 70–200 Hz power averaged across electrodes in each participant is significantly greater during correct compared to incorrect memory retrieval (correct-incorrect .0051 ± .0049, t(20) = 2.39, p = .013), removing the discrete 80–120 Hz ripples eliminates this difference (correct-incorrect: .0075 ± .0053, t(20) = 1.42, p = .086; Figure 1J; Figure 1—figure supplement 6E). HFA power during correct retrieval is significantly reduced when removing the 80–120 Hz ripples (original correct-control correct: .0096 ± .0023, t(20) = 4.21, p = 2.2 × 10^-4). We also examined the correlation across all ROIs between the participant average difference in ripple rate between correct and incorrect trials and the difference in HFA after removal of the discrete ripples and found that this relation is no longer significant (r = -.042, p = .270; Figure 1—figure supplement 6F,G). Finally, to confirm that much of the 70–200 Hz power is driven by relatively band limited 80–120 Hz ripples, we repeated our analysis after detecting discrete ripple events in a higher 120–200 Hz frequency band. Across ROIs from the entire cortex, we did not find a significant correlation between the participant average difference in HFA between correct and incorrect trials and the difference in 120–200 Hz ripple rate (Figure 1—figure supplement 6H).

In a subset of six participants, we had the opportunity to record micro-scale local field potentials (LFPs) and single unit spiking activity from a micro-electrode array (MEA) implanted in the anterior lateral temporal lobe underneath the iEEG electrodes (Figure 2A; Figure 2—figure supplement 1; see Materials and methods). In an example participant, ripples present in a single iEEG recording electrode overlying the MEA clearly occur simultaneously with ripples in the LFP across multiple micro-electrodes within the MEA (Figure 2B; Figure 2—figure supplement 2). These discrete events are accompanied by increases in spiking activity across multiple units, and therefore transient increases in the overall population spiking rate across the MEA (Figure 2—figure supplement 3 and Figure 2—figure supplement 4).

Figure 2 with 5 supplements see all

Download asset Open asset

Figure 2—source code 1

Matlab code of correlations between continuous spiking, LFP and iEEG.

https://cdn.elifesciences.org/articles/68401/elife-68401-fig2-code1-v2.zip

Download elife-68401-fig2-code1-v2.zip

The continuous time data of iEEG, LFP, and spiking activity suggest that 80–120 Hz ripples at both the macro- and micro- spatial scale and single unit spiking activity are related. To examine this relation, we computed $z$-scored 80–120 Hz ripple band amplitude in both the overlaying iEEG electrode and the average $z$-scored ripple band amplitude and spike rate in each of the MEA electrodes during 100 ms non-overlapping windows over all retrieval trials (ripple rate .84 ± .43 Hz across micro-electrodes across all participants). Across all time windows in this participant, the average spike rate across MEA electrodes is significantly correlated with iEEG and LFP ripple band amplitude (spike rate v LFP amplitude r = .61, p < 1 × 10^-18; spike rate v iEEG amplitude r = .12, Pearson correlation; Figure 2C). We found that the relation between spiking activity and ripple band amplitude at both spatial scales is consistent and significant across participants (spike rate v LFP amplitude, Fisher z-transform: r = .751 ± .188; t(5) = 4.00, one-tailed t-test; spike rate v iEEG amplitude: r = .118 ± .049; t(5) = 2.39, p = .031; Figure 2D).

These data demonstrate that the continuous time measure of 80–120 Hz ripple band amplitude is related to spiking activity. However, we were interested in understanding whether the amplitude and duration of discrete ripples may exist on a continuum reflecting underlying neuronal activity. We therefore relaxed our criteria for identifying discrete ripple events in order to detect ripples that are smaller and shorter duration, which are often assumed to be noise (see Materials and methods). We found ripples at both the macro- and micro-scale with a range of amplitudes and durations (Figure 1—figure supplement 2). During every ripple detected in each LFP trace, we collected the average spike rate of units recorded in the respective MEA electrode and computed the Pearson correlation between LFP ripple amplitude and spike rate across all ripples. Across participants, LFP ripple amplitude is consistently and significantly correlated with spike rate (Fisher $z$-transform, r = .10 ± .02; t(5) = 3.62, p = .008, one-tailed $t$-test; Figure 2E; Figure 2—figure supplement 5A). Even when ripples have amplitudes or durations below previously used thresholds, spiking activity is present in the microelectrode recording (Figure 2—figure supplement 5B).

While we found a strong relation between spiking activity and ripple amplitude, this observation could be confounded by a correlation between ripple amplitude and duration (Figure 1—figure supplement 2E). Larger amplitude are larger in duration and therefore may have more opportunity to co-occur with spikes by chance. To account for this, we shuffled the time indices of the detected spike times and computed the correlation between LFP ripple amplitude and the spike rate. The true relation between LFP ripple band amplitude and the spike rate is significantly greater than the shuffled distribution (true-shuffled r = .096 ± .020; t(5) = 3.312, p = .011, paired one-tailed t-test; Figure 2F).

In a similar manner, during every iEEG ripple we determined how many individual units spike as a proportion of the total number of units identified in each experimental session (Figure 2—figure supplement 5C). We computed the Pearson correlation between the percentage of actively firing units and the iEEG ripple band amplitude and duration across all detected iEEG ripples in every participant (Figure 2H; Figure 2—figure supplement 5D). Across participants, this correlation is significant (Fisher $z$-transform, r = .15 ± .03; t(5) = 4.679, p = .003, one-tailed $t$-test). We performed the same shuffling procedure to account for ripple duration, and found that the true correlations across participants are significantly larger than the shuffled data (true-shuffled r = .069 ± .027; t(5) = 2.517, p = .027, paired one-tailed t-test; Figure 2H). Together, these data demonstrate a strong relation between underlying unit spiking activity and ripples observed at the micro- and macro-scale in the human temporal cortex.

Our data suggest that ripples observed at both spatial scales may be related to one another. We hypothesized that the amplitude of the ripple observed in the iEEG signal is related to both the total number of LFP ripples and the extent to which the LFP ripples observed across the underlying MEA electrodes are aligned (Figure 3A).

Figure 3 with 3 supplements see all

Download asset Open asset

Figure 3—source code 1

Matlab code of pairwise phase consistency between LFP ripple signal and iEEG ripple amplitude.

https://cdn.elifesciences.org/articles/68401/elife-68401-fig3-code1-v2.zip

Download elife-68401-fig3-code1-v2.zip

We first examined the relation between the amplitude of the iEEG ripple and the number of LFP ripples simultaneously present in the MEA electrodes. In every participant, we detected ripples in each of the four iEEG electrodes closest to the MEA. For every detected ripple, we computed the mean 80–120 Hz ripple band amplitude across the iEEG electrodes and the number of LFP ripples simultaneously observed across the MEA (Figure 3—figure supplement 1A). In each participant, the iEEG ripple amplitude is positively correlated with the percent of MEA electrodes exhibiting LFP ripples across all detected iEEG ripples (% MEA electrodes with ripples, Fisher z-transform; r = .11 ± .02, p = .005; Figure 3B; Figure 3—figure supplement 1B,C). We accounted for the possibility that the longer durations observed in higher amplitude iEEG ripples may result in a larger number of detected LFP ripples by using a similar shuffling procedure. In this case, during each shuffle we performed a random circular shift of the time indices of the detected LFP ripples. After accounting for these longer durations, we still found that the true relation between iEEG ripple amplitude and the number of simultaneously detected LFP ripples is significantly greater than the shuffled distribution (true-shuffled r = .05 ± .02; t(5) = 2.543, p = .026, paired one-tailed t-test; Figure 3B).

We then examined the relation between the amplitude of the iEEG ripple and the extent to which the LFP ripples in the underlying MEA are synchronized. For every detected iEEG ripple, we extracted the LFP 80–120 Hz ripple band instantaneous phase for all 96 MEA electrodes and computed the maximum pairwise phase consistency (PPC) over all time points within the duration of that iEEG ripple (Figure 3C and D; see Materials and methods). Across participants, the PPC is significantly correlated with the maximum amplitude of the observed iEEG ripples (Fisher z-transform; r = 1.06 ± .43; t(5) = 2.34, p = .033, one-tail t-test; Figure 3E–F; Figure 3—figure supplement 1D). In addition, the correlations across participants are significantly greater than those that would be observed by chance (true-shuffled r = .03 ± .01; t(5) = 2.28, p = .036, paired one-tail t-test; Figure 3E; see Materials and methods). We repeated this analysis using only microelectrodes with detected ripples and found that, across participants, the PPC is still significantly correlated with the maximum amplitude of the observed iEEG ripples (Fisher z-transform; r = .066 ± .02; t(5) = 3.71, p = .014, one-tail t-test). These data together suggest that the iEEG ripple reflects both the aggregate sum and alignment of the underlying LFP ripples.

To further examine the relation between ripples detected in the LFP signal and ripples detected in the iEEG explicitly, we measured the coincidence of ripples detected at the two spatial scales by computing the cross-correlogram of ripples detected in the LFP and iEEG traces. We found that ripples are coincident above chance for all detection parameters tested (Figure 3—figure supplement 2; see Materials and methods). Moreover, the extent to which ripples are coincident between electrodes at the two spatial scales is significantly related to the distance between them (all participants: r = -.572, p = .0035; across each participant, n = 6: r = -.595 ± .271, mean ± SEM, t(5) = -2.40, p = .0615; Figure 3G; Figure 3—figure supplement 3).

Given that synchronization of LFP ripples was associated with increased ripple amplitudes, we also hypothesized that spike timing during ripples would be phase locked (Le Van Quyen et al., 2008). In individual participants, we often observed that unit firing preferentially occurs at the trough of the corresponding LFP ripples (Figure 4A). When the LFP ripples are aligned, spiking activity also appears to preferentially occur at the trough of the overlying iEEG ripple. In all participants, spikes from all units are locked to the trough of the 80–120 Hz ripple band in the LFP signal (p < 10^-4, Rayleigh test across all units in each participant; p = 4.8 × 10^-4, Rayleigh test across six complex means, one from each participant; Figure 4B). We did not observe such phase consistency when examining the extent to which spiking activity is locked to the phase of 80–120 Hz ripple band activity in the macro-scale iEEG signal (p = .12, Rayleigh test across all units and across participants; Figure 4B).

Figure 4 with 2 supplements see all

Download asset Open asset

Figure 4—source code 1

Matlab code of pairwise phase consistency between spiking and LFP.

https://cdn.elifesciences.org/articles/68401/elife-68401-fig4-code1-v2.zip

Download elife-68401-fig4-code1-v2.zip

When we visualized the spike triggered average of the LFP signal in individual participants, we often observed that spiking activity also appeared locked to negative deflections in the LFP (Figure 4C). These negative deflections contain spectral power within low frequencies. We therefore also examined the distribution of 2–10 Hz low-frequency phases present in the LFP signal around each spike and found significant locking to the trough in all participants (p < 10^-4, Rayleigh test; p = 5.4 × 10^-4, Rayleigh test across six complex means, one from each participant; Figure 4D). Spikes from all units appear locked around the trough of the 2–10 Hz low-frequency iEEG signal when pooled, which reflects the negative deflection in the iEEG signal, but when examined separately for each participant, the apparent spike locking to the 2–10 Hz iEEG signal is not consistently at the same phase across participants (p < 10^-4, Rayleigh test across all units in each participant; p = .70, Rayleigh test across six complex means, one from each participant; Figure 4D).

To examine the relation between spiking activity and individual frequencies within the LFP signal, we computed PPC across all spikes within each MEA electrode for each frequency between 2 Hz and 400 Hz (see Materials and methods) (Vinck et al., 2010). Across participants, spiking activity is significantly locked to specific high-frequency bands in the LFP (peak 86.9 Hz, p < .05, permutation test; see Materials and methods Figure 4E). We confirmed that spiking activity is locked to this high-frequency band across participants by also computing the phase-locking value (Figure 4—figure supplement 1A). This observed locking between spikes and this high-frequency band in the LFP signal was robust to different detection thresholds (Figure 4—figure supplement 2). Spikes also appear locked to a low-frequency band, but this likely represents the sharp negative deflections observed in the iEEG and LFP traces that accompany burst of spiking activity. Spikes are significantly more locked to high frequencies when they arise during ripples as compared to between ripples (p < .05, permutation test; Figure 4, Figure 4—figure supplement 1B).

We next examined the relation between the extent to which spiking activity locks to the 80–120 Hz frequency within each ripple and the amplitude of the ripple. Across participants, mean spike-LFP PPC within the 80–120 Hz ripple band is significantly correlated with 80–120 Hz ripple amplitude across all MEA electrodes (Fisher $z$-transform, r = .084 ± .026, t(5) = 3.27, p = .011; Figure 4G). To account for any possible effects of ripple duration on the calculation of PPC, we compared this true distribution to a chance distribution and found that across participants LFP ripple amplitude exhibits a significantly stronger correlation with spike locking to the 80–120 Hz band in the true data as compared to the chance distribution (t(5) = 2.64, p = .023; see Materials and methods). Together with our data examining the relation between spiking activity and ripple amplitude, these data suggest that the amplitude of ripples in the LFP signal may reflect both the sum and alignment of underlying spiking.

Finally, given the observed relation between spiking activity and ripples, we then examined whether ripples themselves also exhibit a phase preference. As with the individual spikes, we considered each LFP ripple as an event and visualized the ripple-triggered average of the iEEG and LFP signal (Figure 4—figure supplement 1C). Ripples appear to exhibit a clear relation with negative deflections in the iEEG and LFP trace. We therefore examined the distribution of 2–10 Hz low frequency phases present in the LFP signal during each LFP ripple and found significant locking to the trough in all participants (p < 10^-4, Rayleigh test across all ripples in each participant; p = .0099, Rayleigh test across complex means, one from each participant; Figure 4—figure supplement 1E). Micro-scale ripples are also locked to the 2–10 Hz low-frequency band in the iEEG signal within individual participants but at variable phases (p < 10^-4, Rayleigh test across all ripples in each participant; p = .72, Rayleigh test across six complex means, one from each participant).

Twenty-one participants with drug resistant epilepsy underwent a surgical procedure in which platinum recording contacts were implanted on the cortical surface as well as within the brain parenchyma. In each case, the clinical team determined the placement of the contacts to localize epileptogenic regions. In all the participants investigated here, the clinical region of investigation was the temporal lobes.

For research purposes, in six of these participants (4 female; 34.8 ± 4.7 years old) we placed one or two 96-channel microelectrode arrays (MEA; 4 × 4 mm, Cereplex I; Blackrock Microsystems, Inc, Salt Lake City, UT) in the anterior temporal lobe (ATL) in addition to the subdural contacts. We implanted MEAs only in participants with a presurgical evaluation indicating clear seizure localization in the temporal lobe and the implant site in the ATL was chosen to fall within the expected resection area. Each MEA was placed in an area of cortex that appeared normal both on the pre-operative MRI and on visual inspection. Across participants, MEAs were implanted 14.6 ± 3.7 mm away from the closest subdural electrode with any ictal or interictal activity identified by the clinical team. Four out of the six participants received a surgical resection which includes the tissue where the MEAs were implanted. One participant had evidence of focal cortical seizure activity and received a localized resection posterior to the MEA site. One participant did not have a sufficient number of seizures during the monitoring period to justify a subsequent resection. Neither participant experienced a change in seizure type or frequency following the procedure, or experienced any noted change in cognitive function. The data captured from these MEAs in these participants were included in a previous study (Vaz et al., 2020).

Data were collected at the Clinical Center at the National Institutes of Health (NIH; Bethesda, MD). The Institutional Review Board (IRB) approved the research protocol (11 N-0051), and informed consent was obtained from the participants and their guardians. All analyses were performed using custom built Matlab code (Natick, MA). Data are reported as mean ± SEM unless otherwise specified.

Each participant performed a paired associates verbal memory task (Yaffe et al., 2014; Jang et al., 2017; Vaz et al., 2020). Previous studies have demonstrated that correct memory retrieval in this task is associated with increases in high-frequency activity (Yaffe et al., 2014; Jang et al., 2017; Vaz et al., 2020). Here, we replicate these previous findings using a subset of participants that were included in these previous studies as well as additional new participants. During the study period, participants were sequentially shown a list of word pairs and instructed to remember the novel associations between each pair of words (encoding). Later during testing, they were cued with one word from each pair selected at random and were instructed to say the associated word into a microphone (retrieval).

A single experimental session for each participant consisted of 25 lists, where each list contained six pairs of common nouns shown on the center of a laptop screen. The number of pairs in a list was kept constant for each participant. Words were chosen at random and without replacement from a pool of high-frequency nouns and were presented sequentially and appearing in capital letters at the center of the screen. We separated the study and test of each word pair by a minimum lag of two study or test items. During the study period, each word pair was preceded by an orientation stimulus (’+’) that appeared on the screen for 250–300 ms followed by a blank interstimulus interval (ISI) between 500–750 ms. Word pairs were then presented stacked in the center of the screen for 4000 ms followed by a blank ISI of 1000 ms. Following the presentation of the list of word pairs, participants completed an arithmetic distractor task of the form A + B + C = ? for 20 seconds.

During the test period, one word was randomly chosen from each of the presented pairs and presented in random order, and the participant was asked to recall the other word from the pair by vocalizing a response. Each cue word was preceded by an orientation stimulus (a row of question marks) that appeared on the screen for 4000 ms followed by a blank ISI of 1000 ms. Participants could vocalize their response any time during the recall period after cue presentation. We manually designated each recorded response as correct, intrusion, or pass. A response was designated as pass when no vocalization was made, when the participants made an unintelligible vocalization like ’umm’, or when the participant vocalized the word ’pass’. During pass trials where no vocalization was present, we assigned a response time by randomly drawing from the distribution of correct response time during that experimental session. We did not include such pass trials where no vocalization was present in our analysis of incorrect trials. We defined all intrusion and other pass trials as incorrect trials. A single experimental session contained 150 total word pairs. Each participant completed between 1 and 3 sessions (2.2 ± .3 per participant). Participants studied 93 ± 8 word pairs, and successfully recalled 30.1% ± 4.1% of words. While patients were presented with 150 words pairs in each experimental session, the number of word pairs they actually studied was reduced if they did not complete the session due to interruptions or participant fatigue.

We collected intracranial EEG (iEEG) data from a total of 1660 subdural and depth recording contacts (79 ± 4 per participant; Figure 1—figure supplement 6). Subdural contacts were arranged in both grid and strip configurations with an inter-contact spacing of 10 mm. We captured iEEG signals sampled at 1000 Hz. For clinical visual inspection of the recording, signals were referenced to a common contact placed subcutaneously, on the scalp, or on the mastoid process. The recorded raw iEEG signals used for analyses were referenced to the system hardware reference, which was set by the recording amplifier (Nihon Kohden, Irvine CA) as the average of two intracranial electrode channels. We used the Chronux toolbox to apply a local detrending procedure to remove slow fluctuations (≤ 2 Hz) from the time series of each electrode and a regression-based approach to remove line noise at 60 Hz and 120 Hz (Mitra and Bokil, 2009). We did not see a noticeable peak at the 180 Hz harmonic when we surveyed the power spectral density of several electrodes for noise and therefore did not remove line noise at that harmonic to avoid introducing artifacts. We implemented additional thresholds to remove movement artifacts and pathological activity related to the patient’s epilepsy.

We quantified spectral power and phase in the iEEG signals by convolving the voltage time series with 200 linearly spaced complex valued Morlet wavelets between 2 and 200 Hz (wavelet number 6). We extracted data from all retrieval periods, beginning 4 s preceding vocalization to 1 s following vocalization and included a 1000 ms buffer on both sides of the clipped data. We squared and log-transformed the continuous-time wavelet transform to generate a continuous measure of instantaneous power for each frequency. To account for changes in power across experimental sessions, we $z$-scored power values separately for each frequency and for each session using the mean and standard deviation of all respective values for that session. When examining the average changes in high frequency activity (70–200 Hz) during memory retrieval across trials, we temporally smoothed the $z$-scored spectrogram for each iEEG channel using a sliding 600 ms window (90 % overlap) as a point of comparison with previous studies of human memory retrieval (Greenberg et al., 2015).

We localized electrodes in each participant by identifying high-intensity voxels in a post-operative CT image, which was co-registered to a pre-operative T1-weighted MRI. Electrode locations were adjusted to account for routine post-operative parenchymal shift by applying a standardized algorithm combining intraoperative photography, electrode spatial arrangement, and dural and pial surface reconstructions (Trotta et al., 2018). Pial surfaces were reconstructed using FreeSurfer (http://surfer.nmr.mgh.harvard.edu) (Fischl, 2012) and were resampled and standardized using the AFNI SUMA package (Cox, 1996). The resulting surfaces each contained 198812 vertices per hemisphere, with vertices indexed in a standardized fashion, such that for any vertex i, the i th vertex is located in an anatomically analogous location across participants. We identified the location of each MEA on each participant’s surface reconstruction. We co-registered the individual participant reconstructions with a standard template brain, and visualized the locations of each participant’s MEA on the template brain.

We aggregated vertices from the surface reconstruction into a standard set of surface-based regions of interest (ROIs) as previously described (Figure 1—figure supplement 6; Trotta et al., 2018). Briefly, we sampled 2400 equally spaced vertices per hemisphere to use as ROI centers. ROI centers were uniformly distributed across the surface at an average geodesic distance of approximately 5 mm. We assigned all vertices within a 10 mm geodesic radius of an ROI center to that ROI, which achieves a coverage of 99.9% coverage or greater of the pial surface in each participant (Trotta et al., 2018). Because ROIs overlap, vertices may be assigned to multiple ROIs. On average, there were 669.44 ± 74.30 vertices per ROI and each vertex mapped to 8.08 ± .90 ROIs. We modeled each electrode as a cylinder with radius 1.5 mm, found the pial vertices closest to it, and then assigned each electrode to the same ROIs as its nearest pial vertices. Due to the overlap between ROIs, each electrode is assigned to multiple ROIs and each ROI may contain more than one electrode. For analyses within ROIs across participants, we only included ROIs that contained electrodes from at least five participants.

We implemented several measures to provide the most conservative sampling of non-pathological signals possible. We implemented a previously reported automated trial and electrode rejection procedure based on excessive kurtosis or variance of iEEG signals to exclude high-frequency activity associated with epileptiform activity (Jang et al., 2017; Wittig et al., 2018; Vaz et al., 2019). We calculated and sorted the mean iEEG voltage across all trials, and divided the distribution into quartiles. We identified trial outliers by setting a threshold, Q3 + w*(Q3-Q1), where Q1 and Q3 are the mean voltage boundaries of the first and third quartiles, respectively. We empirically determined the weight w to be 2.3. We excluded all trials with mean voltage that exceeded this threshold. The average percent removed across all sessions in each participant due to either system-level noise or transient epileptiform activity was 5.17% ± .86% of all electrodes and 2.89% ± .34% of all trials.

In addition, system level line noise, eye-blink artifacts, sharp transients, and inter-ictal epileptic discharges (IEDs) can confound the interpretation of our results. We therefore implemented a previously reported automated event-level artifact rejection (Staresina et al., 2015; Vaz et al., 2019). We calculated a z-score for every iEEG time point based on the gradient (first derivative) and amplitude after applying a 250 Hz high pass filter (for identification of IEDs). All time points within 100 ms of any time point that exceeded a z-score of 5 with either gradient or high-frequency amplitude were marked as artifactual. We visually inspected the resulting iEEG traces and found that the automated procedure reliably removed IEDs and other artifacts. In total, following exclusion of electrodes because of artifact, we retained 1577 electrodes (75 ± 4 per participant) for analysis. We approximated a reference-free montage within each participant by subtracting the common average reference of all retained electrodes from the voltage trace of each individual electrode for that participant.

In six participants, we additionally captured spiking activity and micro-scale local field potentials (LFP) from the MEAs implanted in the anterior temporal lobe. Microelectrodes were arranged in a 10 × 10 grid with each electrode spaced 400 μm apart and extending 1.5 mm into the cortical surface (1.0 mm for one participant). Post-operative paraffin blocks of the resected tissue demonstrated that the electrodes extended approximately halfway into the 3-mm-thick gray matter. We digitally recorded microelectrode signals at 30 kHz using the Cereplex I and a Cerebus acquisition system (Blackrock Microsystems), with 16-bit precision and a range of ± 8 mV.

To extract unit spiking activity, we re-referenced each electrode’s signal offline by subtracting the mean signal of all the electrodes in the MEA, and then used a second order Butterworth filter to bandpass the signal between 0.3 and 3 kHz. Using a spike-sorting software package (Plexon Offline Sorter, Dallas, TX, USA), we identified spike waveforms by manually setting a negative or positive voltage threshold depending on the direction of putative action potentials. The voltage threshold was set to include noise signals used in calculating unit isolation quality (see below). Waveforms (duration, 1.067 ms; 32 samples per waveform) that crossed the voltage threshold were stored for spike sorting. Spike clusters were manually identified by viewing the first two principal components, and the difference in peak-to-trough voltage (voltage versus time) of the waveforms. We manually drew a boundary around clusters of waveforms that were differentiable from noise throughout the experimental session. In this manner, we identified a total of 989 putative single units across all sessions (average of 72 ± 21 units per participant). The average spike rate across all units was 2.82 ± .01 Hz. In addition to the spiking data, we also used a 500 Hz low pass filter to extract the LFP signals from each microelectrode, down-sampled to 1000 Hz, and then performed a similar line noise removal and channel selection procedure to that used for the iEEG channels to exclude artifacts related to epileptiform activity or other system level noise. Across the six participants, after pre-processing we retained recordings from 78 ± 27 MEA electrodes for further analysis.

Due to variability in the signal quality across recordings and the subjective nature of spike sorting, we quantified the quality of each unit by calculating an isolation score and signal to noise ratio (SNR) (Joshua et al., 2007). The isolation score quantifies the distance between the spike and noise clusters in a 32-dimensional space, where each dimension corresponds to a sample in the spike waveform. The spike cluster consisted of all waveforms that were classified as belonging to that unit, and the noise cluster consisted of all waveforms that crossed the threshold that were not classified as belonging to any unit. The isolation score is normalized to be between 0 and 1, and serves as a measure to compare the isolation quality of all units across all experimental sessions and participants. Across participants, the mean isolation score for all units was .93 ± .1.

In addition to isolation quality, we computed the SNR for each unit:

where V_peak and V_trough are the maximum and minimum voltage values of the mean waveform, and $C$ is a scaling factor (set as 5). To obtain $Noise$, we subtracted the mean waveform from each individual waveform for each identified unit, concatenated these waveform residuals, and then computed the standard deviation of this long vector. Therefore, the noise term quantifies the within-unit variability in waveform shape. Across participants, the mean SNR for all units was 1.71 ± .12.

We estimated the instantaneous spike rate for each unit by convolving the spike rasters with a Gaussian kernel (σ = 25 ms). We used the mean and standard deviation of the spike rate over an entire experimental session to generate a z-scored spike rate for each unit.

We detected ripples in both the iEEG and LFP signals as previously reported (Vaz et al., 2019). We first bandpass filtered the voltage time series in the ripple band (80–120 Hz) using a second order Butterworth filter, and then applied a Hilbert transform to extract the instantaneous amplitude and phase within that band. We selected events where the Hilbert envelope exceeded two standard deviations above the mean amplitude of the filtered traces. We only retained events that were at least 25 ms in duration and had a maximum amplitude greater than three standard deviations as ripples for analysis. We did not specify an upper limit for ripple duration. We joined adjacent ripples that were separated by less than 15 ms. We identified every ripple that satisfied these criteria in every electrode contact, and assigned each such identified ripple a start time index and an end time index. The difference between them defined the duration of each ripple.

To assess the overlap between detected ripples and inter-ictal epileptic discharge (IED) artifacts, we computed the joint probability of iEEG and LFP ripples and the identified IEDs for each participant. We found that IEDs overlapped with .79 ± .11% of iEEG ripples and with 1.38 ± .11% of LFP ripples across the six participants with MEAs (Figure 1—figure supplement 5A,B). We excluded all IEDs and high-frequency oscillations associated with IEDs (ripple on spike waveforms, pathologic ripples) and any detected ripple that overlapped with an IED from our analyses. The remaining ripples that we retained for our analyses therefore occurred without an associated IED and are more likely to be physiologic.

To examine the relation between ripple amplitude and spiking activity, as well as to examine the relation between ripples across spatial scales, we used the Hilbert phase and amplitude of the 80–120 Hz ripple band signal extracted from both the iEEG and LFP signals. To assess for a spectrum of ripple amplitudes and durations, we relaxed the detection thresholds to include all events during which the Hilbert amplitude of the LFP signal exceeded only one standard deviation above the mean amplitude of the filtered traces. We designated all such events with a minimum duration of 10 ms and with a maximum amplitude at least two standard deviations above the mean as putative ripples for these analyses.

To account for the possibility that ripples with higher amplitudes and therefore longer durations may be associated with more spiking activity by chance, we compared the true correlation between ripple amplitude and spiking activity to the correlations we would observe by chance. In each of 1000 permutations, we performed a random circular shift of the spike indices in each trial and computed the correlation between LFP ripple amplitude and spike rate across units and MEA electrodes. We compared the true correlation to the mean of the distribution of 1000 shuffled correlations in each participant. We determined that 1000 permutations was sufficient by initially examining the mean correlation as a function of the number of permutations in a single participant, and found that the mean value for the correlation observed by chance converged after only 500 permutations. We performed a similar permutation procedure when examining the relation between iEEG ripple amplitude and the proportion of active units, and between iEEG ripple amplitude and the number of underlying LFP ripples.

To examine the extent to which individual events such as spikes or ripples are aligned to consistent phases in the LFP or iEEG oscillations, we computed the pairwise phase consistency (PPC) (Vinck et al., 2010). Briefly, for each spike or ripple, we extracted the instantaneous phase of the LFP or iEEG signal either of individual frequencies or within low (2–10 Hz) or ripple band (80–120 Hz) frequency bands. For individual frequencies, we used the instantaneous phase extracted by convolving the LFP or iEEG time series with complex valued Morlet wavelets (wavelet number 6) for 60 frequencies logarithmically spaced between 2 and 400 Hz. To extract the instantaneous phase of the two frequency bands, 2–10 Hz and 80–120 Hz, we filtered the LFP and iEEG signal into each frequency band and then extracted the instantaneous phase from the complex time series generated by the Hilbert transform of the filtered time series. Across multiple spikes or ripples, we therefore generate a distribution of phases. To calculate the PPC, we computed the average angular distance, or vector dot product, for all pairs of phases in each distribution. We defined the preferred phase for each distribution as the phase angle of the complex mean of the distributions of these phases. In addition to PPC, we also assessed phase consistency by testing whether each distribution of phases significantly deviated from a uniform distribution using a Rayleigh test of uniformity.

We used PPC to examine the extent to which 80–120 Hz ripple band phases are aligned across all microelectrodes in each MEA during each ripple detected in the larger scale iEEG signal. In this case, during every time point within each iEEG ripple, we collected a distribution of 80–120 Hz ripple band phases from all 96 microelectrodes, and computed the PPC on that distribution. We assigned the maximum PPC computed over the duration of each iEEG ripple as the microelectrode 80–120 Hz PPC for that iEEG ripple. In each participant, we then computed the correlation between iEEG ripple amplitude and 80–120 Hz PPC in the underlying LFP across all iEEG ripples identified from all retrieval trials. We compared these true correlations to chance using a shuffling procedure. In each of 100 permutations, we circularly shifted the time series of LFP phase by a random amount within each detected iEEG ripple and then computed the correlation between iEEG ripple amplitude and LFP PPC. We calculated the average correlation across permutations in each participant as the chance level. We performed an identical procedure when examining the extent to which the alignment of spiking activity to the 80–120 Hz ripple band signal in the LFP is correlated with the 80–120 Hz ripple band amplitude.

To examine the extent to which spiking activity is locked to individual frequencies in the LFP and iEEG signal, we computed PPC using the instantaneous phases of each spike from each unit. In each participant, we computed the average spike PPC across all units in each trial, and then computed the average across trials to generate a spike PPC value for each participant. In order to compare PPC values across participants, we converted the raw PPC to a z-score in each participant by using the mean and standard deviation of a null distribution of 100 spike PPC values generated by randomly shuffling the trial labels associated with the spike indices.

We then assessed whether the distribution of spike PPC values is significant across participants using a non-parametric cluster-based procedure. For each frequency, we compared the distribution of z-scored spike PPC values to zero using a t-test, thus generating a true $t$-statistic and p-value for each frequency. We then randomly permuted the participant-specific values by randomly reversing the sign of z-scored PPC within each participant and recomputing the average value of the distribution of permuted PPC values across participants. For n participants, this results in an empiric distribution of $2^n$ possible values that are all equally probable under the null hypothesis. We generated an empiric distribution from 1000 permutations for each frequency and calculated $t$-statistics for each of the permuted frequencies.

To correct for multiple comparisons across frequencies, we identified clusters of adjacent frequencies that exhibited a significant difference between the average PPC across participants and zero (where in each frequency cluster, p < .05). For each cluster of significant frequencies identified in the true and permuted cases, we defined a cluster statistic as the sum of the $t$-statistics within that frequency cluster. We retained the maximum cluster statistic during each of the 1000 permutations to create a distribution of maximum cluster statistics. We assigned p-values to each identified cluster of the true data by comparing its cluster statistic to the distribution of maximum cluster statistics from the permuted cases. We determined clusters to be significant and corrected for multiple frequency comparisons if their p-value calculated in this manner was less than .05.

We also compared spike PPC between two sets of conditions - PPC for spikes that occurred during an identified LFP ripple as compared to PPC for spikes that occurred outside an LFP ripple, and PPC for spikes that occurred during correct versus incorrect memory retrieval. We only included units for this analysis that exhibited a minimum of 10 spikes in each condition during an experimental session. In addition, because each condition tends to have a low total number of spikes in each trial, we computed PPC in these analyses by aggregating spiking events across trials rather than initially computing PPC within individual trials. Because we are making a direct comparison between PPC values within individual participants, we used the raw PPC rather than the z-scored value for these tests. In all cases, we computed the average PPC across all units separately for each condition in each participant. We then compared the average PPC between conditions by using a similar permutation procedure that corrects for multiple comparisons described above. In this case, in each permutation we randomly switched the label for each condition in each participant. To ensure that lower spike counts in one condition would not bias our results, we identified which condition had the lower total number of spikes, and randomly subsampled the spikes from the other condition. We performed this subsampling 200 times, computed PPC for each iteration, and assigned the average of the PPC from the 200 iterations of subsampling to the condition with the larger number of spikes. We repeated all of these analyses when examining the extent to which ripples are locked to individual frequencies, and to compare the extent of locking between conditions.

In order to obtain a measure of phase locking that does not depend on number of observations, we look at pairs of phases. Phases that are consistently clustered around a mean phase have a small angular distance to each other. The absolute angular distance is expressed as

where $\phi$ represents the phase of spike to a frequency bin and $\omega$ represents the phase of another spike from the same neuron to the same frequency bin. For each neuron, we can compute this for all frequency bins.

We compute the average pairwise circular distance (APCD), or the absolute angular distance between relative phases, which can be expressed as:

The pairwise phase consistency (PPC) is equivalent to the population statistic of the APCD, which is equivalent to the population statistic of the square of the phase-locking value.

We compute the sample estimate of the PPC by evaluating:

where f(φ,ω) = cos(φ)cos(ω) + sin(φ)sin(ω) and N represents the number of spikes.

To efficiently compute the PPC of spikes to one frequency bin of the local field potential, we express each spike phase as a unit vector and evaluate the dot product for all pairs of unit vectors. We compute the spike-LFP PPC from the resulting symmetric matrix by removing the values along the diagonal and then taking the mean.

To measure the phase consistency of ripple oscillations across MEA electrodes, we compute the absolute angular distance using Equation 3 where ${\theta }_j$ represent the phase of the ripple band signal for one MEA electrode, ${\theta }_k$ represents the phase of ripple band signal for a different MEA electrode for one time point, and N represents the number of MEA electrodes. Each time point within a iEEG ripple was treated as an observation for the MEA electrode. In other words, for a 50 ms long iEEG ripple, we evaluate the dot product for the pairs of ripple phases across all pairs of MEA electrodes. To efficiently compute the PPC of ripple oscillations across MEA electrode pairs, we express each ripple oscillation phase as a unit vector and compute the mean dot product for all pairs of unit vectors in a similar manner as spike-LFP PPC.

To measure the extent to which ripples in the anterior temporal lobe (ATL) are coupled with ripples in the medial temporal lobe (MTL), we identified the time index of peak ripple power for each rippled detected in both regions. We then generated cross-correlograms between MTL and ATL ripples (Vaz et al., 2019). For each electrode in the MTL, we computed a cross-correlogram with each electrode in the ATL. We then pooled these cross-correlograms across trials for each electrode pair in each participant. This generates a cross-correlogram for each pair of electrodes that we can compare between conditions and to a chance distribution (see below). To generate a single cross-correlogram representing the relation between the ATL and the MTL in each participant, we computed the average cross-correlogram across all electrode pairs.

For every pair of electrodes, we calculated a shift predictor for the cross-correlogram that characterizes the cross-correlation that would be expected by chance given the presentation of a stimulus (Vaz et al., 2019; Brody, 1999; Steinmetz et al., 2000; Morris et al., 2004). This chance distribution was generated by cross-correlating the time indices relative to the presentation of the stimulus for each ripple in an ATL electrode during an individual trial with the time indices of each ripple in an MTL electrode in every other trial. For n trials, we create n - 1 cross-correlations, which are then averaged to create a chance cross-correlogram (the shift predictor) for that trial. This procedure was repeated for all trials, and the average across all trials represents the average shift predictor for that trial condition. We aggregated these chance cross-correlograms (shift predictors) across all electrode pairs that involve each region of interest to generate a shift predictor for each region.

The ratio between the true cross-correlogram and the shift predictor reflects the extent to which two signals are synchronized greater than would be expected by chance given the presentation of a stimulus. We calculated a normalized synchronization metric by finding the sum of the true cross-correlation values in a ± 50 ms window and then dividing by the corresponding area of the chance distribution. In this manner, our metric directly quantifies how much more synchronized the true case is relative to chance, which would result in a value of 1. To test the effect of a range of detection parameters on the correlation, we detected ripples using duration thresholds ranging from 10 to 40 ms, increasing in increments of 10 ms, and max amplitude thresholds ranging from 2 to 4 SD, increasing in increments of 1 SD. We used the same detection threshold for LFP and iEEG ripple detection. We used this metric to compare synchronization between detection parameters.

To measure the coincidence of LFP and iEEG ripples, we identified the time index of peak ripple power for each rippled detected in each microelectrode (LFP ripple) and iEEG electrode (iEEG ripple). We then generated cross-correlograms between LFP and iEEG ripples. For each participant, we included four iEEG electrodes nearest to the MEA. To generate a single cross-correlogram representing the relation between the LFP and iEEG ripples in each participant, we computed the average across all electrode pairs. For every pair, we calculated chance cross-correlograms by randomly shifting in time each trial of the ripples detected in the microelectrode. We computed the average across trials for each electrode pair. We calculated a normalized synchronization metric by finding the average true cross-correlation values in a ± 50 ms window and then dividing by the corresponding area of the chance distribution. The ratio between the true and chance cross-correlograms quantifies how much more synchronous the LFP and iEEG ripples are relative to chance, with a value of 1 indicating a measurement equal to chance.

To measure the extent to which units spike together in bursts within detected iEEG ripples, we summed the spiking across all units and computed the auto-correlogram of the population spiking within each detected iEEG ripple. We detected ripples using a duration threshold of 10 ms and an amplitude threshold of 1 SD with a maximum of at least 2 SD in four iEEG electrodes nearest to the MEA. To compare this auto-correlogram within ripples to spiking outside of ripples, we generated random duration matched windows between ripples and computed the chance population spiking auto-correlograms. We calculated a burst metric by finding the average of the true auto-correlogram in a ± 25 ms window centered around zero and then dividing by the corresponding area of the chance correlogram.

To assess whether distributions of population spike rate, LFP ripple power and iEEG ripple power are bimodal, we used Hartigan’s dip test. We postulated that these distributions would be bimodal if there were indeed transient bursts of activity and periods of little activity in between. The dip test computes the maximum difference between the empirical distribution function and the unimodal distribution function that minimizes that maximum difference (Hartigan and Hartigan, 1985). To compute the dip statistic, we generated a probability density function (PDF) of samples aggregated across all four second trials in 200 bins over the range of the data. We computed a true dip statistic for spike rate and for LFP ripple power for each microelectrode and for iEEG ripple power for each iEEG channel. We generated a chance distribution of dip statistics for unimodal distributions to quantify the significance of the true dip statistic. For this procedure, we randomly generated 10,000 uniform PDFs and z-scored the true dip statistic using the mean and standard deviation of the chance distribution. The average z-scored dip statistic across all microelectrodes was used for the spike rate and LFP ripple power for each participant. The average z-scored dip statistic across four iEEG channels in the anterior temporal lobe and iEEG channels in the medial temporal lobe were used to compute the z-scored dip statistic for each participant. This analysis was performed on the six participants with a MEA.

We used an independent and previously validated method for detecting transient episodes of narrowband oscillations to assess whether ripples detected using duration and amplitude thresholds in the 80–120 Hz frequency range capture similar events detected using other approaches. For this procedure, we used the continuous-time wavelet transform (wavelet number 6) to compute the mean power spectrum over the trial, which is then used to generate a background 1/f fit. We generated a 1/f fit to the 70–200 Hz range of the power spectrum for each trial and identify narrowband oscillations that exceed it. The signal is then bandpass filtered within the identified narrowband frequency ranges and a Hilbert transform is used to compute the instantaneous power and phase. The instantaneous frequency is estimated using a frequency sliding estimation method previously described (Cohen, 2014). Periods in which the power is below the 1/f fit is removed. Given we perform this for each trial, we identify a unique narrowband oscillation for each trial for each iEEG electrode. For each participant, we aggregate the oscillations across trials across iEEG electrodes to generate a distribution of center frequencies of narrowband oscillations and a distribution of durations of the periods when the oscillations exceeds 1/f background signal.

Given the variability in number of ripples and other characteristics across participants, we quantified within and across participant variability and computed an estimate of the total true correlations. We assessed whether random variation accounts for the observed correlations by performing a meta-analysis where we used restricted maximum-likelihood estimation to fit a random effects model (Viechtbauer, 2010). For each participant, we computed the true correlation and z-scored it using a distribution of correlation values for shuffled data to generate the r equivalent, a measure of effect size. We computed the sampling variance for each participant from the number of samples (Rosenthal and Rubin, 2003). These measures were used to fit the random effects model.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank J Chapeton, V Sreekumar, and Z Xie for helpful and insightful comments on the manuscript. We are indebted to all patients who have selflessly volunteered their time to participate in this study. This work was supported by the Intramural Research Program of the National Institute of Neurological Disorders and Stroke. This work was also supported by NINDS grant F31 NS113400 (APV). Conflicts of Interest: The authors declare no competing financial interests.

Human subjects: Data were collected at the Clinical Center at the National Institutes of Health (NIH; Bethesda, MD). The Institutional Review Board (IRB) approved the research protocol (11-N-0051), and informed consent was obtained from the participants and their guardians.

This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.