Independent dynamics of low, intermediate, and high frequency spectral intracranial EEG activities during human memory formation

A wide spectrum of brain rhythms are engaged throughout the human cortex in cognitive functions. How the rhythms of various frequency ranges are coordinated across the space of the human cortex and time of memory processing is inconclusive. They can either be coordinated together across the frequency spectrum at the same cortical site and time or induced independently in particular bands. We used a large dataset of human intracranial electroencephalography (iEEG) to parse the spatiotemporal dynamics of spectral activities induced during formation of verbal memories. Encoding of words for subsequent free recall activated low frequency theta, intermediate frequency alpha and beta, and high frequency gamma power in a mosaic pattern of discrete cortical sites. A majority of the cortical sites recorded activity in only one of these frequencies, except for the visual cortex where spectral power was induced across multiple bands. Each frequency band showed characteristic dynamics of the induced power specific to cortical area and hemisphere. The power of the low, intermediate, and high frequency activities propagated in independent sequences across the visual, temporal and prefrontal cortical areas throughout subsequent phases of memory encoding. Our results provide a holistic, simplified model of the spectral activities engaged in the formation of human memory, suggesting an anatomically and temporally distributed mosaic of coordinated brain rhythms.


Introduction
Brain rhythms are thought to support memory and cognitive functions by coordinating activities of connected neurons into synchronous interactions ( Singer, 1999 ;Buzsaki, 2006 ;Fries, 2015 ;Klimesch, 1996 ;Fell and Axmacher, 2011 ). The spectrum of rhythmic activities involved in these interactions includes the theta (3-9 Hz) and the gamma (30-120 Hz) frequency bands with specific roles proposed for the low and high frequency activities in memory processing ( Düzel et al., 2010 ;Sauseng et al., 2010 ;Fries et al., 2007 ;Tallon-Baudry and Bertrand, 1999 ;Nyhus and Curran, 2010;Lisman and Jensen, 2013;Herweg et al., 2020) . The intermediate frequency alpha (9)(10)(11)(12) and mation at a defined timescale, engaging only a subset of the spectral activities in a given frequency range. For example, only theta rhythms would be indued at a given site. Both scenarios agree that the brain activities at defined frequency ranges can be conceptualized as "spectral fingerprints " that subserve specific perceptual or cognitive processes at a corresponding temporal scale of neural interactions ( Siegel et al., 2012 ;Hanslmayr and Staudigl, 2014) . Resolving the dynamics of these spectral fingerprints accurately in the anatomical space and time is critical to test these hypotheses and provide a holistic model of the neural activities engaged during cognitive tasks.
Intracranial electroencephalographic (iEEG) recordings present a unique opportunity to probe with high resolution the large scale of neural activities engaged during cognitive and other brain functions in discrete areas of the human cortex ( Lhatoo et al., 2019 ;Johnson et al., 2020 ;Jacobs and Kahana 2010 ;Engel et al., 2005 ). The iEEG signals are sampled from electrode contacts implanted either directly on the surface of the neocortex or into its deeper layers and subcortical structures. Thus, spectral activities generated by local neural populations and recorded from multiple separate contacts can be tracked in discrete cortical areas during the time of memory processing. Previous iEEG studies commonly employed the spectral power induced in high gamma frequency ranges (60-120 Hz) to track cortical processing in memory and cognitive tasks ( Jerbi et al., 2009 ;Lachaux et al., 2012 ;Crone et al., 2006 ;Lundqvist et al., 2018 ). Theta ( Sheehan et al., 2018 ;Lin et al., 2017 ) and alpha ( Staresina et al., 2016 ) were also investigated previously. Oscillatory activities in the theta and alpha frequencies were shown to propagate as local, independent, travelling waves in discrete cortical areas during memory processing ( Zhang et al., 2018 ) , but the relative spatiotemporal dynamics of these low frequency oscillations and the faster beta and gamma activities during memory processing remain largely unexplored.
Here, we took advantage of a large iEEG dataset of 164 participants encoding lists of words for subsequent free verbal recall. Previous studies with this task revealed a distributed network of brain regions associated with sensory processing and with higher-order declarative memory functions ( Burke et al., 2013 ;Burke et al., 2014 a;Kucewicz et al., 2019 ). The induced high gamma power was observed first in the visual areas of the occipito-temporal cortex and then in more anterior areas of the prefrontal cortex, inspiring a two-stage model of memory encoding with an early sensory and a late association stage (Burke et al., 2014a) . A similar sequence of posterior-to-anterior order of induced power was also observed in the theta frequency band ( Burke et al., 2013 ) and confirmed in activities beyond the gamma frequency range ( Kucewicz et al., 2014 ). Our recent study found the sequence of the induced high gamma power during memory encoding to be continuous with increasing latencies along a hierarchy of gradually higher-order association areas ( Kucewicz et al., 2019 ), culminating in the anterior prefrontal cortex late into word encoding. The previous iEEG studies of verbal memory processing were either limited to smaller datasets of participant recordings with sparse electrode coverage of the neocortex or focused only on selected frequency bands. Hence, spectral power across various bands would either be averaged spatially across multiple electrode sites of larger cortical regions, thereby losing the resolution and granularity of individual sites, or would not be studied across the theta, alpha and beta, and gamma brain activities.
Given the wide scope of the frequency spectrum and brain areas analyzed collectively in the previous studies, we hypothesized that specific spectral activities are induced in a granular pattern of discrete cortical locations. Hence, aggregating activity from several locations would give rise to a broadband tilt in power across the entire spectrum ( Kilner et al., 2005 ), which could actually be composed of multiple, specific, "spectral fingerprints" ( Fellner et al., 2019 ). Exact localization of different fingerprints would either be common to the same cortical sites of electrode recording or would be distributed in a mosaic of distinct anatomical sites. Likewise, in terms of the timing, the low, intermediate, and high frequency spectral activities would be induced at different times of memory encoding or all at the same time. In addition to the theta and gamma activities, alpha and beta rhythms were expected to be induced at distinct cortical sites and times of memory encoding. In general, we tested for independent spatiotemporal dynamics of the spectral activities induced during formation of human verbal memory traces. Our goal was to provide a holistic model of brain rhythms engaged across a broad range of frequencies, cortical anatomy, and phases of memory encoding.

