Discrete ripples reflect a spectrum of synchronous spiking activity in human association cortex

Direct brain recordings have provided important insights into how persistent oscillatory activity support human memory retrieval, but the extent to which transient fluctuations in intracranial EEG (iEEG) captures the dynamic coordination of underlying neurons involved in memory processing remains unclear. Here, we simultaneously record iEEG, local field potential (LFP), and single unit activity in the human temporal cortex. We demonstrate that cortical ripples contribute to broadband high frequency activity and exhibit a spectrum of amplitudes and durations related to the amount of underlying neuronal spiking. Ripples in the macro-scale iEEG are related to the number and synchrony of ripples in the micro-scale LFP, which in turn are related to the synchrony of neuronal spiking. Our data suggest that neural activity in the human cortex is organized into dynamic, discrete packets of information.

correspond to the detection of individual ripples ( Figure 1B). Across all trials, ripple rates demonstrate a clear  Figure S2E). Even when ripples have amplitudes or durations below previously used thresholds, spiking activity is present in the microelectrode recording (Supplementary Figure S2F).
While we found a strong relation between spiking activity and ripple amplitude, this observation could be con-143 founded by a correlation between ripple amplitude and duration (Supplementary Figure S2G). Larger amplitude are 144 larger in duration and therefore may have more opportunity to co-occur with spikes by chance. To account for this, 145 we shuffled the time indices of the detected spike times and computed the correlation between LFP ripple amplitude 146 and the spike rate. The true relation between LFP ripple band amplitude and the spike rate is significantly greater 147 than the shuffled distribution (true-shuffled r = .096 ± .020; t(5) = 3.312, p = .011, paired one-tailed t-test; Figure   148 2F).

149
In a similar manner, during every iEEG ripple we determined how many individual units spike as a proportion correlation is significant (Fisher z-transform, r = .15 ± .03; t(5) = 4.679, p = .003, one-tailed t-test). We performed 154 the same shuffling procedure to account for ripple duration, and found that the true correlations across participants 155 are significantly larger than the shuffled data (true-shuffled r = .069 ± .027; t(5) = 2.517, p = .027, paired one-tailed 156 t-test; Figure 2H). Together, these data demonstrate a strong relation between underlying unit spiking activity and 157 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 161 the amplitude of the ripple observed in the iEEG signal is related to both the total number of LFP ripples and the 162 extent to which the LFP ripples observed across the underlying MEA electrodes are aligned ( Figure 3A). 163 We first examined the relation between the amplitude of the iEEG ripple and the number of LFP ripples si-  the PPC is significantly correlated with the maximum amplitude of the observed iEEG ripples (Fisher z-transform; 181 r = 1.06 ± .43; t(5) = 2.34, p = .033, one-tail t-test; Figure 3E; Supplementary Figure S3D). In addition, the 182 correlations across participants are significantly greater than those that would be observed by chance (true-shuffled 183 r = .03 ± .01; t(5) = 2.28, p = .036, paired one-tail t-test; Figure 3E; see Methods). These data together suggest 184 that the iEEG ripple reflects both the aggregate sum and alignment of the underlying LFP ripples.

186
Spiking Activity is Phase-Locked to Ripples

187
Given that synchronization of LFP ripples was associated with increased ripple amplitudes, we also hypothesized 188 that spike timing during ripples would be phase locked (Quyen et al., 2008). In individual participants, we often 189 observed that unit firing preferentially occurs at the trough of the corresponding LFP ripples ( Figure 4A). When 190 the LFP ripples are aligned, spiking activity also appears to preferentially occur at the trough of the overlying iEEG  Figure 4B).

196
When we visualized the spike triggered average of the LFP signal in individual participants, we often observed 197 that spiking activity also appeared locked to low frequency oscillations in the LFP ( Figure 4C). We therefore also 198 examined the distribution of 2-10 Hz low frequency phases present in the LFP signal around each spike and found 199 significant locking to the trough in all participants (p < 10 −4 , Rayleigh test; p = 5.4 × 10 −4 , Rayleigh test across six 200 complex means, one from each participant; Figure 4D). Spikes from all units appear locked around the trough of the  Figure 4D).