Subjects
The dataset of iEEG recordings from a total of 164 participants undergoing surgical evaluation of drug-resistant epilepsy was taken from a large multicenter collaborative study (all de-identified data are available at http://memory.psych.upenn.edu/Electrophysiological_Data ). The recordings were collected from the following centers: Mayo Clinic, Thomas Jefferson University Hospital, Hospital of the University of Pennsylvania, Dartmouth-Hitchcock Medical Center, Emory University Hospital, University of Texas Southwestern Medical Center, and Columbia University Hospital. One common research protocol was approved by the respective Institutional Review Board at each center, and informed consent was obtained from each participant. For the 139 patients included in the final analysis, information regarding number of word trials complete and seizure onset zone are summarized in Suppl. Table 1. Seizure onset was localized to the right hemisphere in 46 patients, the left hemisphere in 55, bilaterally in 26, and undetermined in 12. Demographic information such as age and sex were only available for 85 patients (mean age 35.6 years with standard deviation 11.6 years; 43 Females to 42 Males, Suppl. Table 1) due to limited access to the sensitive health-care information from multiple clinical centers that collected the aggregate dataset. Electrophysiological recordings were collected from standard clinical subdural and depth electrodes (AdTech Inc., PMT Inc.) implanted on the cortical surface and into the brain parenchyma, respectively. Subdural electrode contacts were arranged either in a grid or a strip configuration with 10 mm separation, and depth electrode contacts were separated by 5 to 10 mm. The placement of electrodes was determined by a clinical team with the goal of localizing seizure foci for possible epilepsy resective surgery or implantation of a device for therapeutic electrical brain stimulation.

Anatomic localization and brain surface mapping
Cortical surface parcellations were generated for each participant from pre-implant MRI scans (volumetric T1-weighted sequences) using Freesurfer software (RRID:SCR_001847). The hippocampus and surrounding cortical regions were delineated separately based on an additional 2 mm thick coronal T2-weighted scan using the Automatic Segmentation of Hippocampal Subfields (ASHS) multi-atlas segmentation method. Electrode contact coordinates derived from co-registered postimplant CT scans were then mapped to the pre-implant MRI scans to determine their anatomic locations. For subdural strips and grids, the electrode contacts were additionally projected to the cortical surface using an energy minimization algorithm to account for postoperative brain shift. Contact locations were reviewed and confirmed on surfaces and cross-sectional images by a neuroradiologist. The T1-weighted MRI scans were also registered to the MNI152 standard brain to enable comparison of recording sites in a common space across subjects. Anatomic locations of the recording sites, including Brodmann areas, were derived by converting MNI coordinates to Talairach space and querying the Talairach daemon (www.talairach.org).

Electrophysiological recordings
The iEEG signals were recorded using one of the following clinical electrophysiological acquisition systems (dependent on the institution for data collection): Nihon Kohden EEG-1200, Natus XLTek EMU 128, or Grass Aura-LTM64. Depending on the acquisition system and the preference of the clinical team, the signals were sampled at either 500, 1000, or 1600 Hz and were referenced to a common contact placed either intracranially, on the scalp, or on the mastoid process. For analysis, all recordings using higher sampling rates were filtered by an antialiasing filter and down-sampled to 500 Hz. A common bipolar montage was calculated post hoc for each subject by subtracting measured voltage time series on all pairs of spatially adjacent contacts. This resulted in N-1 bipolar signals in case of the penetrating and the strip electrodes, and N = (i − 1) * j + (j − 1) * i, where i and j are the numbers of contacts in the vertical and horizontal dimensions of the grid. For the data analysis of this study, one "electrode " refers to the bipolar signal from one bipolar pair of contacts. It is important to note that contacts may be used for multiple electrodes, so not all signals are truly independent.

Free recall task
A classic paradigm for probing formation of verbal episodic memory was employed ( Kahana, 2014 ), in which subjects were shown lists of words for subsequent free recall. Participants were instructed to study lists of words presented on a laptop computer screen for a delayed test of vocal recall. Lists were composed of 12 words chosen at random and without replacement from a pool of high-frequency nouns in the subject's native language (either English or Spanish; http://memory.psych.upenn.edu/WordPools). Each word appeared on the screen for 1600 ms, followed by a random jitter of 750 to 1000 ms blank interval between stimuli. At the end of each word list, a distractor task was performed by the subject. This task lasted for 20 s and was composed of a series of simple, arithmetic problems of the format A + B + C, where A, B , and C were random, single-digit integers between 1 and 9. After the distractor task, participants were given 30 s in which to recall as many of the words from the list as possible in any order ( Fig. 1 ). Vocal responses were recorded digitally by the laptop and were later scored manually for analysis. Only subjects who recalled at least 15% of words and completed at least 12 lists of the task were included in further analysis ( Long et al., 2014 ). This left 139 of the 164 original subjects for a grand total of 14,219 electrodes used in this study. Electrophysiological recordings were synchronized to stimulus appearance on the screen through an electric pulse generator operated by the task laptop, which sent pulses to a designated event channel in the clinical acquisition system. The events were timestamped after the recording session using custom-written MATLAB codes and were used to extract specific epochs of interest around word presentation. Each recorded epoch was 3000 ms long and included 1600 ms of word presentation on the screen with 700 ms of blank screen of interstimulus interval before and after each word presentation.
Trial-averaged spectral power of iEEG signals from epochs of word presentation (one word per trial) was integrated across the epoch and used as a feature to classify active electrodes that record from brain regions engaged in verbal memory encoding. Active electrodes were classified independently in six frequency bands of the iEEG spectrum using normalized estimates of the power change (z-score transform). The classification procedure used an unsupervised method based on the Gaussian Mixture Model . Notice induced activity across all six frequency bands in the example spectrogram plotted from trialaveraged activity from one electrode localized in the occipital cortex.

Data analysis
We analyzed the iEEG recording epochs during word presentation for memory encoding. Each presentation epoch was FIR filtered (2000order Barlett-Hanning with zero-phase distortion filtering, bandpass with frequency limits specified for each frequency band) before being spectrally decomposed, normalized, and binned independently into distinct frequency bands between 2 and 120 Hz: low theta (2-4 Hz), high theta (5-9 Hz), alpha (10-15 Hz), beta (16-25 Hz), low gamma (25-55 Hz), and high gamma (65-115 Hz). The boundaries for these frequency bands were determined based on non-overlapping and continuous ranges growing in width along the frequency spectrum. The upper limit of 115 Hz was chosen to ensure at least 4 samples per cycle given the sampling frequency of 500 Hz and to avoid the harmonic of the 60 Hz line noise. To convert filtered signals to a spectrogram with a 2 Hz frequency bin resolution from 2 Hz to 115 Hz, we used 500 ms sliding windows with 5 ms slide lengths and the multi-taper method with 1 taper from the chronux toolbox in MATLAB (http://chronux.org/; Mitra, 2007 ). Spectrograms of the word presentation epoch were averaged together after log-normalization and a z-score transformation of each time-frequency point. Z-scoring was performed according to the following formula: where x is the signal, t is the time bin, f is the frequency bin, f is the mean, and f is the standard deviation for the given frequency. Each trial was normalized separately. Log normalization was implemented in order to correct for the 1/f-power law effect of lower frequencies on estimating power in the higher frequency bands, and the z-score transformation was performed to provide a normalized scale of power change for comparison of signals from different electrodes, sessions and participants ( Kucewicz et al., 2019Saboo et al., 2019 ). This z-score normalization produces positive and negative power changes relative to the mean power within any one word epoch instead of using a "baseline " period. Hence, even small positive or negative deviations of absolute power would effectively be augmented relative to the signal mean, as described before ( Kucewicz et al., 2019 ;Alotaiby et al., 2015 ). Edge artifacts in the power estimate were eliminated by clipping two time bins from either end of the final spectrogram for all frequencies. For each electrode and frequency band, we determined the average spectral power at each time bin across all word presentation epochs. The overall induced power in a given frequency band during the word presentation period was quantified as the area under the absolute value of the corresponding time-series curve, as shown in Fig. 1 . Electrodes were then classified into "active " and "inactive " in each specific band using our unsupervised clustering method based on Gaussian Mixture Modeling (GMM) to identify those that showed task-induced power changes . This method essentially performs a binary machine-learning classification of active and inactive channels based on the induced power value of each channel. The classification was done separately for each frequency band of interest. The entire procedure is summarized in Fig. 1 . We chose a GMM-based method to ensure that different numbers of electrodes would be allowed for the clusters of active electrodes and the inactive electrodes. The number of clusters was set to two because the active-inactive categorization of electrodes is binary: change in the spectral power can either be induced (increased or decreased) or not during the word presentation period. The use of a GMM also enabled classification without a requirement for ground truth data. Electrodes were pooled across subjects and then clustered to identify active electrodes.

Induced power mapping
Most of the active electrodes were localized in the cortical regions associated with processing visual and semantic information, as well as with declarative memory and executive functions. Out of 39 Brodmann Areas, the majority of electrode locations grouped into 9 brain regions each comprising Brodmann Areas: visual (Brodmann areas 18 and 19), inferior temporal (Brodmann areas 20 and 37), precuneus (Brodmann areas 30 and 31), lateral parietal (Brodmann areas 7 and 39), mesial temporal (Brodmann area 28 and Hippocampus), lateral temporal (Brodmann areas 21 and 22), Broca's areas (Brodmann areas 44 and 45), lateral prefrontal (Brodmann areas 9 and 46), and frontal pole (Brodmann areas 10 and 11). Unless stated otherwise, all further analysis was focused on the recalled word epochs, i.e. epochs with successfully encoded words that were subsequently freely recalled, in these nine cortical regions.
Induced power changes in each frequency band were plotted from all active electrodes onto an average brain surface for each time point of the word epoch using a custom MATLAB code for an electrode-locationdependent weighing technique. The brain surface was created as a mesh, where each element of the mesh was assigned a color dependent on the induced power of the surrounding electrodes weighed by the relative distance to the element (greater weight to more proximal electrodes). For each mesh element of the brain surface plotted, only the effects of electrodes within a threshold of 10 mm from the element center point were considered. Hence, the induced power was averaged from all electrodes within the threshold radius that were linearly weighted by their proximity according to the following formula: Weighted Value = I P 1 * (threshold− r 1 ) + I P 2 * (threshold− r 2 ) + … + I P n * (threshold− r n ) ( threshold− r 1 ) where IP is the induced power change value for an electrode, and r is the distance between the electrode and the chosen element of the brain surface mesh.

Statistics
ANOVA tests with multiple comparison corrections for the effects of frequency band, brain region (Brodmann Area), and hemisphere were used to analyze the percentage of electrodes in an area labeled as "active " ( Fig. 2 ). Temporal differences between hemispheres were quantified using a two-tailed, unpaired t -test statistic with Bonferroni correction ( Fig. 3 A). We then broke each word epoch into pre-encoding, early encoding, late-encoding, and post-encoding phases ( Burke et al., 2014 ;Kucewicz et al., 2019 ) and used repeated measures ANOVA (encoding phase as the within subjects factor) with post-hoc Tukey-Kramer test to analyze the effects of frequency band, brain region, and hemisphere on induced power during each phase of the verbal memory task ( Fig. 3 B). Peak latency for each frequency band and area was determined as the time bin with the maximum positive peak of the induced power. Latencies were compared using repeated measures ANOVA (encoding phase as the within subjects factor) and post-hoc Tukey-Kramer multiple comparison tests ( Fig. 4 ). Pearson's correlation test was used to assess the effect of position along each axis of the brain on the time of peak power latency ( Fig. 5 B). The alpha value for all tests was set at 0.05.

Spectral activities induced during memory encoding are heterogeneously distributed across the human cortex
We first identified a subset of electrodes from cortical areas activated during memory encoding. For this purpose, spectral power-inband in six non-overlapping frequency ranges (low theta: 2-4 Hz, high theta: 5-9 Hz, alpha: 10-15 Hz, beta: 16-25 Hz, low gamma: 25-55 Hz, high gamma: 65-115 Hz) was used as features for automated classification of active electrodes  independently in each band ( Fig. 1 ). Out of the total of 14,219 electrodes implanted across all participants, 4738 (33.3%) showed the induced spectral power in at least one of the frequency bands. Hence, a typical electrode site showed induced power in a defined range of either low, intermediate, or high frequency bands ( Fig. 2 A). In agreement with our previous results ( Kucewicz et al., 2019 ;Alotaiby et al., 2015 ), these active electrodes were widely distributed across 39 Brodmann areas that were defined in all cortical lobes ( Fig. 2 B). The relative proportion of electrodes in each Brodmann area that were active varied from 4.3% to 75%, depending on location within the limbic, frontal, prefrontal, parietal, temporal, or occipital lobes (ANOVA F = 99.21, p < 0.001, df = 5), with significantly more active electrodes found in the sensory visual areas of the occipital and the parietal cortex (post-hoc Tukey Kramer, p < 0.05). Within each area, similar proportions were observed in each frequency band, including the alpha and the beta band, with no significant differences across the spectrum (N-way ANOVA F = 0.17, p = 0.9747, df = 5). In summary, regardless of frequency band analyzed, at least 4% of electrode sites within any of the cortical areas analyzed revealed significantly induced power during the task.
To investigate whether these various low and high frequency spectral activities were induced on the same electrodes, we summarized the overlap across the six frequency bands for nine cortical regions of interest (ROI; each composed of two Brodmann areas with > 18 active electrodes from > 8 participants; see Table 1 ). An electrode can be active in one or more bands. We found that the highest proportion of electrodes were active in only one (43.12%) or in two (22.94%) bands, and proportionally fewer in more than two bands with only a minority (5.4%) active in all six ( Fig. 2 C). The highest overlap of induced spectral power in more than two bands was observed for the electrodes in the posterior sensory visual areas of the occipital cortex ( Fig. 2 C; red), where 4-5 bands were activated on average per electrode during word encoding. Other visual areas of the temporal and parietal cortex also showed relatively high distributions of electrodes with more than two bands activated, in contrast to the anterior higher order processing areas in the temporal and prefrontal cortex, where electrodes were active in 1-2 frequency bands ( Fig. 2 C; black). Overall, the majority of electrode sites showed the induced spectral power in one or two specific frequency bands, suggesting that the various activities were spatially segregated in the cortex. Proportions of active electrodes in six frequency bands studied (low theta: 2-4 Hz, high theta: 5-9 Hz, alpha: 10-15 Hz) beta: 16-25 Hz, low gamma: 25-55 Hz, high gamma: 65-115 Hz) show significant differences (ANOVA) in the anatomical distribution (red asterisks; p < 0.05 post-hoc Tukey Kramer test) but not across the frequency bands as summarized in the inset box plots. (C) Distribution of the electrodes active in one or more of the frequency bands reveals the highest overlap in the posterior areas of the occipital and parietal visual cortex, which was gradually decreasing in the more anterior cortical areas as shown in the average brain surface plot and the cartoon summary of all active electrodes (dots colored according the bar chart) and in the violin plots (black and red lines indicate the mean and median, resp. ( Hoffmann, 2015 )). (D) Distribution of all electrodes active in the gamma (high frequency), alpha/beta (intermediate frequency), and theta (low frequency) bands differs across nine selected cortical regions of interest (V -visual; IT -inferior temporal; Pre -precuneus; Par -lateral parietal; MTL -mesial temporal lobe; LT -lateral temporal; Br -Broca's area; PFC -prefrontal cortex; FP -frontal pole) with relatively more prefrontal electrodes active in the gamma bands (red) and more electrodes in the visual areas active across all the bands (black) as shown in the average brain surface plot, the cartoon summary, and the violin plot. Notice that most electrodes overall were active in just one or two frequency bands either of the low, intermediate, or high frequencies, except for the visual areas with more than two frequency bands active per any one electrode (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).
We next asked whether the spectral activities that were induced at any one electrode site overlapped in similar frequency ranges. We grouped the activities into low, intermediate, and high frequency groups (theta, alpha/beta, and gamma, respectively) and compared proportions of active electrodes with different group combinations. More than half of all active electrodes showed induced power in exclusively the low, the intermediate, or the high frequency bands ( Fig. 2 D). Electrodes with high overlap across all groups were predominantly localized in the posterior visual areas ( Fig. 2 D; black). The remaining electrode sites that showed induced power in one of the three frequency groups were uniformly distributed across the cortex with similar proportions in each ROI, apart from the prefrontal areas where relatively more high frequency gamma activities were present during memory encoding ( Fig. 2 D; red). Our results showed that activities of a particular frequency range were induced mostly at specific cortical sites, forming a mosaic-like pattern of distribution consistent with the notion of "spectral fingerprints" (Siegel et al., 2012) .

Laterality effect of the induced power is specific to particular frequencies and cortical areas
Our next question was if the total power-in-band induced during the time of word encoding ( Fig. 1 ) was different depending on the anatomical location, the frequency band, or the hemisphere. Individual active electrodes can show similar magnitudes of the induced power, despite the heterogeneous distribution across the cortical regions ( Fig. 2 ). Overall, there was a significant effect of the frequency band on the induced power (repeated measures ANOVA F = 1346.99, p < 0.001, df = 5). Brain region (ANOVA F = 38.21, p < 0.001, df = 37) and of the hemisphere location of the electrode (ANOVA F = 17.4, p < 0.001, df = 1). Post-hoc Fig. 3. Temporal profiles of induced power during encoding reveal laterality effects in specific frequencies and cortical areas. (A) Each panel shows trial-averaged power changes across the time of word presentation (gray background) estimated from all active electrodes in a given frequency band and brain region plotted separately for the left and right hemisphere (red and blue, resp.). The gray area plots in the background summarize the statistical difference (t-statistic) between the two hemispheres in 50 ms time bins with dark gray indicating significantly more (positive t values) or less (negative t values) power in the left hemisphere (horizontal dashed lines; Student t -test, p < 0.05). Notice that the significant laterality effect localized in specific frequency bands in a given brain region. (B) Summary comparison of the total induced power (repeated measures ANOVA) from all frequency bands across the nine brain regions (see Fig. 2 for labels) and four stages of memory encoding -the within subjects factor (red asterisks indicate significant difference; Tukey-Kramer test, p < 0.05). (C) Memory effect was the largest in the left prefrontal cortex (indicated by black double arrow), whereas, MTL showed the largest effect in the right hemisphere areas (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.). group analysis (Tukey-Kramer, p < 0.05) confirmed significantly smaller power with increasing frequencies of the bands, except between the low and high gamma bands ( Supplementary Fig. 1). Electrodes within the occipital lobe showed significantly higher induced power than in any other area in agreement with our previous study of the high gamma activities ( Kucewicz et al., 2019 ).
In general, the total induced power was significantly different between the left and the right hemisphere. We investigated this laterality effect by estimating the temporal profiles of the mean power change in each ROI. While there was a significant laterality effect in at least one frequency band of every ROI ( Fig. 3 A; t-statistic, p < 0.05), the timing and magnitude of this effect was specific to a particular band and ROI. The temporal profiles of the induced power in the left and the right hemisphere ( Fig. 3 A; red and blue, resp.) provide another spectral signature of particular cortical areas. Hemisphere had an effect on the induced power in time ( Fig. 3 B; repeated measures ANOVA, encoding phase as within subjects factor, F = 17.01, p < 0.001, df = 3). When the total induced power was compared between the brain regions (posthoc Tukey-Kramer, p < 0.05) only the prefrontal, inferior temporal, and visual cortex showed a significant laterality effect ( Fig. 2 C). The effect in these three areas was significant regardless of the frequency band analyzed. Significantly more power in the left prefrontal and inferior temporal areas was expected given the role of these areas in processing verbal information.

Low, intermediate and high frequency activities show distinct spatiotemporal dynamics
In agreement with the previous studies ( Burke et al., 2014 a; Burke et al., 2013 ;Kucewicz et al., 2014Kucewicz et al., , 2019, the spectral power induced in this task followed a temporal sequence of brain regions, showing enhanced or suppressed power at specific times of word encoding ( Fig. 4 ). We compared these sequences across the nine ROI in the six frequency bands with reference to the four phases of word encoding (PRE -before the word appears on the screen, EARLY -the first 800 ms of the word presentation, LATE -the second 800 ms the presentation, and POST -after the word disappears from the screen). First, we confirmed, while accounting for ROI, hemisphere, and frequency band, that the spectral power was significantly modulated by the phase of memory encoding (repeated measures ANOVA, F = 729.45, p < 001, df = 3). In each band, we found that the power alternates between a relative induction and a suppression, but the timing and the ROI sequences were different across the frequency spectrum, as revealed by the different orders of ROI labels in the y-axis of Fig. 4 A. The low (theta), intermediate (alpha/beta), and the high (gamma) frequency activities each showed a distinct pattern ( Fig. 4 A). The intermediate frequency (alpha and beta) power was enhanced mainly in the pre-and the post-encoding phases. In contrast, the low and the high frequency band power was induced in the early encoding phase in two sequences of activation -first at word onset in the theta bands and then following the onset in the gamma frequencies. We saw similar patterns for forgotten word trials, but at lower magnitudes of power change ( Supplementary Fig. 2). There was a significant effect of the frequency band (repeated measures ANOVA F = 18.06, p < 0.001, df = 5), of the brain region (ANOVA F = 15.57, p < 0.001, df = 37), and of the hemisphere location (ANOVA F = 72.72, p < 0.001, df = 1) on the peak power latency. In the gamma bands, earlier latencies were observed in the visual compared to the prefrontal cortical areas (post-hoc Tukey-Kramer, p < 0.05), which was not the case for the theta bands ( Fig. 4 B). The high gamma power peaks occurred significantly later than the theta peaks (post-hoc Tukey-Kramer, p < 0.05) in the prefrontal cortical areas ( Fig. 4 C left) but not in the visual areas ( Fig. 4 C  right), where the alpha/beta peaks followed later than the high gamma peaks. Therefore, the sequences of the induced power were different, both between the low and high frequency activities in the same cortical areas and between the cortical areas in the same frequency bands. On a more granular level, we investigated the peak latency time of each individual electrode plotted on brain surfaces ( Fig. 5 A) to confirm the general posterior-to-anterior sequence (see Fig. 4 A), particularly in the gamma frequency bands. This general pattern only applied to a subset of electrodes in a given cortical area. The mosaic-like pattern (see Fig. 2 ) of various peak latencies was observed across the cortex. Significant correlations between the peak latency and anatomical position (Pearson's correlation) were predominantly observed in the anterior-posterior axis ( Fig. 5 B), providing additional evidence for the anatomical sequences of the induced spectral activities. This anatomical sequence, as well as the mosaic-like pattern of individual electrode activations, was most evident in the low theta and the high gamma bands (Video 1).
Two sequences of induced power during successfully encoded words appear "moving " across the cortical space. Time bar at the top and bottom changes from black to red to indicate word display.
Finally, to provide a complete picture of these spatiotemporal dynamics across the frequency spectrum we interpolated the induced power values from all active electrodes on an average brain surface. Spectral power was predominantly induced in the posterior visual and anterior prefrontal cortical areas at the encoding phases specific for the low, intermediate, and high frequency activities ( Fig. 6 top). Induced power was first observed in the theta activities of the visual areas before word onset and then right after the onset in the anterior prefrontal areas. This early posterior-to-anterior sequence of theta power induction was followed by a secondary induction of the high frequency gamma power in the same sequence later into the word encoding. These two sequences of induced power are "moving " across the anatomical space in both hemispheres (Video 1). There was no analogous sequence observed in the alpha or the beta bands (Suppl. Video 1). The alpha/beta power was induced in the same posterior and anterior areas but predominantly before and after word presentation, suggesting a more preparatory or inhibitory role in our proposed model of human memory encoding ( Fig. 6  bottom). During preparation for encoding of the incoming words, the intermediate frequency oscillations would prevail and be followed by induction of the low frequency power in the theta bands across the brain at the time of word onset. The theta rhythms would in turn entrain the posterior and the anterior areas for final induction of the high frequency Fig. 4. Low, intermediate, and high frequency activities are temporally arranged into distinct sequences of cortical brain regions. (A) Mean power changes estimated from correct recall trials of all active electrodes in a given brain region of the left hemisphere are ordered from the earliest to the latest peak latency (dashed lines separate different phases of memory encoding). Notice consistent sequences of activation in the theta, the alpha/beta, and the gamma frequency bands. Areas listed on the y-axes are sorted in order of the earliest peak and thus are specific to each frequency band. (B) Summary of the mean latency ranges (post hoc ANOVA mean comparison) in the time of word presentation (0 corresponds to onset) shows significantly earlier gamma power induction in the most anterior FP areas (black) relative to the more posterior areas (red; p < 0.05, Tukey-Kramer test), and an opposite pattern of the theta power induced earlier in the posterior V area relative to the more anterior brain regions (see Fig. 2 for labels). (C) Anterior and posterior brain regions show analogous patterns of power induced significantly earlier in the low frequencies of the prefrontal areas on the left, and significantly later in the alpha/beta frequencies of the visual areas on the right (red; p < 0.05 Tukey-Kramer test) relative to the high gamma band (black) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.). gamma activities during the time of word encoding, first in the visual, and then progressively later in the higher-order prefrontal cortical areas ( Burke et al., 2014 a;Kucewicz et al., 2019 ). The model is proposed to link the dynamics of particular spectral activities with the associated cognitive processes involved in anticipation, attention, and encoding of the remembered stimuli.

Discussion
Separate anatomical and temporal distributions of the low, intermediate, and high frequency spectral activities in the human cortex are the main findings of this work, along with the resultant simplified holistic model of the spectral power dynamics associated with human memory encoding. Because of the large number of participants and electrodes, implanted throughout all cortical regions, recording from discrete and spatially distinct sites across 39 Brodmann areas was possible. With exception of the visual cortex, within any one cortical area a majority of the electrode sites showed task-induced activity only in one or two bands of the frequency spectrum. At the same time, activities in all six bands were found among all sites within a given area, producing a "mosaiclike " pattern of theta, alpha, beta, and gamma activities recorded at specific cortical sites. Each site could be viewed as one small tile of the mosaic with a specific color corresponding to a particular band as visualized in the Fig. 2 brain surface plots. Our results are in agreement with the spectral fingerprint view (Siegel et al., 2012) that processing in a given cortical site is characterized by specific neural activities, as previously demonstrated in the MEG and iEEG studies ( Fellner et al., 2019 ;Keitel and Gross, 2016 ). Various spectral activities were sparsely distributed in any one cortical region in more or less equal proportions, despite subtle differences, e.g., relatively more gamma band activities in the prefrontal cortex. This general pattern of heterogeneous, sparse distribution would explain the observation of a broadband "spectral tilt" of power across all frequencies Voytek and Knight 2015 ;Kilner et al., 2005 ;Miller et al., 2014 ;Herweg et al., 2020 ) when activities from multiple sites in a given brain area were averaged together. As a result, low frequency power ( < 30 Hz) is typically decreased and high frequency power ( > 30 Hz) is increased upon activation of a cortical area in a given task. This broadband tilt in the spectral power can be resolved into specific spectral fingerprints at the level of individual electrode sites. One could go even further to ask how resolved this spectral distribution is and whether it can still be observed on the level of micro-electrode sites. Our study was limited to multiple macro-contacts that were grouped from different subjects -high-density recordings from multiple macro-and micro-electrodes implanted in individual subjects (Kucewicz, Berry, Worrell chapter in ( Lhatoo et al., 2019 )) would be required to address these questions further and elucidate the neuronal activities underlying our reported mosaic of spectral fingerprints.
The six frequency bands were similarly represented in the cortex, as quantified by the relative proportion of electrode sites recording theta, alpha/beta, and gamma activities. Even though cortical task activation was typically associated with spectral changes in the theta and gamma frequency bands ( Miller et al., 2014 ;Greenberg et al., 2015 ;Osipova et al., 2006 ;Solomon et al., 2017 ;Burke et al., 2013 ;Kucewicz et al., 2014 ), we found the same or higher number of active electrode sites in the alpha and beta bands. Oscillations in these frequency bands were previously implicated with anticipatory or even inhibitory processes preceding task activations ( Spitzer and Haegens, 2017 ;Engel and Fries, 2010 ). Compared with post-stimulus gamma activities related to sensory processing, alpha and beta activities were reported before stimulus onset and modulated by attention ( van Ede et al., 2014 ;Bauer et al., 2014 ). Our results confirm the different timing of alpha and beta power induction, occurring mainly before and after the period of word presentation on the screen. These were observed in anatomical areas that were overlapping with the gamma and theta power induced during word presentation. Hence, the spatiotem-poral pattern of spectral power confirmed distinct roles played by the theta and gamma, as well as the alpha and beta activities.
Theta frequency bands were also found to show a profile of spatiotemporal changes that was distinct from the gamma activities. Taskinduced power in the high gamma band is known to arrange into a hierarchical sequence of brain regions activated first in the posterior sensory areas and then in the more anterior associational temporal and prefrontal cortex (Kucewicz et al., 2019(Kucewicz et al., , 2014. This posterior-to-anterior sequence of high gamma power inspired a two-stage model of memory encoding ( Burke et al., 2014 a) to distinguish the early sensory and the late semantic phases of processing words. Theta power was also reported in the same areas preceding the gamma activation ( Burke et al., 2013 ). Our recent analysis showed that this hierarchical activation from the sensory posterior areas to the anterior semantic systems propagates continuously along the ventral visual and semantic processing streams, ending in the ventrolateral prefrontal cortex of the Broca's speech area and in the frontal pole ( Kucewicz et al., 2019 ). Here, we observed a "mosaic-like " pattern of induced power observed across the cortex in the posterior-to-anterior sequence that was different for the theta and the gamma activities (Video 1). First, a temporal sequence of theta peaks was induced at the time of stimulus onset, in contrast to the sequence of gamma peaks following the onset. In addition, the shortest latencies of the induced theta power peaks were localized in the Broca's and prefrontal areas even before the visual cortex. The gamma power peaks, on the other hand, revealed the shortest latencies in the visual cortex and the longest in the Broca's area, prefrontal cortex and the frontal pole. Hence, both the timing and the anatomical sequence of cortical areas was different for the power induced in the theta and the gamma bands. This suggests that the low and high frequency activities are independently induced, playing different roles in stimulus processing. It remains an open question whether this sequence of peak power in our study is in any way related to the actual travelling waves evident in the spectral phase analysis (not power) of theta and alpha oscillations ( Zhang et al., 2018 ) , which confirmed the general posterior-to-anterior directionality of the local, cortical waves. That study focused on more local discrete cortical locations of either slow theta, high theta, or alpha phase propagation, which is different from more global sequences of power across a broad frequency range studied here. Another possible way of analyzing the directionality of spectral amplitude or power changes would be to use cross-correlations between pairs of electrode sites ( Adhikari et al., 2010 ). Here, we focused on the independent mosaic-like pattern of sequential power inductions in different frequency bands for the benefit of a simplified holistic model.
We summarized this holistic picture of the low, intermediate, and high frequency neural activities in a simplified model of the relative spatiotemporal dynamics during memory encoding (see Fig. 6 ). Alpha and beta oscillations dominated in the initial pre-stimulus presentation phase potentially related to anticipatory and attentional processes in both the posterior sensory and the anterior association areas. Theta activities are then induced in these areas, we propose, in preparation for the expected stimulus processing around the time of presentation. Gamma activities are finally induced in response to the stimulus presentation in a sequence of the processing stream from the visual cortex to the prefrontal areas. This comprehensive view of the relative dynamics of these activities can be best captured using all dimensions of the spectral scale, anatomical space and time of stimulus processing (see Suppl. Video 1). Functional roles of these distinct spectral activations are not the subject of this study, although they are congruent with anticipatory, attentional, preparatory, and perceptual processes proposed in the literature. In general, the lower frequency spectrum dominates in early phases of preparation and receiving stimulus information, whereas the late phase of encoding the presented information is characterized by gamma activities when the low frequency power is back at baseline or even suppressed. The observation that the low and high frequency activities were not induced at the same time would not need to preclude a possible cross-frequency interac-tion. For instance, bursts of induced gamma power that were previously described in humans and non-human primates ( Kucewicz et al., 2014Lundqvist et al., 2016Lundqvist et al., , 2018, understood in frames of various models like the theta-gamma memory buffer ( Lisman and Jensen, 2013 ), still would be reconciled with our reported independent spatiotemporal dynamics. Our results show no evidence for amplitude-amplitude interactions. However, phase-phase or phase-amplitude interactions between the low and the high frequency activities ( Canolty et al., 2006 ) are still possible and could explain the capacity limit of the remembered items ( Kami ń ski et al., 2011 ), which was actually observed also in our study of free recall. In conclusion, our model shows that spectral power in the three frequency ranges is induced independently in anatomical space and time of memory formation, but this does not preclude phase interactions between activities in these ranges, which were beyond the scope of this study.
Neural dynamics differed between the two hemispheres, revealing different profiles of spectral power, particularly in the prefrontal and the temporal cortex. The lateral prefrontal and inferotemporal cortices showed the greatest laterality effect and SME with more power in the left hemisphere, while the visual cortex had more power in the right hemisphere. In general, this hemispheric asymmetry was found to various degrees and in at least one frequency band in each brain region. Spectral fingerprints were thus specific to the hemisphere. One possible reason for this asymmetry is a differential engagement in the encoding and retrieval processes. According to the HERA (hemispheric encoding/retrieval asymmetry) model, we would expect brain activities during encoding to be stronger in the left hemisphere than in the right hemisphere and for the opposite to be true for retrieval ( Buckner et al., 1996 ;Rugg et al., 1996 ;Fletcher et al., 1998Fletcher et al., b, 1998Grady et al., 1998) . Therefore, it is not surprising that higher induced power was found in the left lateral prefrontal and inferior temporal cortex, especially in a verbal task with words. While a complimentary analysis of the induced power during retrieval of recalled words was beyond the scope of this study, others reported in the recall epoch of the same task more theta power in the right temporal and more gamma power in the left prefrontal and temporal cortex areas ( Burke et al., 2014 a). Previous studies with different tasks ( Tulving et al., 1994 ;Fletcher et al., 1997 ) associated the right prefrontal cortex with episodic memory retrieval, and the left prefrontal cortex with encoding. Therefore, the preferential activation of the left prefrontal cortex during encoding in our task can be explained both within the HERA model and also by the verbalizability of the stimuli used ( Golby, 2001 ). The latter could be further studied with determination of language dominance for every subject, which was not possible in our study with information available for only 4 subjects (Supplementary Table 1). Obtaining and sharing this information in large, multi-center projects is challenging, and future studies will have to address this limitation. Nevertheless, in terms of the hypotheses tested in our study, these effects provide further support for independent spectral activities in the two hemispheres of particular cortical areas.
Influence of the pathophysiology of epilepsy on the laterality effect or the estimates of spectral power in general can only be minimized but not completely eliminated in this patient population. Electrode channels that showed epileptiform activities in the seizure-generating areas were assumed to be excluded from the analysis during the automated selection of active electrodes -an added benefit of the selection method . Although occurrence of epileptiform activities was shown to be related with memory processing at distinct phases and anatomical locations ( Matsumoto et al., 2013 ;Horak et al., 2017 ), it is highly unlikely that these pathological activities would consistently occur at specific phases of memory encoding and thus bias the power estimates used for selecting active electrodes. Hence, the laterality or the memory effects reported here are presumed to be only minimally affected.
Our results showed that it is not enough to generalize the neural activities to a given cortical area but a further distinction between the area in the left and the right hemisphere is needed, at least in cases of verbal memory tasks. In our verbal task, the cortical areas with the greatest laterality effect were also the ones to reveal a high subsequent memory effect. As discussed previously, the strong laterality effect may be due to language lateralization in the left hemisphere, as would be expected in the majority of the general population. In the left temporal and prefrontal cortex, spectral power induced on trials with words that were subsequently recalled was greater than on the ones with words that were forgotten. This subsequent memory effect was gradually increasing in magnitude with successively more anterior cortical areas. Previous neuroimaging and electrophysiological studies found this memory effect to be the highest in the lateral prefrontal and the temporal cortices during similar, verbal memory tasks ( Wagner et al., 1998 ;Long et al., 2010Long et al., , 2014Kucewicz et al., 2019 ;Kim, 2011 ;Burke et al., 2013 ;Sederberg et al., 2003 ). Our recent investigation of the memory effect in the same task localized the highest magnitude in the left ventrolateral prefrontal cortex close to the Broca's area, and in the left occipitotemporal cortex ( Kucewicz et al., 2019 ). Altogether, the cortical areas identified with the greatest induced power and hemispheric laterality in the task overlap with localization of the memory effect in this verbal task. This pattern may be specific to verbal memory tasks, suggesting it is driven by the lateralization of language processing but not necessarily memory formation.
Despite its known involvement in declarative memory function ( Beason-Held et al., 1999 ;Zola-Morgan et al., 1994 ;Parkinson et al., 1988 ;Squire and Zola-Morgan, 1991 ), MTL and the hippocampus had a relatively low number of active electrodes, magnitude of the induced spectral power, and laterality effect. One reason for this finding is a sparser distribution of hippocampal activities than in the neocortex, which could result in effectively less induced power recorded on any one macro-contact electrode. Sparse micro-contact neural activities would be averaged out by the surrounding non-active areas, compared to more topographically mapped cortical areas. Another reason is that the task itself may engage the MTL less than a more spatial or episodic memory task. Remembering the sequence of words on the list would be expected to increase hippocampal involvement in the task. A greater neocortical contribution to this task would explain our reports of enhanced performance in this task with stimulation in the lateral temporal cortex Ezzyat et al., 2018 ). Stimulation in MTL was Fig. 6. Anatomical spread of low and high frequency activities reveals independent posterior-to-anterior temporal sequences of spectral power. (Top) Average brain surface plots summarize interpolated power from all active electrodes at selected four timepoints of the pre-, early, late, and post-word encoding periods (gray background indicates word presentation on the screen). Notice analogous patterns of the anatomical spread for the low (theta), intermediate (alpha/beta), and high (gamma) frequency activities, especially in posterior visual and the anterior prefrontal areas. (Bottom) Model of independent activation of the low and high frequency power "moving " from the posterior to the anterior cortical areas during word presentation, preceded and followed by the intermediate frequency alpha and beta band activation outside of the word encoding period. Red arrows indicate directions of the sequences (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.). found to have an opposite effect ( Kucewicz et al., 2018 b;Jacobs et al., 2016 ), suggesting different roles for the two structures in this task. And given reports of the contrary, where performance improved with MTL stimulation ( Suthana and Fried, 2014 ;Fell et al., 2013 ), the roles and outcomes may differ with subtle changes within the structure as well.
Mapping neural activities that are critical for memory encoding guides the development and targeting of brain modulation, e.g. with the direct electrical stimulation, for therapeutic and research purposes. Out of the multitude of spectral activities that are induced at various times and anatomical locations it is necessary to identify prospective targets for modulation. Our holistic map and model of independent spectral activities and their various derivatives, including the laterality or the memory effect, offer simple and accessible biomarkers for new braincomputer interface technologies. Such biomarkers can be rapidly computed to map the brain regions, times and states for modulating brain activities underlying memory and other cognitive functions. For instance, they can be used to determine the neural activities correlated with memory enhancement ( Kucewicz et al., 2018 b), to predict memory states for brain stimulation and rescue poor encoding trials ( Ezzyat et al., 2017 ), or trigger modulation online in a closed-loop design for responsive electrical stimulation . Hence, the spectral biomarkers can be useful for developing and optimizing the emerging technologies and new therapies for mapping and restoring memory functions.
VSM wrote the manuscript and analyzed the results. KVS preprocessed the data and developed the methodology. CT contributed to data processing and analysis. ML and TPT processed and analyzed the brain mapping imaging data. PN and VK contributed to data processing and analysis. GAW designed the study. MTK designed the study and wrote the manuscript. All authors contributed to manuscript edition and data analysis and/or interpretation.

Data availability statement
This work was supported by the Defense Advanced Research Project Agency grant called: "Restoring Active Memory (RAM) " under Cooperative Agreement N66001-14-2-4032, as part of the BRAIN initiative (Brain Research through Advancing Innovative Neurothechnologies). The raw data used for this analysis can be requested on the project website (http://memory.psych.upenn.edu/RAM). The views, opinions, and/or findings contained in this material are those of the authors and should not be interpreted as representing the official views or policies of the Department of Defense or the U.S. Government. Code used for identifying active electrodes can be found at https://github.com/kvsaboo/TaskActiveElectrodeIdentification. Violin plots generated using code originally written by Holger Hoffmann and available on MATLAB's file exchange (https://www.mathworks.com/matlabcentral/fileexchange/45,134violin-plot). Other code available upon request.
Video 1. An early posterior-to-anterior sequence of low theta power induction is followed by a secondary induction of the high gamma power during word encoding.

Acknowledgments
Cindy Nelson and Karla Crockett assisted in participant recruitment and data collection. This work would not be possible without the dedicated effort and participation of participants and their families. Openaccess datasets were originally collected as part of a BRAIN Initiative project called Restoring Active Memory (RAM) funded by the Defense Advanced Research Project Agency (DARPA; Cooperative Agreement N66001-14-2-4032). This research was supported from the First Team grant of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund (Grant No. POIR.04.04.00-00-4379/17 ).