205
To examine the relation between spiking activity and individual frequencies within the LFP signal, we computed 206 PPC across all spikes within each MEA electrode for each frequency between 2 Hz and 400 Hz (see Methods) (Vinck 207 et al., 2010). Across participants, spiking activity is significantly locked to specific low and high frequency bands 208 in the LFP (peak at 3.43 Hz and 86.9 Hz, p < .05, permutation test; see Methods Figure 4E). We confirmed that 209 spiking activity is locked to these two frequency bands across participants by also computing the phase-locking value   Brief window around one iEEG ripple showing ripple band iEEG signal (blue), instantaneous phase of a pair of MEA electrodes (purple; out of many pairs, not shown) and instantaneous phase difference of the ch 1 and ch 2 pair (black). Maximum of iEEG ripple indicated with small black square in the shaded window above the ripple band iEEG signal. Polar histogram of all pairwise phase differences during a detected iEEG ripple is centered around 0 (blue). Polar histogram all pairwise phase differences outside of a iEEG ripple is more uniform (black). (E) Fisher z-transformed Pearson correlations between maximum pairwise phase consistency across all MEA electrode pairs and maximum amplitude of iEEG ripples. Group level results, reported as mean ± SEM, persists when duration of the iEEG ripple is accounted for by shuffling ( * p < 0.05). Each data point represents a participant.  (Supplementary Figure S4A). Spikes are significantly more locked to low and high frequencies when they arise during 211 ripples as compared to between ripples (p < .05, permutation test; Figure 4F, Supplementary Figure S4B). We were 212 also interested in whether the locking between spiking activity and LFP oscillations is related to successful memory 213 retrieval. We examined only spikes that occur during identified LFP ripples since ripples have been previously linked 214 with successful memory retrieval (Norman et al., 2019;Vaz et al., 2019). Across participants, we found that locking 215 between spiking activity and low frequency oscillations appears greater during correct compared to incorrect memory 216 retrieval, but the difference is not consistent across participants (p > .05, permutation test; Supplementary Figure   217 S4C,D). 218 We next examined the relation between the extent to which spiking activity locks to the 80-120 Hz frequency 219 within each ripple and the amplitude of the ripple. Across participants, mean spike-LFP PPC within the 80-120 Hz 220 ripple band is significantly correlated with 80-120 Hz ripple amplitude across all MEA electrodes (Fisher z-transform, 221 r = 0.084 ± 0.026, t(5) = 3.27, p = .011; Figure 4G). To account for any possible effects of ripple duration on the 222 calculation of PPC, we compared this true distribution to a chance distribution and found that across participants 223 LFP ripple amplitude exhibits a significantly stronger correlation with spike locking to the 80-120 Hz band in the true 224 data as compared to the chance distribution (t(5) = 2.64, p = .023; see Methods). Together with our data examining 225 the relation between spiking activity and ripple amplitude, these data suggest that the amplitude of ripples in the 226 LFP signal may reflect both the sum and alignment of underlying spiking.

227
Finally, given the observed relation between spiking activity and ripples, we then examined whether ripples 228 themselves also exhibit a phase preference. As with the individual spikes, we considered each LFP ripple as an event

236
In the same manner as with spiking activity, we computed locking between LFP ripples and oscillations at each 237 logarithmically spaced frequency between 2 Hz and 400 Hz. In each participant, LFP ripple activity is significantly 238 locked to low frequencies in the LFP signal (p < .05, permutation test; Supplementary Figure S4F). We did not find 239 a consistent difference in the locking of LFP ripples to LFP low frequency oscillations when comparing correct and 240 incorrect memory retrieval (Supplementary Figure S4G). However, when we examined the extent to which microscale 241 LFP ripples are locked to the larger iEEG signal, we found that LFP ripples that are present during iEEG ripples are 242 significantly more locked to the lower frequency oscillations during correct compared to incorrect memory retrieval

Discussion
Despite significant advances over the past several decades, how to accurately interpret the various fluctuations and 247 dynamics observed through direct recordings of the human brain has remained challenging. Our results demonstrate 248 that many of the changes in broadband high frequency activity captured through iEEG reflect discrete 80-120 Hz 249 ripple oscillations that in turn reflect underlying bursts of spiking activity. These discrete and short bouts of neuronal 250 activity may constitute packets of information that underlie the neural substrates of human cognition.

251
One of the challenges in examining the role of ripple oscillations in cognition, especially in the human brain, has 252 been determining whether any particular event does or does not qualify as a ripple. Many of the criteria used for  Our data demonstrate that cortical ripples captured through human brain recordings exist on a continuum of 260 amplitudes and durations. By recording neural activity across spatial scales, we find that even ripples with smaller in our data therefore suggests that strictly adhering to predefined criteria for what constitutes a ripple may run the 269 risk of overlooking functionally meaningful events across these spatial scales.

270
By examining both iEEG and LFP recordings in the human brain for the presence of discrete ripples, our data 271 also support the hypothesis that many of the dynamics observed in broadband high frequency activity captured from 272 the human brain are driven by well-defined and brief bouts of neural activity. A common approach for investigating 273 the neural correlates of human cognition has been to average neural activity over multiple trials and over broad  The relation between band limited 80-120 Hz ripples and broadband high frequency activity that we observe in our 278 data suggests that many of the interpretations regarding the neural substrates of human memory may be better 279 served by considering these discrete events. It is important to recognize, however, that this relation is not absolute 280 and appears less robust outside of the MTL and ATL. Even within these brain regions, this relation is clearer only 281 during correct compared to incorrect memory retrieval. Hence, while discrete 80-120 Hz ripples may underlie many of the phenomena observed through broadband high frequency activity, there are likely other neural mechanisms 283 that contribute to the dynamics observed in the iEEG signal.

284
The possibility that information is neurally encoded through packets of activity has been relatively under-explored 285 in human brain recordings. Recent evidence captured through animal recordings related to both memory and 286 perception, however, supports this possibility (Luczak et al., 2009(Luczak et al., , 2015Lundqvist et al., 2016). These advances are 287 partly due to the more sophisticated tools that are available for in vivo recordings of large populations of spiking 288 neurons in animals. By recording spiking activity from a population of neurons in the human temporal cortex through 289 microelectrode arrays, we find direct evidence that 80-120 Hz ripples that we observe in our data is accompanied by 290 a burst of neuronal spiking. Hence, our data demonstrate that neural activity in the human temporal cortex may 291 be temporally organized into discrete packets of information. Our data focus on these packets as participants form 292 and retrieve memories since the information contained within these packets has been recently linked with memory 293 retrieval (Pfeiffer, 2020; Vaz et al., 2020). However, our data cannot address whether the relation between ripples 294 and underlying bursts of synchronized spiking is unique to just the temporal lobe or just to memory. Ripples have 295 been most studied in the MTL in both animals and humans, but appear to jointly occur in brain regions that either . It is possible that this relation between ripples and spiking activity is unique to brain regions 298 that communicate directly with the MTL or that are directly involved in memory. Our results, however, raise the 299 possibility that such packets may be a general feature of neural coding in the human brain.

300
Given previous evidence demonstrating that spiking activity within ripples appears locked to the trough of each 301 cycle, it is not surprising that we observe similar locking in our data (Nitzan et al., 2020;Quyen et al., 2008). We find 302 more consistent locking of spiking activity to higher frequencies in the micro-compared to macro-scale. This may be 303 because synchronous spiking can occur within local neuronal ensembles while varying across ensembles. However, we 304 also find that spikes, and consequently ripples themselves, are also locked to low frequency oscillations. Such locking 305 of both spiking activity and ripples to low frequencies can account for several phenomena that have been previously 306 described in human brain recordings. For example, phase amplitude coupling between low frequency oscillations and  More broadly, the locking to lower frequencies also suggests a possible mechanism by which packets of information 312 may be conveyed from one brain region to another. Oscillations observed in the iEEG have been hypothesized to to another. Recent evidence has also suggested that higher frequency oscillations that are synchronous across brain regions may also facilitate communication, although the evidence for this still remains unclear (Bosman et al., 2012; coherence, as two brain regions that each exhibit bursts of spiking activity, either conveyed from one to the other 321 directly or driven by a third region, can each generate high frequency ripple oscillations. If the underlying neuronal 322 interactions in each brain region are similar, the ripples may appear coherent and at the same high frequency.

323
Conversely, if the underlying architecture of each region is different, then any resulting higher frequency oscillations 324 may differ in morphology and frequency, and therefore appear desynchronized.

325
Together, our data offer insights into the dynamic fluctuations observed in direct recordings from the human brain 326 and suggest that neural activity may be organized into discrete 80-120 Hz ripple events that reflect underlying bursts 327 of neuronal spiking. Our data argue against using fixed criteria to identify these ripples, and instead demonstrate 328 that these ripples exist on a continuum of activity. As each of these 80-120 Hz ripples reflects bursts of neuronal 329 spiking with varying degrees of activity, our data more broadly suggest that these discrete packets of activity may    During the test period, one word was randomly chosen from each of the presented pairs and presented in random 370 order, and the participant was asked to recall the other word from the pair by vocalizing a response. Each cue 371 word was preceded by an orientation stimulus (a row of question marks) that appeared on the screen for 4000 ms 372 followed by a blank ISI of 1000 ms. Participants could vocalize their response any time during the recall period after 373 cue presentation. We manually designated each recorded response as correct, intrusion, or pass. A response was 374 designated as pass when no vocalization was made or when the participant vocalized the word 'pass'. During pass 375 trials where no vocalization was present, we assigned a response time by randomly drawing from the distribution of 376 correct response time during that experimental session. We defined all intrusion and pass trials as incorrect trials.   Figure S1). Subdural contacts were arranged in both grid and strip configurations with 383 an inter-contact spacing of 10 mm. We captured iEEG signals sampled at 1000 Hz. For clinical visual inspection of 384 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 386 was set by the recording amplifier (Nihon Kohden, Irvine CA) as the average of two intracranial electrode channels. 387 We used the Chronux toolbox to apply a local detrending procedure to remove slow fluctuations ( 2Hz) from the 388 time series of each electrode and a regression-based approach to remove line noise at 60Hz and 120Hz (Mitra & Bokil,389 2009). We implemented additional thresholds to remove movement artifacts and pathological activity related to the 390 patient's epilepsy (SI Appendix).    there were 669.44 ± 74.30 vertices per ROI and each vertex mapped to 8.08 ± 0.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 421 ROIs as its nearest pial vertices. Due to the overlap between ROIs, each electrode is assigned to multiple ROIs and 422 each ROI may contain more than one electrode. For analyses within ROIs across participants, we only included ROIs 423 that contained electrodes from at least five participants.    classified as belonging to that unit, and the noise cluster consisted of all waveforms that crossed the threshold that 474 were not classified as belonging to any unit. The isolation score is normalized to be between 0 and 1, and serves 475 as a measure to compare the isolation quality of all units across all experimental sessions and participants. Across 476 participants, the mean isolation score for all units was 0.93 ± 0.1.

477
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 479 factor (set as 5). To obtain N oise, we subtracted the mean waveform from each individual waveform for each 480 identified unit, concatenated these waveform residuals, and then computed the standard deviation of this long 481 vector. Therefore, the noise term quantifies the within-unit variability in waveform shape. Across participants, the 482 mean SNR for all units was 1.71 ± 0.12. 483 We estimated the instantaneous spike rate for each unit by convolving the spike rasters with a Gaussian kernel 484 20 (σ = 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.

486
Ripple Detection 488 We detected ripples in both the iEEG and LFP signals as previously reported (Vaz et al., 2019). We first bandpass 489 filtered the voltage time series in the ripple band (80-120 Hz) using a second order Butterworth filter, and then 490 applied a Hilbert transform to extract the instantaneous amplitude and phase within that band. We selected events 491 where the Hilbert envelope exceeded 2 standard deviations above the mean amplitude of the filtered traces. We 492 only retained events that were at least 25 ms in duration and had a maximum amplitude greater than 3 standard 493 deviations as ripples for analysis. We joined adjacent ripples that were separated by less than 15 ms. We identified 494 every ripple that satisfied these criteria in every electrode contact, and assigned each such identified ripple a start time 495 index and an end time index. The difference between them defined the duration of each ripple. As described above, 496 we excluded all IEDs and high frequency oscillations associated with IEDs (ripple on spike waveforms, pathologic 497 ripples) from our analyses. The remaining ripples that we retained for our analyses therefore occurred without an 498 associated IED and are more likely to be physiologic.

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

506
To account for the possibility that ripples with higher amplitudes and therefore longer durations may be associ-507 ated with more spiking activity by chance, we compared the true correlation between ripple amplitude and spiking 508 activity to the correlations we would observe by chance. In each of 1000 permutations, we performed a random

557
To correct for multiple comparisons across frequencies, we identified clusters of adjacent frequencies that exhibited 558 a significant difference between the average PPC across participants and zero (where in each frequency cluster, 559 p < 0.05). For each cluster of significant frequencies identified in the true and permuted cases, we defined a cluster 560 statistic as the sum of the t-statistics within that frequency cluster. We retained the maximum cluster statistic 561 during each of the 1000 permutations to create a distribution of maximum cluster statistics. We assigned p-values 562 to each identified cluster of the true data by comparing its cluster statistic to the distribution of maximum cluster 563 statistics from the permuted cases. We determined clusters to be significant and corrected for multiple frequency 564 comparisons if their p-value calculated in this manner was less than 0.05. 565 We also compared spike PPC between two sets of conditions -PPC for spikes that occurred during an identified 566 LFP ripple as compared to PPC for spikes that occurred outside an LFP ripple, and PPC for spikes that occurred In order to obtain a measure of phase locking that does not depend on number of observations, we look at pairs of 584 phases. Phases that are consistently clustered around a mean phase have a small angular distance to each other.

585
The absolute angular distance is expressed as where ϕ represents the phase of spike to a frequency bin and ω 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. 588 We compute the average pairwise circular distance (APCD), or the absolute angular distance between relative 589 phases, which can be expressed as: The pairwise phase consistency (PPC) is equivalent to the population statistic of the APCD, which is equivalent 591 to the population statistic of the square of the phase-locking value. 592 We compute the sample estimate of the PPC by evaluating: where f (ϕ, ω) = cos(ϕ)cos(ω) + sin(ϕ)sin(ω) and N represents the number of spikes.

594
To efficiently compute the PPC of spikes to one frequency bin of the local field potential, we express each spike 595 phase as a unit vector and evaluate the dot product for all pairs of unit vectors. We compute the spike-LFP PPC 596 from the resulting symmetric matrix by removing the values along the diagonal and then taking the mean.      Group level statistics are shown as mean ± SEM across six participants. Given the concern that spikes in the signal may generate spectral artifacts in the ripple band, we performed two control analyses to confirm that the significant correlation between LFP ripple amplitude and spike rate was not due to bandpass filtering over spikes. In the first, we removed spikes from the LFP data and in the second we restricted our analysis only to MEA electrodes that exhibited no spiking. We compared the true correlation between spike rate and LFP ripple amplitude with the correlations observed after spike removal and interpolation (top) and channels without spiking (bottom; paired t-test, t(5) = 1.29, p = .254; orig, t(5) = 5.32, p = .003; spike removal, t(5) = 4.78, p = .005). (F) Relation between number of spikes and LFP ripple duration (top) and amplitude (bottom) for all ripples in all participants. Each color represents a participant (n = 6). (G) Pearson correlation between amplitude and duration of LFP ripples (left) and iEEG ripples (right) compared to Pearson correlations after random circular shifts of ripple indices by trial. The true relation between ripple band amplitude and duration is significantly greater than the shuffled distribution (true r = .372 ± .033; true-shuffled t(5) = 9.07, p = 1.4 × 10 −4 , paired one-tailed t-test).   PLV confirms that spiking activity is locked to low and high frequency activity across participants (peak at 2.6 Hz and 86.9 Hz, p < .05, permutation test). (B) Spike-LFP pairwise phase consistency (PPC) for spikes that occur with LFP ripples and for spikes that do not occur with LFP ripples for each participant, shown as mean ± SEM across MEA electrodes. (C) Spike-LFP PPC for correct and incorrect memory retrieval for spikes that co-occur with LFP ripples, shown as mean ± SEM across MEA electrodes for each participant (top). Difference in spike-LFP PPC between correct and incorrect memory retrieval shown as mean ± SEM across participants (bottom). (D) Difference in spike-LFP PPC between correct and incorrect memory retrieval for spikes that co-occur with LFP ripples in all MEA electrodes from all six participants across all frequencies. Grand average ± SEM of the difference in spike-LFP PPC between correct and incorrect memory retrieval across all MEA electrodes (bottom). (E) LFP ripple triggered average (RTA) for LFP ripples that co-occur with iEEG ripples, in purple, with the bandpass filtered signal, in black, for one representative participant. Also shown is the bandpass filtered RTA during correct memory retrieval (green) and incorrect trials (orange). Distribution of phases of LFP low frequency (lower left histogram) and iEEG low frequency (lower right histogram) signals across LFP ripple times for all MEA electrodes. Complex mean of the distribution of phases for each participant is depicted in a polar plot with circles filled with a star if the distribution shows significant phase-locking (Rayleigh test, p < 0.001). Black line shows the average of six distributions across participants. (F) LFP ripple-LFP PPC for each frequency shown as mean ± SEM across participants. LFP ripples are locked to low frequency activity (peak at 2.39 Hz, p < .05, permutation test). (G) LFP ripple-LFP PPC for correct and incorrect memory retrieval, shown as mean ± SEM across MEA electrodes for each participant. (H) LFP ripple-iEEG PPC for correct and incorrect memory retrieval for LFP ripples that co-occur with iEEG ripples, shown as mean ± SEM across MEA electrodes for each participant (top). Difference in LFP ripple-iEEG PPC between correct and incorrect memory retrieval shown as mean ± SEM across participants (bottom). LFP ripples that occur within iEEG ripples are preferentially locked to low frequency activity across participants (peak at 5.50 Hz and 9.97 Hz, p < .05, permutation test). (I) Difference in LFP ripple-iEEG PPC between correct and incorrect memory retrieval for LFP ripples that co-occur with iEEG ripples across all MEA electrodes from all six participants. Grand average ± SEM of the difference in LFP ripples-LFP PPC between correct and incorrect memory retrieval across all MEA electrodes. 32