Planar, spiral, and concentric traveling waves distinguish behavioral states in human memory

Domain: neuroscience
Published: 2026-04-11
Editorial outcome: minor_revision

Abstract
A fundamental challenge in neuroscience is explaining how widespread brain regions flexibly
interact to support behaviors. We hypothesize that traveling waves of oscillations are a key
mechanism of neural coordination, such that they propagate across the cortex in distinctive
spatial patterns that control how different regions interact. To test this hypothesis, we used direct
brain recordings from humans performing multiple memory experiments and an analytical
framework that flexibly measures the propagation patterns of traveling waves. We found that
traveling waves propagated along the cortex in not only plane waves, but also spirals, sources
and sinks, and more complex patterns. The propagation patterns of traveling waves correlated
with novel aspects of behavior, with specific spatial patterns reflecting particular cortical
processes and even individual remembered items. Our findings show that large-scale cortical
patterns of traveling waves reveal a diverse range of spatial patterns in the brain that are relevant
for neural decoding.

Introduction
Neurons in the human cortex have rich dendritic arbors that integrate diverse inputs and axons
that project outputs to multiple, distributed areas (Swanson, 2003). However, recent
developments in brain imaging have suggested a dynamic interplay of distributed brain regions
that underlie complex human behaviors (Sporns et al., 2004). How do individual neurons or
regions reorganize their activity so that they are selectively and dynamically linked to certain
other areas even as the structure of these neurons does not change on the timescale of behavior?
We hypothesize that propagating patterns of brain oscillations, or “traveling waves,” underlie
this selective reorganization. Traveling waves propagate across the cortex in specific directions
during behavior, so that certain spatial arrangements of neurons are active at similar times, thus
flexibly and efficiently communicating between the given brain regions that are relevant for a
task (Mohan et al., 2024). Traveling waves are present at various frequencies (Ermentrout &
Kleinfeld, 2001), including the theta/alpha (~ 4-12 Hz) range in humans (Zhang et al., 2018).
Because propagating oscillations correlate with underlying neuronal activity (Jacobs et al.,
2007), the spatiotemporal organization of traveling waves at each moment may indicate which
brain areas are active and how they correspond to particular computational processes and
memory representations in the cortex (Pinotsis et al., 2023; Pinotsis & Miller, 2023). The
propagation of traveling waves may therefore show how the brain’s spatial organization and
connectivity flexibly adapt to task demands and behaviors.

Traveling waves are known to play a critical role for behaviors such as visual processing and
spatial navigation in rodents and non-human primates (Agarwal et al., 2014; Davis et al., 2021;
Lubenov & Siapas, 2009; Muller et al., 2014; Zanos et al., 2015), and recently in human
cognition (Alamia & VanRullen, 2019; Kleen et al., 2021; Mahjoory et al., 2020; Stolk et al.,
2019; Zhang et al., 2018). Traveling waves also appear in spirals during sleep in humans (Muller
et al., 2016). However, our current understanding of traveling waves in cognition is largely
focused on planar waves, where a consistent direction of wave propagation is maintained across
a large region of cortex. It is possible that other types of traveling waves with more complex
spatial patterns are also relevant functionally (Bhattacharya, Brincat, et al., 2022; Denker et al.,
2018; Freeman, 2009; Muller et al., 2016), so this focus on plane waves in humans may have
limited our understanding of the role of traveling waves and oscillations in modulating cortical
representations to support behavior.

We recently found that the direction of propagation by planar traveling waves during behavior
correlated with distinct memory processes (Mohan et al., 2024). Extending this idea, we tested
here the broader idea that the spatial patterning of wave propagation changed shape according to
the current behavior.We designed a flexible analytical framework for measuring general patterns
of traveling wave propagation and applied this procedure to direct brain recordings from
neurosurgical patients performing multiple memory experiments. Our results showed an array of
complex traveling wave propagation patterns including spirals, concentric waves (sources and
sinks) (also known as target waves (Hagan, 1981; Xu et al., 2015; Zhang et al., 2003)), and
heterogeneous directional propagation patterns that extended beyond those seen previously in
humans.

Results
Theta/alpha and beta frequency oscillations are widespread and organized as spatiotemporally
stable traveling waves

To probe the role of brain oscillations and traveling waves in organizing the spatial and temporal
structure of cortical activity during memory, we examined human electrocorticographic (ECoG)
recordings from surgical epilepsy patients as they performed spatial and verbal memory tasks
(Methods, Figure 1, Supplementary Figure 1). We sought to identify brain oscillations that were
spatially organized into traveling waves and test whether they reorganize into different
directional patterns to distinguish separate behaviors. To examine this hypothesis, we developed
a flexible framework for measuring traveling waves in ECoG recordings, quantifying their
instantaneous spatial structure, and identifying spatial patterns that differentiate individual
behavior.

To flexibly detect spatial patterns of traveling waves in each subject, we first identified the
frequencies where groups of contiguous electrodes showed common oscillations (Methods,
Figures 1C–E). Identifying a single common frequency is crucial because, by definition, a
traveling wave involves an oscillation at a single frequency that progressively propagates across
a region of cortex, thus making it possible to detect the traveling wave when it passes by these
electrodes (Das et al., 2022). Overall, oscillations were most often present in the theta (~ 4-8
Hz), alpha (~ 8-12 Hz), and beta (~ 12-30 Hz) frequency bands (Figure 1C, Supplementary
Tables 1, 2), and ~86% of all electrodes showed a narrowband oscillation in at least one of these
bands. Interestingly, sometimes the same ECoG grid showed multiple peaks in the power
spectrum, i.e., they had peaks at both the theta/alpha and beta frequency bands (Supplementary
Tables 1, 2, also see Figure 8), which suggests the presence of multiple overlapping oscillatory
networks.

We next distinguished the patterns of oscillations that were traveling waves. Because a traveling
wave involves an oscillation at a single frequency that is present simultaneously on multiple
electrodes, we identified the contiguous clusters of electrodes with oscillations at the same
frequency. A traveling wave requires the presence of a systematic propagation of phase across a
cluster of electrodes. To test for this pattern, we filtered the signals of the electrodes in each
cluster at their peak frequencies and extracted the instantaneous phase at each electrode and
timepoint using the Hilbert transform (Figure 1F). Visual observation of the propagation of the
absolute phase (Hilbert-transformed phase) across time showed clear traveling waves, with phase
advancing in consistent spatial patterns across neighboring electrodes (see Figure 1H, for
example).

Next, to probe this phenomenon systematically, we implemented a quantitative framework for
measuring the spatial pattern of each traveling wave. Here, for each electrode in a cluster, we
used a windowed circular–linear regression approach to identify the relative phases that exhibit a
linear relationship with electrode locations locally (Das et al., 2022) (Methods). This regression
also measured the properties of each wave’s propagation, including its phase velocity and
direction at each timepoint (Figure 1H, Methods).

Recent evidence has shown that oscillations and traveling waves in humans often appear in
transient bursts, but that they nonetheless are functionally relevant (Freeman & Rogers, 2002;
Muller et al., 2018; Roberts et al., 2019; Schmidt et al., 2023; van Ede et al., 2018; van Vugt et
al., 2007). We thus designed a single-trial method that identified stable periods of wave
propagation (Methods, Figure 1G) and used this procedure to identify individual epochs where
waves progressed across the cortex in a consistent spatial arrangement (Methods, Figures 1G, H,
Supplementary Videos 1-7). This procedure revealed that the stable periods of individual
traveling waves lasted ~80–180 msec on average (Figures 1EG, H). These wave epochs were
statistically robust (all p’s < 0.001, Methods, Figures 1IR, J, Supplementary Figure 2,
Supplementary Tables 3, 4) ruling out the possibility that these stable epochs could arise due to
chance or noise, and their period of stability was shorter for oscillations at high frequencies than
low (median 107 versus 134 msec, p = 0.002, Mann-Whitney U-test, Supplementary Table 5).
Consistent with earlier work (Das et al., 2022; Halgren et al., 2019; Stolk et al., 2019; Zhang et
al., 2018), the propagation speed of these traveling waves was ~0.3-3.3 m/s, and faster for
oscillations at higher temporal frequencies than low (median 1.1 m/s for 12-26 Hz versus 0.5 m/s
at 5-12 Hz, p < 0.001, Mann-Whitney U-test; (Zhang et al., 2018)).

Independent component analysis (ICA) reveals the most dominant traveling wave patterns

Using this approach, we observed a diverse range of robust traveling wave patterns across
subjects (Figures 2, 3). These patterns included plane waves as seen previously, but also spirals,
concentric waves, and spatially heterogeneous, complex waves (for examples, see
Supplementary Figures 9-18). We also saw that, over time, the traveling waves at individual
electrodes often shifted between different patterns. To quantitatively distinguish the full diversity
of spatial wave patterns over time, we used an algorithm based on independent component
analysis (ICA) (Fu et al., 2015; Li & Adalı, 2010) to identify the range of spatial patterns of
traveling waves that appeared on each electrode cluster (see Methods). Using this algorithm,
we labeled each category of spatial wave pattern that appeared on each cluster, which we refer to as
a “mode” (Figure 2, Methods). We classified each mode according to its shape, including plane,
spiral, and concentric waves, as well as other waves with complex shapes (Figure 3). This
algorithm provided a series of weights that quantify the contribution of each mode to the epoch’s
current shape. We call these weights the epoch’s “activation function”, and they quantify the
instantaneous magnitude and direction of a mode for each epoch. Thus, by examining the modes
and activation functions, we quantitatively showed the types of traveling wave patterns that were
most strongly present at each moment in the recording, with the magnitude of the activation
function indicating how strongly each mode was represented.

In addition to applying our methods on the experimental data, we also carried out extensive
simulations to statistically validate our proposed ICA approach to confirm that it was working as
expected in terms of sensitivity and specificity for extracting true traveling wave patterns and
demonstrating negative results in the absence of such patterns. In particular, we carried out four
different types of simulations, which validate our approach (see Methods for details) and identify
robustness against noise, the ability of ICA to extract multiple types of traveling waves, and
show negative results in the absence of true traveling waves (Supplementary Figures 4, 5).

This ICA algorithm effectively quantified the spatial wave patterns that appeared visually in the
recordings, as the dominant wave patterns that were evident from visual inspection generally
corresponded with the modes with the strongest weights from ICA. For example, Figure 2D
shows three different wave patterns seen in the raw data, which match the corresponding
activation functions for the top 6 modes shown in Figure 2Ei, and the corresponding modes
shown in Figure 2Eii (also see Supplementary Figures 7, 8).

Using this ICA procedure to identify the modes that were present on each cluster throughout the
task, we identified a wide range of spatial traveling wave patterns across all clusters (Methods,
Figure 2, Supplementary Figures 9-18). Individual clusters showed a mean of 6 significant
modes, indicating that there were generally multiple spatial patterns of traveling waves within
individual clusters. On average the first three modes explained ~34%, ~19%, and ~10% variance
respectively, and the first six modes combined explained >80% variance (Figure 4E). Thus, ICA,
can identify the diverse traveling wave patterns that are present at each oscillation cluster, and
reveal how the spatial pattern of wave propagation varies throughout the task.

Diverse spatial patterns of traveling waves and their characteristics

To explain the diversity of the spatial wave patterns in the data, we classified each mode based
on its shape into one of the following categories: “planar”, “rotational”, “concentric”
(“expanding” or “contracting”), or “complex” (Methods, Figure 3) (Townsend et al., 2015).
Complex waves were those that showed a robust spatial propagation pattern (i.e., p’s < 0.001)
but did not meet the criteria for the other categories. Some complex waves showed a
combination of multiple patterns, such as separate subsets of electrodes showed planar, rotating,
or expanding/contracting waves (Figure 9, Supplementary Figures 9-18).

One previous study showed rotational and expanding traveling waves during human sleep
spindles, supporting the idea that these types of waves may support memory consolidation and
plasticity (Muller et al., 2016). Based on this, we hypothesized that rotational and expanding
(sources) traveling waves would be more prevalent in the spatial episodic memory task which
involves processing a broad range of inputs including spatial and non-spatial information.

The spatial organization of traveling waves varied across regions and temporal frequencies.
For clusters with oscillations at low temporal frequencies, plane waves were most dominant
(median wave strength 59% across tasks), and this pattern was present in both the spatial (2(3) = 165.0, p
< 0.001, Figure 4A) and verbal (2(3) = 288.7, p < 0.001, Figure 4C) memory tasks. At higher
frequencies, although plane waves were again most dominant overall (median wave strength
43% across tasks), there was an increased prevalence of rotational waves (median wave strength
37% across tasks), especially in the spatial memory task (2(3) = 162.3, p < 0.001, Figure 4A).
When we compared the strength of these different traveling wave patterns as indexed by variance
explained (see Methods), we found that planar waves were strongest, followed by rotating and
expanding/contracting waves (Figures 4E, F). Complex waves, although highly statistically
significant and strongly relevant to behavior (see below), were less strong and explained less
variance in the raw data compared to the other wave types.

We also compared these effects between hemispheres, in light of work on lateralized oscillations
in hippocampus and neocortex (Das et al., 2022; Miller et al., 2018). Although plane waves were
dominant in both hemispheres, they were stronger on the right hemisphere (~64%) than the left
(~42%; 2 test p < 0.01, Figures 4 B, D). We next compared this lateralization between the two
different memory tasks. In the left hemisphere, during the spatial memory task wave patterns
were preferentially rotational rather than concentric, (2(3) = 150.0, p < 0.001, Figure 4B).
Conversely, concentric wave patterns were more common in the verbal memory task (2(3) =
288.1, p < 0.001, Figure 4D). But, in the right hemisphere, we found statistically similar patterns
across both tasks. This suggests that the shape of traveling waves across the cortex during
behavior reveals specific task-related neural assemblies that are lateralized across the cortex.

Given that we saw both counter-clockwise and clockwise rotational wave patterns, we examined
whether the orientation of wave propagation was functionally relevant. One previous study
showed that during sleep, human cortical spindles are rotational traveling waves that propagated
in a temporal→parietal→frontal (TPF) direction (Muller et al., 2016). Therefore, we examined
whether rotating waves in our dataset had a directional preference by labeling each one’s
direction (TPF or TFP). If there were multiple rotational modes in a single oscillation cluster,
then we labelled each separately. We did not find a significant preference for a specific rotational
direction in either task. In the spatial memory task, there were 7 TPF waves versus 9 TFP waves
(p > 0.05, binomial test) and, similarly, in the verbal memory task (12 TPF versus 15 TFP, p >
0.05, binomial test; Supplementary Table 6).

We also compared the propagation of concentric traveling waves, comparing the prevalence of
inward (sink) versus outward (source) propagation. In the spatial task, there were more sources
compared to sinks, as 78% (7 of 9) concentric waves were sources. Inversely, in the verbal
memory task, there were more sinks than sources, as 68% (15 of 22) concentric waves were
sinks (Supplementary Table 6). Thus, the brain shows different types of concentric waves (2
test, p < 0.02) between spatial and verbal memory.

Traveling waves can distinguish both broad and specific behavioral states of human memory
representations

Next, to identify the potential functional role of traveling waves, we examined how the
prevalence of traveling waves with different shapes shifted between stages of memory. First, we
examined encoding, retrieval, navigation, etc. in the spatial memory task. To assess statistical
significance, we used multivariate analysis of variance (MANOVA) (Methods) to test for
changes in activation functions between conditions.

Traveling waves showed different propagation patterns between stages of memory processing
(Figure 5, Supplementary Figures 9-18, Supplementary Videos 1-4). As an example, Figures 5A-
C show a cluster of electrodes that showed three distinct modes of traveling waves that changed
their propagation direction between the stages of memory. In this cluster, the plane wave pattern
(mode 1) was nearly absent during navigation and distractor stages (as indicated by its low
magnitude) but was strongly present in other task phases. However, its direction of propagation
shifted between behaviors, with a different direction for retrieval compared to the encoding,
confidence, and feedback periods. Similarly, this cluster’s mode 2 also differed in propagation
direction between memory stages and had the highest wave strength during retrieval and
feedback. These state-specific traveling waves were prominent across our dataset, as all 13
oscillation clusters showed significant shifts in traveling wave direction or strength between
stages of the spatial memory task.

We also examined whether the directions and strength of traveling waves shifted to represent the
specific item that each subject was viewing, extending item-specific gamma oscillations seen
previously (Jacobs & Kahana, 2009). In the verbal memory task, we found that many clusters
showed traveling waves with different propagation directions and strengths for the encoding of
individual items. For example, Figures 6A-C show a cluster of electrodes with planar traveling
waves (on mode 1) that propagated consistently posteriorly when the subject viewed any letter
except for “G”. This same cluster also showed a direction-shifting traveling wave in mode 3
which propagated in an antero-superior direction for letters “G” and “Q” and the opposite
direction for letters “D” and “J”. These item-specific wave shifts were common, as 13 of 26
clusters in this task had traveling waves that shifted direction or strength for individual letters.
Thus, traveling waves often distinguish both broad and specific human memory representations
(Supplementary Figures 9-18).

Waves that distinguish behavioral states are of diverse patterns and widespread across the
human brain

Beyond plane waves, we were also interested in testing the behavioral significance of rotational,
concentric, and complex patterns of waves. Similar to plane waves, rotational waves were also
stable at the individual epoch level, distinguished both broad and specific behavior, and were
widespread across the frontal, temporal, and parietal lobes (Figures 7A-D, Supplementary Video
5, Supplementary Figures 9-18). These patterns were evident in both the spatial and verbal
memory tasks (Figures 7E, G, also see Supplementary Figures 9-18). For rotational waves, the
wave strength was higher during the distractor period compared to the other task periods,
especially in the left hemisphere (2(5) = 67.0, p’s < 0.001, Figure 7F).

Concentric and complex wave patterns also distinguished both broad and specific behavioral
states across brain regions and tasks (Figures 8, 9, Supplementary Videos 6, 7). Interestingly,
during the distractor phase of the task, the strength of complex waves was significantly higher
compared to other task periods across all frequencies and hemispheres (2(5) = 43.2, p’s < 0.001,
Figure 9F), however, this pattern was not seen for the concentric waves (2(5) < 19.6, p’s > 0.01,
Figure 8F).

Finally, we quantified the number of modes that were significant for each type of wave.
The percentage of modes that were significant in each category of wave did not differ from each
other, in any frequency or hemisphere (2(3) < 6.4, p’s > 0.05, Supplementary Figures 19A-D).
Moreover, the percentage of significance in each mode did not differ from each other (2(5) <
4.2, p’s > 0.05, Supplementary Figures 19E, F).

Behavioral states of human memory representations can be decoded from traveling waves

Since the spatial patterns of traveling waves varied reliably with behavioral states, we
hypothesized that we would be able to use these signals to support prediction and brain-computer
interfacing. To decode each subject’s behavior from their spatial pattern of traveling waves, we
trained multilayer neural networks, with cross-validation, for classifying pairwise behavioral
states. We performed this decoding separately for the spatial and verbal memory tasks (Methods,
Figure 10, Supplementary Figure 20). We used the extracted activation functions (weights) as
inputs for the neural network classifiers.

We found that behavioral states can be decoded at the individual cluster level and also at the
group level (Figure 10C, p < 0.001, one-sided sign test), for the spatial memory task.
Interestingly, we found that some behavioral states can be more easily decoded than others
(Figures 10C, F). For example, Figure 10C shows that for cluster 7, the spatial wave patterns for
the encoding and distractor periods were distinct from each other and hence were more easily
distinguishable compared to those from the confidence and feedback periods. Overall, the
navigation and distractor states were generally decodable from other states, but not from each
other, which indicates that these behaviors generally exhibit similar traveling wave patterns
(Figures 10F, G).

Similarly, we found that the identity of the specific memory item that a subject was viewing can
be decoded from the shape of their traveling waves in the verbal memory task, and also at the
group level (Supplementary Figure 20A, p < 0.001, one-sided sign tests across clusters versus
chance). Overall, we often found that letter “H” was the most decodable, followed by the letters
“J” and “Q”, compared to the other letters (Supplementary Figures 20D, E). These results mean
that letters “H”, “J”, and “Q” often elicited distinctive patterns of traveling waves (Jacobs &
Kahana, 2009). Visually, behavioral states that showed lower decoding accuracy had more
similar wave patterns compared to the ones that had higher decoding accuracy (Figure 10C,
Supplementary Figure 20A), thus the results of this network decoding analysis were closely
aligned with the results from the MANOVA analysis (Figures 5, 6). This similarity between the
results of these analyses was robust, as we found mostly positive correlations between the neural
network decoding accuracy and the dissimilarity of traveling wave patterns (calculated as the
mean Euclidean distance between the ICA weights for the given pair of behavioral states)
(Figures 10D, E, Supplementary Figures 20B, C).

Methods
Human subjects
We examined direct brain recordings from 24 patients with pharmaco-resistant epilepsy who
underwent surgery for removal of their seizure onset zones. All patients consented to having
their brain recordings used for research purposes and all research was approved by Institutional
Review Boards. The patients who performed a treasure hunt (TH) spatial memory task (N=9, 4
females, minimum age = 20, maximum age = 57, mean age = 36.6, see below for details) were
part of a larger data collection initiative that can be downloaded from the UPENN-RAM
consortium (URL: http://memory.psych.upenn.edu/RAM), shared by Kahana and colleagues
(Jacobs et al., 2016). Recordings from the patients who performed the Sternberg (ST) verbal
working memory task (N=15, 7 females, minimum age = 20, maximum age = 58, mean age =
35.9, see below for details), were collected at four hospitals (Thomas Jefferson University
Hospital, Philadelphia; University of Pennsylvania Hospital Philadelphia; Children’s Hospital of
Philadelphia, and University Hospital Freiburg). The direct recordings of these patients can be
downloaded from https://memory.psych.upenn.edu/Data (Jacobs & Kahana, 2009).

Electrophysiological recordings and preprocessing
Patients were implanted with different configuration of electrodes based on their clinical needs,
which included both electrocorticographic (ECoG) surface grid and strips as well as depth
electrodes. In this work, we only examined the ECoG grid electrodes on the cortical surface.
ECoG recordings were obtained using subdural grids (contacts placed 10 mm apart). Anatomical
localization of electrode placement was accomplished by co-registering the postoperative
computed CTs with the postoperative MRIs using FSL (FMRIB (Functional MRI of the Brain)
Software Library), BET (Brain Extraction Tool), and FLIRT (FMRIB Linear Image Registration
Tool) software packages. Preoperative MRIs were used when postoperative MRIs were not
available. From these images, we identified the location of each recording contact on the CT
images and computed the electrode location in standardized Talairach coordinates.

Original sampling rates of ECoG signals in the TH task were 500 Hz, 1000 Hz, and 1600 Hz.
ECoG signals in the TH task were downsampled to 500 Hz, if the original sampling rate was
higher, for all subsequent analysis. Original sampling rates of ECoG signals in the ST task were
400 Hz, 512 Hz, and 1000 Hz. Therefore, ECoG signals in the ST task were downsampled to 400
Hz, if the original sampling rate was higher, for all subsequent analysis. We used common
average referencing (ECoG electrodes re-referenced to the average signal of all electrodes in the
grid), similar to our previous studies on traveling waves (Das et al., 2022; Zhang et al., 2018).
Line noise (60 Hz) and its harmonics were removed from the ECoG signals. For filtering, we
used a fourth order two-way zero phase lag Butterworth filter throughout the analysis.

Cognitive tasks
(a) Treasure Hunt spatial episodic memory task
The patients performed multiple trials of a spatial memory task in a virtual reality environment
(Miller et al., 2018; Tsitsiklis et al., 2020) on a laptop computer at the bedside. They controlled
their translational and rotational movements through the virtual environment with a handheld
joystick. In each task trial, subjects explored a rectangular arena on a virtual 3D beach to reach
treasure chests that revealed hidden objects, with the goal of encoding the location of each item
encountered (Figure 1A). The virtual beach (100 × 70 virtual units, 1.42: 1 aspect ratio) was
bounded by a wooden fence on each side. The locations of the objects changed over the trials,
but the shape, size and appearance of the environment remained constant throughout the
sessions. The task environment was constructed so that the subject would perceive one virtual
unit as corresponding to approximately 1 foot in the real world. Subjects viewed the environment
from the perspective of cycling through the environment and the elevation of their perspective
was 5.6 virtual units. Each end of the environment had unique visual cues to help the subjects
navigate.

Each trial of the task begins with the subject being placed on the ground at a randomly selected
end of the environment. The subject initiates the trial with a button press, then navigates to a
chest using a joystick. Upon arrival at the chest, the chest opens and either reveals an object,
which the subject should try to remember, or is empty. The subject remains facing the open chest
for 1.5 sec (encoding period) and then the object and chest disappear, which indicates that the
subject should navigate (navigation period) to the next chest that has now appeared in the arena.
Each trial consists of four chests; two or three (randomly selected, so that subjects could not
predict whether the current target chest contained an object, which served to remove effects of
expectation and to encourage subjects to always attend to their current location as they
approached a chest) of the chests contain an object, and the others are empty. Each session
consists of 40 trials, and in each session, subjects visit a total of 100 full chests and 60 empty
chests. The chests are located pseudorandomly throughout the interior of the environment. After
reaching all four chests of a trial, subjects are transported automatically so that they have an
elevated 3/4 overhead perspective view of the environment at a randomly selected end of the
environment. They then perform a distractor task (distractor period), a computerized version of
the “shell game”, before entering the retrieval period. During recall, subjects are cued with each
of the objects from the trial in a random sequence and asked to recall the location of the object.
In each recall period, they first indicate their confidence (confidence period) to remember the
location of the object (“high”, “medium”, or “low”). Next, they indicate the precise location of
the object by placing a cross-hair at the location in the environment that corresponds to the
location of the cued item. After the location of each object of the trial is indicated, the feedback
stage (feedback period) of each trial begins. Subjects are shown their response for each object
and given feedback on their response accuracy, following a point system where they receive
greater rewards for accurate responses. A response is considered correct if it is within 13 virtual
units of the true object location. Mean accuracy across subjects was ~41%.

We analyzed the 1.5-sec long trials from the encoding periods of the task. For the navigation
periods, we analyzed 1.5 sec long segments approximately corresponding to the middle of the
navigation trial. Similarly, for the distractor periods, we analyzed 1.5 sec long time segments
approximately corresponding to the middle of the distractor trial. For the confidence periods, we
analyzed 1.5 sec recording immediately following the presentation of the visual cues. For the
retrieval periods, we analyzed 1.5 sec immediately prior to the retrieval of the objects. For the
feedback periods, we analyzed 1.5 sec time segments immediately following the feedback.

The 9 subjects in the spatial task completed 16 sessions in total. On an average, each subject
completed ~ 128 encoding trials, ~ 80 confidence trials, ~ 32 navigation trials, ~ 80 retrieval
trials, ~ 32 distractor trials, and ~ 75 feedback trials, in each session of the spatial task. Such a
relatively large number of trials for the different task periods ensured a high degree of statistical
power to yield robust results.

(b) Sternberg verbal working memory task
Patients performed multiple trials of a Sternberg verbal working memory task (Jacobs & Kahana,
2009; Sternberg, 1966). In each trial of the task, first the patients viewed a fixation cross and
then they were sequentially presented with one to three English letters on the screen of a bedside
laptop computer (Figure 1B). Each letter appeared on screen for 1 sec. Patients were instructed to
closely attend to each stimulus presentation and to silently hold the identity of each item in
memory. After the presentation of each list, the subject viewed a probe item and was asked to
respond by pressing a key to indicate whether the probe was present in the just-presented list or
whether it was absent. After the key press, the computer indicated whether the response was
correct, and then a new list was presented. Across all letter presentations, patients viewed an 8 or
16 letter subset. Mean accuracy across subjects was ~92%. We analyzed the 1 sec long trials
from the encoding period of each letter.

The 15 subjects in the verbal task completed 32 sessions in total. Each subject completed ~ 72
trials of each letter in a session.

Identification of oscillations
The overall goal of our analysis framework was to identify neural oscillations that were present
across the cortex, and to measure their spatial patterns of propagation across the cortex. To this
end, our first goal was to identify the oscillations that were present on multiple electrodes, and
then (as described in subsequent sections), to map their spatial pattern of propagation as traveling
waves. We focused our analysis on oscillations in the theta/alpha and beta bands. To identify the
oscillations present on multiple channels, we adopted methods similar to our previous approach
(Das et al., 2022; Zhang et al., 2018). This subject-specific algorithm accounts for individual
differences in electrode positions and oscillation frequencies across individuals. To identify
traveling waves, a single oscillation frequency per electrode group is crucial since, by definition,
a traveling wave involves a single frequency and whose phase progressively propagates through
a contiguous region of cortex.

In this algorithm, to identify the frequency of interest for each electrode we first used Morlet
wavelets to compute the power of the neural oscillations at 200 frequencies logarithmically
spaced from 3 to 40 Hz. To identify narrowband oscillations, we fit a line to each patient’s mean
power spectrum (mean across trials) in log–log coordinates using robust linear regression (Das et
al., 2022) (Figure 1D). We then subtracted the actual power spectrum from the regression line.
This normalized power spectrum removes the 1/f background signal and emphasizes narrowband
oscillations as positive deflections (Figure 1E). We identified narrowband peaks in the
normalized power spectrum as any local maximum greater than one standard deviation above the
mean (Figure 1E). We repeated this procedure for each phase of both tasks. The mean of the
peak frequencies of the electrodes of a ECoG grid was defined as the grid’s cluster frequency
(CF). The mean peak frequency across all clusters was ~12 Hz. Since by definition, a traveling
wave involves a single oscillation frequency, we only included, in these analyses, oscillation
clusters where at least 2/3rd of the electrodes in a grid had a narrowband peak in the power
spectrum. We excluded 5 patients in the TH task and 2 patients in the ST task from further
analyses since these patients did not meet this criterion. This resulted in 9 patients in the TH task
and 15 patients in the ST task to be included in all subsequent analyses below. Overall, ~87% of
ECoG electrodes in the TH task and ~86% of ECoG electrodes in the ST task had a narrowband
oscillation in one of the theta/alpha (< 12 Hz) or beta (> 12 Hz) frequency bands, indicating that
an overwhelming number of ECoG electrodes show oscillatory activity.

Note that in many cases, electrodes in the same ECoG grid showed multiple peaks in the power
spectrum, i.e, they had peaks at both the theta/alpha and beta frequency bands. We estimated
traveling waves for each of these narrowband frequencies separately. Moreover, since all
electrodes in an oscillation cluster will not have identical oscillation frequency, we also allowed
for electrodes to show oscillations at nearby but nonidentical frequencies (within 5 Hz of cluster
frequency). All our analysis such as traveling waves and independent component analysis below
were carried out based on the individual frequency of each electrode to account for the variation
of frequencies across electrodes.

Identification of traveling waves
Having identified the electrode clusters that showed oscillations, we next identified and
characterized traveling waves. A traveling wave can be described as an oscillation that moves
progressively across a region of cortex. Quantitatively, a traveling (phase) wave can be measured
as a set of simultaneously recorded neuronal oscillations at very similar frequencies whose
instantaneous phases vary systematically with the locations of the recording electrodes. To
identify these traveling waves we used a localized circular-linear regression approach, assuming
that the relative phases of the oscillation clusters exhibit a linear relationship with electrode
locations locally (Das et al., 2022; Das, Wu, et al., 2025; Das et al., 2023; Das, Zhang, et al.,
2025). This locally circular-linear fitting of phase-location can detect complex patterns
(Ermentrout & Kleinfeld, 2001; Muller et al., 2016) of traveling waves in an oscillation cluster in
addition to planar traveling waves.

To identify traveling waves from the phases of each oscillation cluster, we first measured the
instantaneous phases of the signals from each electrode of a given cluster by applying a 4th-order
Butterworth filter at the cluster’s oscillation frequency (bandwidth [f ×.85, f / .85] where f is
p p p
the peak frequency). We used the Hilbert transform on each electrode’s filtered signal to extract
the instantaneous phase.

Next, we used circular statistics to measure the spatial propagation of phase between neighboring
electrodes in each cluster to identify traveling waves at each time point (Fisher, 1993). To
simplify visualizing and interpreting the data, we first projected the 3-D Talairach coordinates
for each cluster to a 2-D plane using principal component analysis (PCA).

To identify traveling waves, we used a series of two-dimensional (2-D) localized circular–linear
regression that were fit to model the direction of wave propagation in a local subcluster in the
25-mm neighborhood of each electrode. The local regression determines the direction of local
wave propagation in the subcluster surrounding each electrode, by measuring whether the local
phase pattern varies linearly with the electrodes’ coordinates in 2-D. Thus, this regression fits a
single direction of wave propagation locally around each electrode at each moment, and these
measures are then subsequently combined across subclusters to identify larger-scale spatial
patterns of propagation.

In this regression, for each nearby electrode, we first identified the neighboring electrodes that
were located nearby, constituting a sub-cluster of the given cluster. Here, let x and y represent
i i
the 2-D coordinates and θ the instantaneous phase of the ith electrode in a sub-cluster. We used a
i
2-D circular-linear model
𝜃̂ = (𝑎𝑥 + 𝑏 𝑦 + 𝜗)mod 3600,
𝑖 𝑖 𝑖
where 𝜃̂ is the predicted phase, a and b are the phase slopes corresponding to the rate of phase
𝑖
change (or spatial frequencies) in each dimension, and ϑ is the phase offset. We converted this
model to polar coordinates to simplify fitting. We define α = atan2(b, a) which denotes the angle
of wave propagation and 𝜉 = √𝑎2 + 𝑏2 which denotes the spatial frequency. Circular–linear
models do not have an analytical solution and must be fitted iteratively (Fisher, 1993). We fitted
α and 𝜉 to the distribution of oscillation phases at each time point by conducting a grid search
over 𝛼 ∈ [00, 3600] and 𝜉 ∈ [0, 18] . Note that 𝜉 = 18 corresponds to the spatial Nyquist
frequency of 180/mm corresponding to the spacing between neighboring electrodes of 10 mm.

In order to keep the computational complexity tractable, we used a multi-resolution grid search.
We first carried out a grid search in increments of 50 and 10/mm for α and 𝜉, respectively. The
model parameters (a= ξcos(α) and b= ξsin(α)) for each time point are fitted to most closely
match the phase observed at each electrode in the sub-cluster. After having relatively coarse
estimates of α and 𝜉, we then carried out another grid search in increments of 0.050 and
0.050/mm around a  2.50 and  0.50/mm neighborhood of the coarse estimates of α and 𝜉,
respectively, to have refined estimates of α and 𝜉. We computed the goodness of fit as the mean
vector length 𝑟̅ of the residuals between the predicted (𝜃̂ ) and actual (θ) phases (Fisher, 1993),
𝑖 i
2 2
𝑟̅ = √[ 1 ∑𝑛 cos(𝜃 −𝜃̂ )] + [ 1 ∑𝑛 sin(𝜃 −𝜃̂ )] ,
𝑛 𝑖=1 𝑖 𝑖 𝑛 𝑖=1 𝑖 𝑖
where n is the number of electrodes in the sub-cluster. The selected values of α and ξ are chosen
to maximize 𝑟̅. This procedure is repeated for each sub-cluster of a given oscillation cluster. To
measure the statistical reliability of each fitted traveling wave, we examined the phase variance
that was explained by the best fitting model. To do this, we computed the circular correlation 𝜌
𝑐𝑐
between the predicted (𝜃̂ ) and actual (θ) phases at each electrode:
𝑖 i
𝜌 = ∑𝑛 sin(𝜃𝑖−𝜃̅)sin(𝜃̂ 𝑖−𝜃̂̅ )
𝑐𝑐 √∑𝑛 𝑠𝑖𝑛2(𝜃𝑖−𝜃̅)∑𝑛 𝑠𝑖𝑛2(𝜃̂ 𝑖−𝜃 ̅̂)
𝑖=1 𝑖=1 𝑖=1
where bar denotes averaging across electrodes. We refer to 𝜌2 as the wave strength (Das et al.,
𝑐𝑐
2022) as it quantifies the strength of the traveling wave (note that 𝜌2 has been referred to as the
𝑐𝑐
phase gradient directionality (PGD) in some prior studies (Muller et al., 2016; Rubino et al.,
2006; Zhang et al., 2018)). Other features of traveling waves such as the wavelength (2π/spatial
frequency) and the speed (wavelength × frequency) can be readily derived from the parameters
of the above 2-D model. Note that traveling waves in some prior studies were detected and
analyzed by calculating the spatial gradient of the phases of the recordings from ECoG
electrodes (Halgren et al., 2019; Muller et al., 2016), however, phase gradients can only be
calculated in two directions (forward and backward), so only a subset of neighboring electrodes
of a given electrode are included in these analyses of spatial gradient. Since our approach
directly includes all possible neighboring electrodes (termed as a sub-cluster in our analysis)
in the circular-linear regression model, thus calculating phase gradients in all possible directions,
it results in a more efficient estimate of the traveling waves parameters.

We note that because a few of the ECoG electrodes did not have a narrowband oscillation on
certain trials, to allow us to estimate the overall shape of the traveling wave across the cluster,
we estimated the traveling waves for those electrodes by an extrapolation procedure where, for
the given electrode, we substituted the mean of the traveling waves of all electrodes within a 25
mm radius (subcluster) of the electrode under consideration. This extrapolation step was
necessary for classification of each oscillation cluster as one of the wave categories using curl
and divergence analysis (rotational or expanding/contracting or complex, see section below on
Identification of planar, rotational, and concentric traveling waves). This is reasonable since as
mentioned above, an overwhelming number of ECoG electrodes showed oscillatory activity
(~87% of ECoG electrodes in the TH task and ~86% of ECoG electrodes in the ST task had a
narrowband oscillation). Nevertheless, we reran our complete series of analyses without the
extrapolation step, and found that the same modes and clusters were still statistically significant
in both the TH and ST tasks, and also the decoding accuracy was very similar to that without the
extrapolation step, indicating that this step did not have a large effect on our decoding results.

Identification of stable epochs
Since a traveling wave is composed of phase patterns that vary relatively smoothly across space
and time, we sought to characterize the spatiotemporal stability of the traveling waves that we
detected in our localized circular-linear regression approach above. We defined stability as the
negative of the mean (across all electrodes in an ECoG grid) of the absolute values of the
difference between the direction and strength of traveling waves observed at consecutive time-
points, defined as
𝑠(𝑡) = −1 ∑𝑛 |𝜌 (𝑡+1)𝑒𝑗𝛼𝑖 (𝑡+1) − 𝜌 (𝑡)𝑒𝑗𝛼𝑖(𝑡)|,
𝑖=1 𝑖 𝑖
𝑛
where, 𝜌 denotes the strength of the traveling wave for the ith electrode, 𝛼 denotes the direction
𝑖 𝑖
of the traveling wave for the ith electrode, 𝑛 denotes the number of electrodes in an ECoG grid, 𝑠
denotes stability, 𝑡 denotes time, and 𝑗 denotes square root of minus one. We repeated this
procedure for each pair of consecutive time-points and for each trial of the encoding, retrieval,
navigation, etc. periods in the TH task and each trial of the English letters in the ST task.
Stability for each trial was z-scored. We identified stable epochs as those for which all stability
values were above the mean of the stability values for a given trial. We ran our stability analysis
across all trials and detected stable epochs from data recorded from all trials of both tasks
(Figures 1F, G).

Previous research in rabbit field potential recordings (Freeman & Rogers, 2002; Freeman &
Schneider, 1982) have found spatial patterns of theta traveling waves last ~80-100 msec in
duration. Moreover, large-scale, whole-brain computational modeling in humans using neural
field theory have shown that spatiotemporally stable traveling waves last ~50-60 msec in
duration (Roberts et al., 2019). Therefore, we focused our analyses on the stable epochs that
were more than 25 time-points long for further analysis; as there was evidence that shorter wave
patterns may be less stable (Roberts et al., 2019). This corresponded to 50 msec in the TH task
(fs = 500 Hz) and 62.5 msec in the ST task (fs = 400 Hz). Overall, ~48% of stable-epochs were
long enough to be included based on this criteria, yielding ~4 stable-epochs per trial. We also
examined how temporal modulation of these stable epochs varied across task periods
(Supplementary Figure 3).

Identification of modes using complex independent component analysis (CICA)
We next used a large-scale dimensionality reduction framework to measure the spatial patterns
of wave propagation across each cluster. This framework identified the different types of spatial
wave patterns that formed traveling waves at each cluster and modeled how they summed to
determine the signal on each epoch. We thus applied complex ICA (CICA) to the data on each
epoch. Within each stable epoch identified as described above, the direction and strength of the
traveling waves remain almost consistent across time-points. Thus, we averaged the direction
and strength of the traveling waves within each stable epoch to find one spatial wave pattern
associated with each stable epoch, which was input to the CICA algorithm.

Next, to analyze how the spatial patterns of waves changed over time and related to behavior, we
concatenated the wave patterns (i.e., direction and strength) for all stable epochs across
encoding, retrieval, navigation, etc. periods into a single matrix and then passed this matrix as
input to the complex version of the independent component analysis (CICA) (Fu et al., 2015; Li
& Adalı, 2010). We used this complex version of the ICA, as compared to real ICA, to
incorporate the 2-D directions of the traveling waves, weighted by the strength (𝜌cos (𝛼) and
𝜌sin (𝛼)), defined for each electrode in the ECoG grid. We then extracted the independent
activation functions (or, weights) (each activation function corresponds to one of the stable
epochs) and the corresponding modes (“raw modes” in Figure 2) as the output from the CICA
(Fu et al., 2015; Li & Adalı, 2010). Multiplication of each of these modes with the mean of the
weights across all stable epochs corresponding to that specific mode results in a unique wave
pattern associated with that mode (the “mean mode” in Figure 2). At the individual epoch level, a
higher CICA weight for that epoch corresponding to a specific mode indicates higher
representation of that wave pattern in that specific epoch and a lower CICA weight for an epoch
corresponding to a specific mode indicates lower representation of that wave pattern in that
specific epoch (Figure 2). Moreover, the higher the variance explained by a given mode, the
higher will be its representation across the trials (see Figure 4E). Generally, we observed more
complex wave patterns in higher order modes (modes beyond 6), which suggests that complex
wave patterns explain relatively less variance in the data compared to spirals and concentric
wave types. In this way, we can extract the CICA weights for each of the encoding, retrieval,
navigation, etc. periods in the TH task and each of the English letters in the ST task and directly
compare them statistically; see Statistical analysis section for details.

Identification of planar, rotational, and concentric traveling waves
After extracting the mean modes from the CICA procedure above, we next sought to classify
each of the mean modes into one of the following categories “planar”, “rotational”, “concentric”
(“expanding” or “contracting”), or “complex” (Figure 3, Supplementary Figures 9-18)
(Townsend et al., 2015), to identify global patterns of traveling waves in particular shapes that
might be meaningful for cortical processing. We will be using the terms “modes” and “wave
patterns”/“spatial patterns” interchangeably throughout our manuscript. This is because each of
the modes corresponded to one unique and independent spatial pattern extracted from the CICA
algorithm.

For classifying each of the modes into a specific wave pattern, we carried out robust statistical
tests. We constructed appropriate surrogate distributions by shuffling the coordinates of the
electrodes which destroys the specific structure associated with a given wave pattern. Please note
that the raw modes are not unique since the activation functions and raw modes can be scaled to
an arbitrary complex factor to yield the same output. Therefore, all our statistical procedures are
done on the mean modes for accurate results.

For classifying the specific shape of each wave, we created a series of scores, which identified
plane waves, and rotational/concentric waves. Clusters that did not meet our criterion for
planar/rotational/concentric waves were labeled as ‘complex’. For detecting planar traveling
waves, we calculated the mean wave direction (weighted by the strength of the wave) of all
electrodes in an oscillation cluster. We then constructed surrogate distributions by shuffling the
coordinates of the electrodes and recalculating the mean wave direction for each shuffling,
therefore constructing a null distribution against which the empirical data can be compared.
From the surrogate distribution, the threshold that corresponded to significance at p = 0.001 was
~ 0.6 across all clusters. Thus, we declared a wave as a planar wave for which the mean wave
direction was higher than 0.6. Visually, this procedure of threshold selection identified individual
mean modes that showed clear planar waves (Figure 3) and small changes in this threshold did
not substantially change the results.

Similarly, we classified the mean modes that represented rotational and concentric traveling
waves by measuring the shape of their wave propagation using the curl and divergence metrics
respectively (Muller et al., 2016). Curl measures rotational patterns (for example, clockwise or
counter-clockwise rotation) in wave dynamics and divergence measures expanding/contracting
patterns (for example, sources or sinks) in wave dynamics (Muller et al., 2016). We first
calculated the curl and divergence metrics for each electrode and then calculated the mean curl
and divergence across all electrodes in an oscillation cluster.

The curl of a traveling wave at a particular electrode can be calculated as follows. For a 2D
traveling wave 𝐹(𝑥,𝑦) = 𝐹 (𝑥,𝑦)𝑒̂ + 𝐹 (𝑥,𝑦)𝑒̂ at electrode locations (𝑥,𝑦), where 𝐹 (𝑥,𝑦)
𝑥 𝑥 𝑦 𝑦 𝑥
and 𝐹 (𝑥,𝑦) denote the traveling wave components in x and y directions respectively and 𝑒̂ and
𝑦 𝑥
𝑒̂ denote the unit vectors in the x and y directions respectively, curl can be defined as
𝑦
𝑐𝑢𝑟𝑙 𝐹 = ∇×𝐹 = (𝜕𝐹𝑦− 𝜕𝐹𝑥)𝑒̂
𝑧
𝜕𝑥 𝜕𝑦
Similarly, the divergence of a traveling wave at a particular electrode can be calculated as
follows
𝑑𝑖𝑣 𝐹 = ∇.𝐹 = (𝜕𝐹𝑥+ 𝜕𝐹𝑦)𝑒̂
𝑧
𝜕𝑥 𝜕𝑦
We used the same approach to assess significance of these labels as above, constructing
appropriate surrogate distributions by shuffling the coordinates of the electrodes and then
recalculating the curl and divergence. For p = 0.001, we identified curl and divergence thresholds
of ~ 0.75 and ~ 0.4 for the rotational and concentric waves, respectively, and these values were
very similar across clusters. We thus declared a wave as a rotational or concentric wave for
which the curl or divergence value was higher than 0.75 or 0.4 respectively. Visually, this
procedure of threshold selection functioned as we expected, as the waves detected with these
thresholds showed clear rotational and concentric dynamics visually (Figure 3) and small
changes in this threshold did not substantially change the results reported here.

Modes that were not classified as one of the planar, rotational, or concentric, were designated as
complex waves. Complex waves were also stable at the individual epoch level as shown in
Figure 9 and also across epochs, as these waves were able to distinguish behavior, for example,
see Supplementary Figure 18. Even though there was no global pattern associated with these
complex traveling waves, many of these complex waves showed interesting local patterns, where
a subset of electrodes in an ECoG grid showed planar, rotational, or concentric waves (Figure 9,
Supplementary Figures 9-18).

If the size of the ECoG grid is small, such as the 64 or 84 grids in our dataset, there is a
possibility that the detected planar waves could represent segments of larger rotational waves.
Therefore, the number of rotational waves that we have detected in our dataset could potentially
be undercounted. Future studies with much larger ECoG arrays across many individuals are
needed for a more accurate assessment of the number of rotational waves and their functional
role during various memory tasks.

Validating CICA using simulations
In addition to statistical validation on the experimental data using shuffling, we have also carried
out four types of simulations to validate our proposed CICA approach to confirm its sensitivity
to extract expected wave patterns and demonstrate specificity to reject negative results in the
absence of such patterns.

Simulation 1, Sensitivity test, single wave patterns: In this simulation, we test the sensitivity
of our framework where we feed-in known traveling wave patterns and then test the robustness
of the CICA framework against noise. We first simulated a plane wave using an 8*8 electrode
grid and added Gaussian random noise to generate multiple trials (100 in total) and then passed
these into our CICA algorithm to test how well our algorithm is able to estimate the wave
pattern. We generated activation functions or weights from a uniform distribution. To quantify,
we calculated the root-mean-squared-error (RMSE) between the ground-truth pattern and the
estimated pattern by the CICA approach and we repeated this procedure from low to high signal-
to-noise ratio (SNR). RMSE for this simulation and all subsequent simulations below was
calculated by taking the root of the sum of squared differences of the activation functions
estimated by CICA and the activation functions in the original data that we simulated. This
demonstrates the ability of CICA to not only reconstruct the modes present in the original data,
but also the degree of representativeness of the individual modes timepoint-by-timepoint. If the
activation functions estimated by CICA match the activation functions in the original data
timepoint-by-timepoint, then it can be claimed that CICA has high sensitivity. Additionally, we
also tested the robustness of our CICA approach using rotational wave, source wave, and
complex wave patterns.

Results from this analysis are summarized in Supplementary Figure 4A, which showed that the
performance of the CICA approach improved as SNR increased, with near-perfect reconstruction
of the underlying traveling wave patterns at higher SNR values. Moreover, the performance of
the different wave shapes such as planar, source, rotational, and complex waves were very
similar, demonstrating the ability of the CICA approach to robustly estimate a diverse range of
wave patterns. Simulation 1 demonstrates the sensitivity of the CICA method for a single wave-
type and CICA is able to reconstruct the waveforms for each individual timepoints, in addition to
reconstructing the shape of the waveform, i.e., the mode.

Simulation 2, Sensitivity test, multiple wave patterns: In this simulation, we test the
sensitivity of our framework for identifying and distinguishing combinations of traveling wave
patterns, i.e., multiple modes, in the presence of noise. For this, we simulated planar+source,
planar+rotational, planar+complex, and planar+source+rotational wave patterns. We generated
the activation functions for each type of wave shape from a uniform distribution. We again added
Gaussian random noise to generate multiple trials and then passed these into our CICA algorithm
to test how well our algorithm is able to separate multiple wave patterns. To quantify, we
calculated the root-mean-squared-error between the ground-truth pattern and the estimated
pattern by the CICA approach for all wave shapes and we repeated this procedure from low to
high SNR.

Results from this analysis are summarized in Supplementary Figure 4B, which showed that the
performance of the CICA approach improved as SNR increased, similar to results from
Simulation 1, with near-perfect reconstruction of the underlying traveling wave patterns at higher
SNR values, demonstrating the ability of the CICA approach to robustly estimate the underlying
wave patterns when multiple traveling waves are simultaneously present. For example, if the
degree of representativeness of the rotational wave-type is high at a specific time-point, then
CICA is accurately able to reconstruct this waveform at that individual time-point, and therefore
the RMSE is close to zero at reasonably high SNR in our simulations.

Simulation 3, Specificity test: We also tested the specificity of the CICA method in the
presence of spatially autocorrelated signal whose spatial autocorrelation matches to those in the
spatial and verbal experimental data. We considered a rotational wave. Note that these rotational
waves were ubiquitous in our experimental data and the spatial autocorrelation of the rotational
wave matches to those in the experimental data. We also created spatially circularly shifted
versions of the original rotational wave: vertically shifted, horizontally shifted and both
vertically and horizontally shifted versions (Supplementary Figure 5). Note that the spatial
autocorrelation of all these circularly shifted waves is the same as the original rotational wave,
although the specific structure is now destroyed. We separately gave these circularly shifted
waves and also the original rotational wave as inputs to the CICA method. For the circularly
shifted waves, we created 500 permutations of the circular shift and for each permutation, we
generated 50 trials by adding Gaussian random noise to the traveling waves and passing them to
the CICA method as input. This analysis revealed that CICA cannot reconstruct the original
rotational wave anymore after using the circularly shifted versions. The error magnitude when
the original wave is given as input to the CICA method is significantly lower than the
distribution of the error magnitudes when the circularly shifted waves are given as input to CICA
(p < 0.001). This analysis strongly demonstrates the ability of CICA to show negative results in
the absence of any specific shape of wave patterns.

Simulation 4, Joining tests of sensitivity and specificity: This test assesses the joint ability of
our framework to distinguish individual trials that showed real versus no waves, within an
individual dataset. For this, we simulated rotational waves and planar+source+rotational wave
patterns. The performance of other wave patterns such as planar waves and source waves was
similar to the rotational wave pattern and hence omitted. For the input, half of the trials consisted
of the given wave pattern (rotational or planar+source+rotational) and the other half consisted of
no wave (zero signal, corresponds to zero weights or activation functions). We then calculated
the detection percentage of our CICA approach, i.e., we calculated the percentage of trials our
CICA approach is able to correctly identify the presence or absence of the wave patterns and we
repeated this procedure from low to high SNR.

This analysis revealed that the CICA correctly estimates both the presence and absence of wave
patterns (Supplementary Figure 4C). For the rotational wave shape, CICA is able to achieve
near-perfect detection accuracy at > ~ 15 dB SNR and for the planar+source+rotational case,
CICA is able to achieve near-perfect detection accuracy at > ~ 20 dB SNR, for both the
sensitivity and specificity tests, demonstrating the ability of the CICA to extract the expected
wave patterns and also demonstrate negative results in the absence of such patterns. Simulation 4
demonstrates both the sensitivity and specificity of the CICA method. Simulation 4 demonstrates
that the time-points where there was a wave present, CICA is able to accurately reconstruct both
the shape and the degree of representativeness of that waveform at that individual time-point. On
the other hand, CICA is also able to estimate the value of zero activation function with high
accuracy at a given time-point where there was no wave, indicating that CICA can produce null
result when there are no waves. Furthermore, in the case where we simulated multiple
waveforms in equal contribution such as the “planar+source+rotational” case, the CICA method
was able to reconstruct all of these three waveforms in equal contribution. Therefore, Simulation
4 strongly demonstrates both the sensitivity and specificity of the CICA method.

We would like to clarify that the individual time-points in our simulations correspond to the
individual stable epochs in the experimental data. The shape of the traveling waves changes from
one stable epoch to the next stable epoch, therefore the wave shapes also change from one task
(for example, encoding) to the other (for example, retrieval). The CICA method is able to capture
the variability of the changes in the wave shapes by using the activation functions whose value
dynamically changes across different stable epochs and task periods, indicating the sensitivity
and specificity of the CICA method. Similarly, in the simulations, the CICA method is able to
capture the presence or absence of waveforms or the degree of representativeness of the
waveforms by estimating the dynamic or time-varying activation functions. Since, in our
simulations, we simulated planar, spiral, and concentric waves, which were ubiquitous in our
experimental data and whose spatial autocorrelation is highly similar to the experimental data,
our simulations therefore very strongly support the claims in our experimental data as well.

Statistical analysis of the relation between traveling waves and behavior
(i) Identifying significant modes in CICA
To estimate the number of significant modes in the CICA output for each cluster, we compared
the variance explained for each mode with the theoretical variance threshold 100/n, where n is
the number of electrodes in an ECoG grid. This theoretical variance corresponds to the variance
of each mode if the total variance (100%) is equally distributed among all modes. Because of the
spatial structure of traveling waves, some modes will explain more variance compared to the
other modes in the empirical data (Figure 4E). We additionally shuffled the electrodes in each
cluster and recalculated the variance distribution across modes and confirmed that the shuffled
variance for all modes converged to the theoretical variance threshold of 100/n.

(ii) Multivariate analysis of variance (MANOVA)
To identify how traveling wave patterns related to behavior, we directly compared the weights
estimated from the CICA procedure above between encoding, retrieval, navigation, etc. periods
in the TH task and between the English letters in the ST task, for each mode using multivariate
analysis of variance (MANOVA). MANOVA was used because the weights were complex. We
statistically distinguished weights corresponding to different behavioral states using the
following model: Real + Imag ~ States, where Real and Imag are the real and imaginary parts of
the weights respectively and States are encoding, retrieval, navigation, etc. in the TH task or
letters “B”, “G”, etc. in the ST task. We used this model for each mode and applied FDR-
corrections for multiple comparisons (p < 0.05) across all modes and oscillation clusters.
Statistical significance from the MANOVA would indicate that traveling waves shift their
direction and/or strength to form distinct directional patterns which can distinguish different
behavioral states. We designated an oscillation cluster to be significant if at least one of the
modes from the CICA for that cluster showed statistical significance in MANOVA. Moreover,
we also repeated our MANOVA analysis by shuffling the labels of the different behavioral states
such as encoding, retrieval, etc. and built a histogram from these shuffles (f-statistics) against
which we compared the empirical data. This analysis revealed that all our MANOVA results
were still statistically significant (p’s < 0.05).

(iii) Robustness of traveling waves at the individual epoch level
We conducted surrogate analysis to test the significance of the estimated stable epochs (see
Identification of stable epochs section above) and whether the observed stable epochs are beyond
chance levels. For each time-point, we shuffled the trial labels so that the temporal contiguity
within a given trial is destroyed and also the electrodes, so that the spatial topography for the
corresponding behavioral state is destroyed, and then ran the stable epoch analysis using
identical methodology as above. In this way, we built a surrogate distribution by aggregating all
time-points corresponding to these shuffled stable epochs against which we then compared the
aggregated time-points from the empirical stable epochs (p < 0.05).

Additionally, we also carried out statistical tests based on spatial autocorrelation-preserving
surrogates to ensure that the traveling waves corresponding to the different behavioral states are
statistically robust. To do this, we circularly rotated the coordinates of the electrodes for each
stable epoch and passed these shuffled wave patterns into our SCICA approach and subjected the
corresponding results to MANOVA to estimate the f-statistics for these shuffled wave patterns.
Surrogates tests are powerful non-parametric tests that are adaptive to the data in consideration,
and highly robust to the specific distribution of the data, therefore reducing Type I and Type II
statistical errors usually found in parametric tests. F-statistics from empirical data are then
compared to this histogram of f-statistics and then significance is assessed based on whether the
empirical data is outside 95% (p < 0.05), 99% (p < 0.01), or 99.9% (p < 0.001), of the histogram.
Note that this shuffling preserves the spatial autocorrelation of the wave patterns, but destroys
the specific structure associated with them, therefore destroying the systematic relationship
between the wave direction and task state. We then constructed the histogram of the f-statistics
of these surrogates and then compared the empirically estimated f-statistics to this surrogate
distribution. This analysis revealed that our results were highly significant across all clusters (all
p’s < 0.001).

(iv) Assessing group-level significance
For assessing statistical significance for the group level results related to fraction of wave
strength and fraction of significance (Figures 4, 7-9, Supplementary Figure 19), we used chi-
squared tests with FDR-corrections for multiple comparisons (p < 0.05) across frequencies
(theta/alpha and beta) and hemispheres (left and right).

(v) Control analysis to assess the effects of oscillation frequency on behavior
To rule out the possibility that the behavioral results could be driven by changes in spectral
content (e.g., oscillation frequency), we ran the CICA analysis on the frequency contents of each
episode. Specifically, we estimated the peak frequency for each trial and electrode separately,
and then passed these measures to the same CICA and MANOVA analysis as in our wave
analysis approach. However, this analysis revealed that the spectral contents do not differ
between any of the task periods in the spatial task or any of the letters in the verbal task (all p’s >
0.05), indicating that the effects that we observed cannot be explained by differences in spectral
contents, rather they are related to the traveling waves.

We also carried out additional analysis by passing the peak frequencies during different tasks to
train our neural network decoder to distinguish pairwise behavioral states. As before, we used
our shuffling analysis procedure by shuffling the labels of the behavioral states such as encoding,
retrieval, etc. and then running the same decoding analysis on this shuffled data to generate a
distribution of shuffled data against which the actual decoding performance can then be
compared. However, we could not decode the behavioral states by using the peak frequencies,
neither in the spatial task nor in the verbal task (all p’s > 0.05). This suggests that the decoding
performance is driven by the different wave patterns that we observed across behaviors, and not
driven by a difference in peak frequencies across different task conditions. These new results
reinforce the fact that the effects that we observed cannot be explained by differences in peak
frequencies, rather they are related to the traveling waves.

(vi) Control analysis to assess the effects of filtering
Since our oscillation detection algorithm is based on identifying peaks in the power spectrum
which is averaged across trials, we carried out a signal-to-noise ratio (SNR) analysis to ensure
that these oscillations are actually present at the individual stable epoch level. Specifically, we
concatenated all the filtered recordings at the oscillation frequency and raw recordings across all
the stable epochs for each cluster and then calculated the Pearson correlation between the filtered
and raw data. This analysis revealed that the correlation values are statistically significant across
all clusters (median Pearson r ~ 0.6, all p’s < 0.001). These high correlation coefficients indicate
that the narrowband filtered signals are capturing a large amount of signal in the raw recording,
showing that the filtering procedure is capturing a dominant amount of variance in the
recordings.

Moreover, we also ran our SNR analysis for the different wave modes such as planar, rotational,
concentric, and complex, to examine whether the SNR differs for the different wave modes. We
first classified each stable epoch as one of the modes (planar/rotational/concentric/complex)
based on the magnitude of the activation function, i.e., a stable epoch was classified as a planar
wave if the magnitude of the activation function corresponding to the planar mode was the
highest, etc. We then concatenated the filtered and raw recordings across all the stable epochs for
each mode and calculated the Pearson correlations between the filtered and raw data for each
mode. We then compared the correlations for the pairwise modes (planar vs. rotational, planar
vs. concentric, rotational vs. concentric, etc.) for each oscillation cluster. This analysis revealed
that the correlation values for the pairwise modes do not differ from each other in any cluster (all
p’s > 0.05, Mann-Whitney U-test). This suggests that the different wave patterns that we
observed in our data are not driven by a difference in SNR.

To rule out the possibility that the effects could be driven by changes in SNR, we calculated the
Pearson correlation between the filtered and raw data for each of the behavioral states separately,
both for the spatial memory and verbal memory tasks. We then used these correlation values in
our neural network decoding analysis to distinguish pairwise behavioral states. As before, we
used our shuffling analysis procedure by shuffling the labels of the behavioral states such as
encoding, retrieval, etc. and then running the same decoding analysis on this shuffled data to
generate a distribution of shuffled data against which the actual decoding performance can then
be compared. However, we could not decode pairwise behavioral states by using these
correlation values (all p’s > 0.05). This suggests that the decoding performance is driven by the
different wave patterns that we observed across behaviors, and not driven by a difference in SNR
across different task conditions.

Decoding analysis of the relation between traveling waves and behavior
Our final goal was to test whether we could decode the broad and specific behavioral states from
the diverse traveling wave patterns in our datasets. Decoding behavioral states in our datasets is a
multiclass classification problem. We converted this problem into a series of binary classification
tasks, as these can be solved straightforwardly with various multivariate algorithms. We trained
multilayer neural networks (Bernardi et al., 2020), with cross-validation, for classifying pairwise
behavioral states (for example, encoding versus retrieval, letter “D” versus letter “J”, etc.),
separately for the spatial and verbal memory tasks, training a neural network classifier for each
oscillation cluster separately. We used a PyTorch-based four-layer neural network for this binary
classification problem (Paszke et al., 2019). We used the extracted weights from the CICA
procedure as features for training our neural network classifiers. We observed that increasing the
number of modes for training increased the network test decoding accuracy (Figure 10B).
Therefore, we included all the modes for classifying the behavioral states. Also please note that
the higher order modes (modes beyond 6) were mostly complex wave shapes and, these higher
order modes progressively contributed less to the test decoding accuracy (Figure 10B), indicating
lack of specific wave structure in these higher order modes compared to the lower order modes
(modes lower than 6).

The multilayer network architecture comprised an input layer, two hidden layers (32 and 16
neurons respectively), and an output layer with a sigmoid activation function, similar to the
neural network architectures previously used for classifying behavioral states (Bernardi et al.,
2020). The inputs for all our decoding analysis were the activation functions or the weights since
for a given mode, these were changing across different task periods. Please note that for a grid
with “N” number of electrodes and hence “N” number of dimensions in CICA, a 2N dimensional
neural network was trained (the factor 2 comes from the fact that the weights are complex
numbers and hence real and imaginary parts were considered separately). ReLU activation
functions were applied to the hidden layers to introduce nonlinearity and improve generalization
(Dahl et al., 2013). The network architectures were kept the same across the spatial and verbal
memory datasets, and individual models were trained for each oscillation cluster separately to
optimize neural network hyperparameters. To further enhance generalization, dropout
regularization was implemented after each hidden layer (Srivastava et al., 2014). For model
optimization, we used Binary Cross Entropy as the loss function, ideal for binary classification
tasks (Ruby & Yendapalli, 2020), and the Adam optimizer during training (Paszke et al., 2019).
To ensure robust optimization and generalization, we employed cross-validation for iterative
training over multiple epochs. For each pair of behavioral states, we first rebalanced the data
such that we have an equal number of epochs per behavioral state (van Gerven et al., 2013). We
used a five-fold cross-validation technique in which the network model was fitted on 80% of the
data and then its performance was tested on the remaining 20%. We used a subsampling
procedure where each epoch was randomly assigned to a particular fold, subject to the constraint
that all behavioral states are evenly represented (van Gerven et al., 2013), and then averaged the
network decoding accuracies across folds.

To access group level statistical significance for the decoding accuracy results, we used one-
sided sign tests versus chance (p < 0.05). Moreover, we also carried out additional statistical
analysis by shuffling the labels of the behavioral states such as encoding, retrieval, etc. and then
running the same decoding analysis on this shuffled data to generate a distribution of shuffled
data against which the actual decoding performance can then be compared. This analysis
revealed that our decoding results are statistically significant across all clusters (all p’s < 0.05),
suggesting that structure of wave activity is the source of the classification accuracy.

Please note that all our traveling waves analysis were carried out for each cluster and subject
separately since each cluster showed specific patterns of wave strength and directions, which
differed across clusters and subjects. Due to the variability of the strength and directions of the
traveling waves, our decoder was also run in each cluster separately rather than collapsing across
subjects (Figure 10, Supplementary Figure 20). Group level results are finally shown in Figure 4
and Supplementary Figure 19 to compare the wave-strength for different types of waves and task
periods.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors relate traveling wave propagation in the human brain (ECoG, epilepsy patients) to behavior (spatial and verbal\nmemory tasks). Traveling waves were defined as propagated oscillations at fixed frequency. Procedure: 1) Find an\noscillation at a single frequency that is present simultaneously on multiple electrodes; 2) Identify a contiguous cluster of\nelectrodes with oscillations at the same frequency; 3) Find a systematic propagation of phase across this cluster of\nelectrodes. Their approach identifies waves that last ~80-180 ms and form complex patterns, such as spirals, sources and\nsinks, besides simpler plane waves. They also used Independent Component Analysis (ICA) to identify and quantify the\nwave patterns into modes and detect how much each mode was present at any given time. Results were validated using\nsimulations.\nResults:\nAverage of 6 significant modes, explaining > 80% of variance (first three explained ~34%, ~19%, ~10% = ~64%). Found\ndifferential presence of patterns across hemispheres during tasks. Concentric waves also differed between spatial and\nverbal tasks. Waves had significant shifts in direction or strength between stages of the spatial memory task. Waves also\noften changed direction or strength for individual letters in the verbal memory task. Behavioral states could be \"reliably\"\ndecoded (~60% accuracy, chance at 50% -- achieved ~70% with ~20 modes) for the spatial memory task. Some states were\nbetter distinguishable by the decoder than others, and this correlates with how distinguishable their wave patterns were.\nLetters were also decoded (~50-70% accuracy, depending on the letters being compared; chance at 50%), and again this\nwas correlated with how distinguishable their wave patterns were.\nCritique:\nThe paper presents a study on traveling waves in humans during task and show that diverse patterns of traveling waves are\nnot only present but also correlate with behavior. The conclusions are backed up by the results and methods look sound.\n1. It is worth noting that a lot of the work is based on previous studies, including some from the same authors. This raises the\nproblem that the novelty of the current study simply arises from applying their established approach to their new datasets.\n2. I would like to see the figure legends be more descriptive. Often, I could not fully understand what was being plotted (e.g.,\nFig. 1i and SFig. 4c).\n3. Decoding numbers are relatively low, but I'm not sure if one would expect much better given the info fed and model used.\nHowever, the authors should soften the claims of \"reliable\" and \"robust\" decoding.\n4. The description of modes and their relation to wave patterns could be clarified more, it feels like modes are essentially\ncapturing pattern types, but modes and pattern types are used on different results back and forth.", "reviewer_2": "This is a technically rigorous, data-rich study that provides a detailed taxonomy of wave patterns observed during memory\ntasks. A comprehensive set of controls was implemented to identify traveling waves and link them to behavioral states.\nWhile the analyses are thorough and address potential confounds in depth, I believe the manuscript would benefit from a\nbroader examination of the functional role of these patterns to better suit the broad readership of Nature Communications.\nLimitations:\nThe reliance on pairwise classification appears methodologically restrictive. Although the authors suggest that their analysis\nof traveling waves may yield novel insights into the brain’s spatial organization, it is not evident how binary discrimination\nalone substantiates this claim. Complex neural datasets frequently exhibit nonlinear decision boundaries and latent\nstructure, such as the hierarchical organization observed in episodic memory tasks. Capturing such complexity may require\ntargeted feature engineering and model architectures that explicitly accommodate multi-class or structured outputs.\nMoreover, performance in a binary framework does not imply that the learned feature representations will generalize to the\nmulti-class setting necessary for comprehensive behavioral decoding, which appears to be the broader aim of this work.\nApproaches employing deeper architectures or ensemble methods may be more appropriate for uncovering meaningful\nstructure in the data. In the absence of such methods, it is difficult to assess the strength of the assertion that “behavioral\nstates of human memory representations can be decoded from traveling waves.” Given that MANOVA has already\ndemonstrated statistically significant differences between states, the added value of the multi-layer network analysis remains\nunclear.\nAlso, a key limitation is the omission of spike rate analyses. Claims that waves “dynamically strengthen” or “transiently\nreorganize” functional connectivity or route spiking should be supported by spike–field coherence or spike rate comparisons\nacross the new wave patterns (e.g. \"complex waves\") considered here.\nI would appreciate a more detailed discussion of the potential mechanisms underlying the predominance of one wave class\nover another. For instance, the observation that source waves are more prevalent in the spatial memory task is intriguing and\ncould potentially reflect changes in inhibitory dynamics and reorganization of the theta cycle between encoding and\nretrieval. Similarly, rotating waves might arise from spatial inhomogeneities, but it is unclear why such factors would\nspecifically influence the spatial memory task more than the verbal task (Fig. 4).\nLast, the potential impact of array positioning on the observed results is unclear. For instance, could the detected plane\nwaves represent segments of larger rotating waves, and if so, how was this possibility addressed? Likewise, could the\nobserved complex wave patterns be artifacts rather than planar or rotating waves?", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.\nThe paper studies how widespread brain regions coordinate their activity to support complex behaviors, which is a\nfundamental question in neuroscience. The authors demonstrate the existence of different types of travelling waves in the\nhuman cortex during memory tasks (planar, spiral, concentric, and complex traveling waves), and show how these specific\nwave patterns correlate with and can be used to decode distinct behavioral states. This is a major contribution.\nThe methodology is clear and robust, especially after the revisions, where various methodological issues were raised.\nMoreover, the authors showed that neither peak oscillation frequency nor SNR could be used to decode behavioral states,\nstrengthening their conclusion that the observed effects are genuinely related to the patterns of traveling waves.\nGiven the novelty of the findings, the methodological rigor, and the comprehensive response to previous critiques, this paper\nmakes a strong case for publication. It is likely to be of broad interest to the neuroscience community and to stimulate further\nresearch into the functional role of complex traveling waves in cognition.\nAs a minor comment, in the the section about the generative mechanisms, the authors should also refer to computational\nmodels that describe the appearance of frequency specific activation patterns (static or travelling) due to axonal time-delays,\nsuch as the normalized connectome by Petkoski & Jirsa 2022 Netw Neurosci, and modelling works for whole-brain activity\nsuch as: Cabral et al. 2023 Nat Comm Biol and Petkoski et al 2019 PTRSA.\nSimilarly, important are the modeling works on the thalamo-cortical circuits which focus on the cortical geometry and were\ndeveloped by the group of Peter Robinson and lately applied by the lab of Michael Breakspear.\nFinally, as a suggestion to both, the authors and the editor, I find it very difficult to read the work in this kind of format, with the\nfigures and captions separated. I strongly suggest including a reader-friendly file with the figures and captions in the text.\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_4": "This paper is devoted to measuring and classifying waves of activity in the human brain performing tasks. The key results\nare that there are many different kinds of such waves, and that these kinds, which can be complex, change in a predictable\nway in the course of carrying out a task. The authors interpret their results to imply that these waves have structure that can\ngive insights into cognition, such as memory.\nMuch of the paper is highly technical and, to this reader, very difficult to get through. The basic kinds of simple patters (plane\nwaves, sources, sinks, spirals inward, spirals outward) are easy to understand. But the data present more complex\nphenomena, and a lot of effort is given to choosing data that can represent “stable epochs” and then using ICA to classify the\nstable epochs in terms of simpler phenomena. I did not understand the details of what they were trying to do (e.g., with\n“activation functions”). It didn’t help that static pictures don’t fully capture movements in a dynamic pattern. I did try to look up\nsome videos with links in the paper; I got a message that they were too big to run on my computer.\nHowever, their big message did come through: the dynamic structure of these patterns contains information that reflect\ncognitive states. That some aspects can reflect broad cognitive states has been argued in previous papers of this group and\nthat of Earl Miller and collaborators. For example, it is argued in previous papers that some plane waves associated with\nlearning go in the opposite direction from waves associated with recall. This makes perfect sense functionally, since the\nregions of the brain participating in learning or recall might be involved in temporally opposing manners in the two functions.\nHowever, this paper makes the further claim that these waves can contain information about very specific information in\nmemory, such as telling one letter of the alphabet from another. They also suggest that other types of waves (non-planar)\nmay also be associated with function. As I explain below, these assertions seem unlikely to me and suggests to me that the\nphenomena associated with such specific items and some kinds of non-planar waves may be epiphenomena: there are\nassociations, but the associations do not reflect cognitive causality. The very interesting result of this paper that the nature of\nthe waves changes over the course of a single trial of an experiment does not, for me, constitute a support for the functional\nuse of all of the waves. However, I should note that the paper does not claim that the waves are important for coding, only for\ndecoding, which may be consistent with epiphenomena and interesting by itself.\nMy intuition about the use of waves such as “targets” come from the chemical experiments in which such traveling waves\nwere found, namely the “Belousov Zhabotinsky” solution, a chemical solution that oscillates between red and blue when it is\nstirred. In a thin layer (and covered to prevent convection) it creates patterns of outwardly going targets that continue to\nchange in time. For me, the essential point here is that the centers of these targets are random in space. If the target patterns\nin the brain are produced by a similar mechanism, it is hard to see how such an unregulated pattern could be of use in any\nhighly regulated function. The spiral patterns in the BZ solution have centers that correspond to some heterogeneity, so\nmaybe it is possible for them to be somewhat regulated. However, I’m finding it very hard to imagine a situation in a wave\n(even a planar wave) corresponds to a letter of the alphabet. Since the latter is learned, what kind of learning rule would set\nup a wave?\nWhat I would have liked to see is a deeper discussion of why the authors think that the new patterns they found, (e.g., the\ntargets and spirals) actually contribute to cognitive processing, as opposed to being associated with processing. One of the\nauthors is a world class expert on such patterns and may have thought more about this than is in the paper. I would also like\nto have seen a deeper discussion of the cited modeling paper by Bhattacharya et al, since the authors of the current paper\nsay that they reproduce the results of that paper experimentally. Do they think their waves come from thalamocortical\ninteractions, as in the Bhattacharya et al paper? These might to be very different from the BZ patterns and perhaps subject to\nmore specific regulation than the latter.\nThe data they found should certainly be published somewhere, but lack of a link, even hypothetical, to some possible\nmechanisms for producing the patterns and linking those mechanisms to cognition makes the paper less exciting. If they\nthink that possible decoding functions constitute the importance of the paper, that should be expanded upon."}, "author_response": "Reviewer #1:\nThe paper presents a study on traveling waves in humans during task and show that diverse\npatterns of traveling waves are not only present but also correlate with behavior. The conclusions\nare backed up by the results and methods look sound.\nThank you for the encouraging comments!\n1.1. It is worth noting that a lot of the work is based on previous studies, including some from the\nsame authors. This raises the problem that the novelty of the current study simply arises from\napplying their established approach to their new datasets.\nThank you for this comment. Please note that even though some of the methodology in our paper is\nbased on our prior work, there are several significant contributions beyond our prior work. First, our\nearlier studies (Mohan et al., 2024; Zhang et al., 2018) described only planar traveling waves, severely\nlimiting their ability to rapidly and dynamically reorganize to adapt to complex behaviors in humans.\nUsing a new analytical framework based on the localized circular-linear regression model and\nindependent component analysis, here we showed that human traveling waves are a much broader\nphenomenon that show a range of shapes, including spiral waves, source and sink waves, and more\ncomplex heterogeneous wave patterns, with more complex links to behavior that were not seen\npreviously. Our results thus provide an important contribution because they show for the first time the\nwide array of complex spatial wave patterns in the cortex that support cognition, extending beyond\nplanar waves. These results further suggest that the cortex exhibits heterogeneous patterns of neuronal\nactivity beyond the previously known anatomical and functional hierarchies, a significant contribution\nover earlier works. Second, in addition to distinguishing broad behavioral states, the complex spatial\npatterns of traveling waves also distinguished specific behavioral states such as English letters “G” vs.\n“H”, suggesting that traveling waves also correlate with specific item states. Showing traveling-wave\ncorrelates of specific neural network states is a major contribution because it shows that detailed patterns\nof neuronal network activity are differentiated by different large-scale patterns of electrical activity and\noscillations. And third, using deep learning models, we were able to robustly predict a subject’s\nbehavioral state at the single trial level based on the direction and strength of the observed patterns of\ncomplex traveling waves, suggesting that they can be useful for decoding and brain-computer\ninterfacing. Thus, the approach we presented here provide a general framework for analyzing a wide\nrange of spatial patterns of traveling waves and their association with behavior, and these methods are\nalso applicable to other data modalities such as scalp EEG, magnetoencephalographic, optical, as well as\nfMRI.\n1.2. I would like to see the figure legends be more descriptive. Often, I could not fully understand\nwhat was being plotted (e.g., Fig. 1i and SFig. 4c).\nWe have now added more descriptions to the figure legends wherever appropriate.\n1.3. Decoding numbers are relatively low, but I'm not sure if one would expect much better given\nthe info fed and model used. However, the authors should soften the claims of \"reliable\" and\n\"robust\" decoding.\nThank you for this suggestion. As suggested, we have softened our claims of “reliable” and “robust”\ndecoding throughout the manuscript.\n1.4. The description of modes and their relation to wave patterns could be clarified more, it feels\nlike modes are essentially capturing pattern types, but modes and pattern types are used on\ndifferent results back and forth.\nThank you for this suggestion. We would like to clarify that we have used the terms “modes” and “wave\npatterns”/“spatial patterns” interchangeably throughout our manuscript. This is because each of the\nmodes corresponded to one unique and independent spatial pattern extracted from the ICA algorithm.\nWe have now clarified this on page 23: “We will be using the terms “modes” and “wave\npatterns”/“spatial patterns” interchangeably throughout our manuscript. This is because each of the\nmodes corresponded to one unique and independent spatial pattern extracted from the CICA algorithm.”\nReviewer #2:\nThis is a technically rigorous, data-rich study that provides a detailed taxonomy of wave patterns\nobserved during memory tasks. A comprehensive set of controls was implemented to identify\ntraveling waves and link them to behavioral states. While the analyses are thorough and address\npotential confounds in depth, I believe the manuscript would benefit from a broader examination\nof the functional role of these patterns to better suit the broad readership of Nature\nCommunications.\nWe would like to thank the reviewer for the very encouraging and positive comments. Please find the\nresponses to your specific comments below.\n2.1. The reliance on pairwise classification appears methodologically restrictive. Although the\nauthors suggest that their analysis of traveling waves may yield novel insights into the brain’s\nspatial organization, it is not evident how binary discrimination alone substantiates this claim.\nComplex neural datasets frequently exhibit nonlinear decision boundaries and latent structure,\nsuch as the hierarchical organization observed in episodic memory tasks. Capturing such\ncomplexity may require targeted feature engineering and model architectures that explicitly\naccommodate multi-class or structured outputs. Moreover, performance in a binary framework\ndoes not imply that the learned feature representations will generalize to the multi-class setting\nnecessary for comprehensive behavioral decoding, which appears to be the broader aim of this\nwork. Approaches employing deeper architectures or ensemble methods may be more appropriate\nfor uncovering meaningful structure in the data. In the absence of such methods, it is difficult to\nassess the strength of the assertion that “behavioral states of human memory representations can\nbe decoded from traveling waves.” Given that MANOVA has already demonstrated statistically\nsignificant differences between states, the added value of the multi-layer network analysis remains\nunclear.\nWe thank the reviewer for raising these issues. The main goal of our decoding analysis was not to\npropose targeted neural network architectures, rather to demonstrate a proof-of-concept framework for\ntraveling waves which can be used as a basis for single trial neural decoding, and therefore can be\npotentially used as a tool to develop a novel brain-computer-interface. Please note that our decoding\nanalysis using multilayer neural networks, even though carried out for pair-wise cortical states, employs\nrigorous cross-validation procedures to ensure robust decoding and avoid any overtraining. We used the\nextracted weights from the ICA procedure as features for training our neural network classifiers and we\nincluded all the modes for classifying the behavioral states. Multilayer neural networks are capable of\nhandling highly nonlinear and high-dimensional decision boundaries and latent structure, such as the\nhierarchical organization observed in our complex spatial episodic memory tasks.\nAs per the reviewer’s suggestion of multi-class classification, we have now carried out additional\nanalyses where we have run our neural network analysis for simultaneous decoding of all behavioral\nstates, where we performed training on all the behavioral states together and then classified the\nbehavioral states on the test data using the same five-fold cross-validation technique. We were able to\nsuccessfully decode the behavioral states for both the spatial (p < 0.001, one-sided sign test, median\ndecoding accuracy ~ 40% versus chance at 16.66%) and verbal (p < 0.001, one-sided sign test, median\ndecoding accuracy ~ 15% versus chance at 12.5%) memory tasks in our multiclass decoding analysis.\nThese multiclass decoding results strongly support our conclusion that behavioral states of human\nmemory representations can be decoded from traveling waves.\nThe reviewer also asked about our rationale for the pairwise decoding approach rather than a multi-class\nclassification approach. Although these two approaches can give analytically similar results (Lorena et\nal., 2008), our rationale for focusing on pairwise decoding analysis in the manuscript is that, some pairs\nof cortical states were more decodable compared to other pairs. For example, Figure 10C shows that for\ncluster 7, the spatial wave patterns for the encoding and distractor periods were distinct from each other\nand hence were more easily distinguishable compared to those from the confidence and feedback\nperiods. Moreover, the navigation and distractor states were generally decodable from other states, but\nnot from each other, which indicates that these behaviors generally exhibit similar traveling wave\npatterns (Figures 10F, G). These important results will be missed from a multiclass decoding analysis\nwhere all the behavioral states are decoded together or, from the MANOVA significance tests where all\nthe behavioral states are simultaneously included for statistical analysis. Therefore, the focus of our\nmanuscript’s decoding is the pairwise approach, with rigorous cross-validation procedures to avoid\novertraining, which again strongly supports our conclusion that behavioral states of human memory\nrepresentations can be decoded from traveling waves. Moreover, to further justify that the binary\napproach is robust, we have now carried out additional statistical analysis by shuffling the labels of the\nbehavioral states such as encoding, retrieval, etc. and then running the same decoding analysis on this\nshuffled data to generate a distribution of shuffled data against which the actual decoding performance\ncan then be compared. This analysis again revealed that our decoding results are statistically significant\nacross all clusters (all p’s < 0.05), suggesting that structure of wave activity is the source of the\nclassification accuracy and behavioral states of human memory representations can be decoded from the\nstructure of the traveling waves. We have included these new results in the revised version of our\nmanuscript (page 31).\n2.2. Also, a key limitation is the omission of spike rate analyses. Claims that waves “dynamically\nstrengthen” or “transiently reorganize” functional connectivity or route spiking should be\nsupported by spike–field coherence or spike rate comparisons across the new wave patterns (e.g.\n\"complex waves\") considered here.\nThank you for this comment. We agree that spike rate or spike-field coherence analysis is interesting\nand would directly show how traveling waves modulate and are modulated by the spiking activity of\nneurons. However, such an analysis will require us to collect entirely new experimental data\n(simultaneous spike data and LFP) from new patients. We believe that this is beyond the scope of the\ncurrent manuscript.\nNevertheless, we have now acknowledged this in the Discussion section of our paper (page 14): “Future\nstudies with simultaneous single neuron and intracranial EEG recordings from neurosurgical patients\nperforming spatial and verbal memory tasks are needed to directly examine the functional relevance of\nthe spike-wave coupling. This will reveal how the complex spatial patterns of traveling waves modulate\nand are modulated by the spiking activity of the neurons.”\n2.3. I would appreciate a more detailed discussion of the potential mechanisms underlying the\npredominance of one wave class over another. For instance, the observation that source waves are\nmore prevalent in the spatial memory task is intriguing and could potentially reflect changes in\ninhibitory dynamics and reorganization of the theta cycle between encoding and retrieval. Similarly,\nrotating waves might arise from spatial inhomogeneities, but it is unclear why such factors would\nspecifically influence the spatial memory task more than the verbal task (Fig. 4).\nWe would like to thank the reviewer for this suggestion. As suggested, we have now added detailed\ndiscussion related to the mechanisms that could underlie different types of traveling waves and their\nfunctional roles.\nThis new text, in pages 10-11, states: “An especially distinctive aspect of our results concerns concentric\nsource and sink waves. A previous study showed that source waves were common than sinks during\nsleep spindles in humans, indicating that sources putatively support complex cortical processes such as\nmemory consolidation and synaptic plasticity (Muller et al., 2016). Because they involve waves\npropagating outward from a single location, we hypothesize that the source waves we observed indicate\nthe presence of a small local group of neuronal assemblies that dominate information flow by routing\ntheir information outward from one local area in a direction towards widespread brain areas. Given\nthat we found an especially large prevalence of source waves in our spatial memory task, it suggests\nthat spatial cognition involves directed information flow propagating outward from local nodes. Spatial\nepisodic memory (navigating or recalling spatial layouts) requires broad integration of rich sensory\ninformation across multiple cortical and subcortical areas, which might favor more sources of waves\nsince these sources involve widespread initiation of oscillatory activity to broadcast information\noutward across multiple networks. In comparison to the spatial task, the Sternberg verbal working\nmemory task, which involves maintaining a few items on each trial that each come from a relatively\nsmall dictionary of possible items, primarily recruits fronto-temporal circuits rather than the larger-\nscale brain networks that are important for richer spatial memories. Therefore, the verbal memory task\nmight rely more on focused convergence of waves onto sinks for more localized cortical processing. In\nthe verbal memory task, inhibitory control is stronger to limit irrelevant spread; thus, in this context our\nresults suggest that the oscillations we observed are relevant because they suggest that neural activity\nconverges onto fewer stable sinks during encoding to help maintain fixed items in buffers. Our findings\nand these predictions regarding the functional role of traveling waves and spatial oscillation patterns in\nmemory are, of course, predictions based on our empirical findings, but they provide useful hypothesis\nfor computational modeling that may probe the importance of these source nodes for cortical processing\nand spatial cognition.\nMuller et al., 2016 also found the presence of rotational waves during sleep spindles. Rotational waves\nsuch as spirals revisit the same brain areas in multiple cycles, therefore sequentially scanning\ndistributed brain areas to help enhance synaptic plasticity for memory consolidation. We observed a\nhigher number of rotational waves during the spatial memory task compared to the verbal memory task.\nSpatial episodic memory requires bidirectional loops between the cortical areas and the subcortical\nareas. Rotational waves naturally emerge in recurrent loop topologies where excitatory signals re-enter\nthe system with delays. These spirals can integrate multimodal spatial signals such as landmarks,\ntrajectories, and contexts into coherent phase patterns, therefore providing a temporal scaffold for\nmultiplexing different feature streams around the spiral core and supporting flexible scene construction\nand navigation. These features might be less critical for short term verbal working memory which is\nmore maintenance-oriented, therefore might favor more simple, stable oscillatory attractors such as\nplanar waves with minimal dynamical complexity.”\n2.4. Last, the potential impact of array positioning on the observed results is unclear. For instance,\ncould the detected plane waves represent segments of larger rotating waves, and if so, how was\nthis possibility addressed?\nWe thank the reviewer for raising this important point. Since the larger 8´8 grids in our dataset span\nalmost one of the entire hemispheres, it’s highly unlikely that a detected plane wave is part of a larger\nrotational wave. But if the size of the ECoG grid is small, such as the 6´4 or 8´4 grids in our dataset,\nit is indeed possible that the detected plane waves could represent segments of larger rotating waves.\nThis is an inherent limitation of our methodology. Therefore, it is theoretically possible that the number of\nrotational waves that we detected in our dataset could potentially be undercounted, and we do not\nbelieve there is any analytical method for overcoming this issue.\nNonetheless, to address this issue transparently, we have now acknowledged this limitation on pages 24-\n25: “If the size of the ECoG grid is small, such as the 6´4 or 8´4 grids in our dataset, there is a\npossibility that the detected planar waves could represent segments of larger rotational waves.\nTherefore, the number of rotational waves that we have detected in our dataset could potentially be\nundercounted. Future studies with much larger ECoG arrays, or noninvasive whole brain recordings,\ncould be useful to provide a more accurate assessment of the number of rotational waves and their\nfunctional role during various memory tasks.”\nLikewise, could the observed complex wave patterns be artifacts rather than planar or rotating\nwaves?\nPlease note that complex waves were those which could not be classified as one of the following\ncategories: “planar”, “rotational”, or “concentric” (“expanding” or “contracting”). Complex waves were\nalso stable at the individual epoch level as shown in Figure 9 and also across epochs since these waves\nwere able to distinguish behavior, for example, see Supplementary Figure 19. Also please note that some\ncomplex waves showed a combination of multiple patterns, such as separate subsets of electrodes\nshowed planar, rotating, or expanding/contracting waves (Figure 9, Supplementary Figures 9-18). Group\nlevel statistics related to complex waves are in Figure 4 and Supplementary Figure 18.\nTo ensure that complex waves were statistically robust, only those complex waves that exceeded the\nvariance explained threshold were included in our analysis (Figure 4 and Supplementary Figure 18).\nSpecifically, we compared the variance explained for each mode with the theoretical variance threshold\n100/n, where n is the number of electrodes in an ECoG grid (please see Methods for details). We\nadditionally shuffled the electrodes in each cluster and recalculated the variance distribution across\nmodes and confirmed that the shuffled variance for all modes converged to the theoretical variance\nthreshold of 100/n. Also, please note that complex waves were subject to the same rigorous statistical\ntests as other wave shapes. These complex waves were statistically robust (all p’s < 0.001, Methods,\nFigures 1I, Supplementary Figure 2) ruling out the possibility that complex waves are artifacts.\nReviewer #3:\nThe paper studies how widespread brain regions coordinate their activity to support complex\nbehaviors, which is a fundamental question in neuroscience. The authors demonstrate the\nexistence of different types of travelling waves in the human cortex during memory tasks (planar,\nspiral, concentric, and complex traveling waves), and show how these specific wave patterns\ncorrelate with and can be used to decode distinct behavioral states. This is a major contribution.\nThe methodology is clear and robust, especially after the revisions, where various methodological\nissues were raised.\nMoreover, the authors showed that neither peak oscillation frequency nor SNR could be used to\ndecode behavioral states, strengthening their conclusion that the observed effects are genuinely\nrelated to the patterns of traveling waves.\nGiven the novelty of the findings, the methodological rigor, and the comprehensive response to\nprevious critiques, this paper makes a strong case for publication. It is likely to be of broad\ninterest to the neuroscience community and to stimulate further research into the functional role\nof complex traveling waves in cognition.\nThank you so much for appreciating the significance of our contributions and your very encouraging\nwords. Please find responses to your specific comments below.\n3.1. As a minor comment, in the the section about the generative mechanisms, the authors should\nalso refer to computational models that describe the appearance of frequency specific activation\npatterns (static or travelling) due to axonal time-delays, such as the normalized connectome by\nPetkoski & Jirsa 2022 Netw Neurosci, and modelling works for whole-brain activity such as:\nCabral et al. 2023 Nat Comm Biol and Petkoski et al 2019 PTRSA.\nSimilarly, important are the modeling works on the thalamo-cortical circuits which focus on the\ncortical geometry and were developed by the group of Peter Robinson and lately applied by the\nlab of Michael Breakspear.\nWe thank the reviewer for pointing out these important references, which we have now cited at\nappropriate places in our manuscript including on page 11: “These locally generated waves can\npropagate across widespread regions between locally connected neurons and interact with other locally\ngenerated waves, to generate complex patterns of propagating oscillations (Bao & Wu, 2003; Cabral et\nal., 2022; Huang et al., 2010; Jeong et al., 2002; Petkoski & Jirsa, 2019, 2022; Schiff et al., 2007).”,\nand also on page 13: “Subcortical generation of traveling waves also hints at the possibility of three-\ndimensional traveling waves coordinating neural activity between the cortical and subcortical brain\nareas, rather than waves propagating only along the cortical surface. Recent fMRI studies using a new\ntype of wave mode namely the geometric eigenmode have demonstrated that three dimensional\nthalamocortical traveling waves play a critical role in many higher order brain functions (Pang et al.,\n2023; Robinson, 2025). Denser sampling of depth electrodes in three-dimensions may help to better\nunderstand the full complexity of the generation and propagation of traveling waves across different\nshapes of cortex, perhaps answering the question about how waves propagate across sulci.”\n3.2. Finally, as a suggestion to both, the authors and the editor, I find it very difficult to read the\nwork in this kind of format, with the figures and captions separated. I strongly suggest including a\nreader-friendly file with the figures and captions in the text.\nFor the reviewer’s convenience, we have now uploaded an additional article file namely\n“Reader_Friendly_File”, where the figures and captions are in the text.\nReviewer #4:\nWhat I would have liked to see is a deeper discussion of why the authors think that the new\npatterns they found, (e.g., the targets and spirals) actually contribute to cognitive processing, as\nopposed to being associated with processing. One of the authors is a world class expert on such\npatterns and may have thought more about this than is in the paper. I would also like to have seen\na deeper discussion of the cited modeling paper by Bhattacharya et al, since the authors of the\ncurrent paper say that they reproduce the results of that paper experimentally. Do they think\ntheir waves come from thalamocortical interactions, as in the Bhattacharya et al paper? These\nmight to be very different from the BZ patterns and perhaps subject to more specific regulation\nthan the latter.\nThe data they found should certainly be published somewhere, but lack of a link, even\nhypothetical, to some possible mechanisms for producing the patterns and linking those\nmechanisms to cognition makes the paper less exciting. If they think that possible decoding\nfunctions constitute the importance of the paper, that should be expanded upon.\nWe thank the reviewer for the suggestion related to enhancing the depth of the Discussion on the\nmechanisms of the waves and their association with cognitive processing such as the alphabets.\nThe reviewer wants us to highlight the novel contributions related to the complex wave patterns,\ntherefore, we want to briefly highlight the novel contributions and the functional relevance of the\ncomplex spatial patterns of traveling waves. In our paper, we have shown the presence of new types of\ncomplex spatial patterns of traveling waves which are associated with both broad and specific behavioral\nstates, with a specific spatial pattern correlated with a particular behavioral state. Therefore, and as we\nhave shown in our paper using deep learning models, the behavioral states can be decoded using these\ncomplex spatial patterns of traveling waves. This is a significant advancement over prior studies which\nused planar waves and their association with human behavior. Since planar waves travel in only one\ndirection, it severely limits their ability to rapidly and dynamically reorganize to adapt to complex\nbehaviors in humans. Our results show that the cortex exhibits heterogeneous patterns of neuronal\nactivity beyond the previously known anatomical and functional hierarchies from connectomics.\nWe have now modified our text to further clarify the generative mechanisms of the traveling waves on\npage 11 of our revised manuscript: “Our findings converge well with theoretical predictions from\nbiologically plausible neural models based on weakly coupled oscillators (Bhattacharya et al., 2021;\nSato, 2022). These models hypothesize that complex patterns of waves can be generated locally based\non the initial spatial activation of neurons, where each neuron is connected to a few of its neighbors,\nwith distance dependent axonal delays in the order of conduction along unmyelinated horizontal fibers\n(Davis et al., 2021; Destexhe, 1994; Ermentrout & Kleinfeld, 2001). These locally generated waves can\npropagate across widespread regions between locally connected neurons and interact with other locally\ngenerated waves, to generate complex patterns of propagating oscillations (Bao & Wu, 2003; Cabral et\nal., 2022; Huang et al., 2010; Jeong et al., 2002; Petkoski & Jirsa, 2019, 2022; Schiff et al., 2007).\nIn addition to the coupling functions, a key determinant of the shape of these wave patterns is the\npresence of local shifts in the frequency of local oscillations. In locally coupled oscillator networks,\nwaves tend to propagate away from the cortical locations with the fastest intrinsic oscillation\nfrequencies, following a gradient to the locations with the slowest oscillations (Kopell & Ermentrout,\n1990). Thus, local activations of strong or fast oscillations can have a strong influence on the global\ntopography of these traveling waves (Bhattacharya et al., 2021; Huang et al., 2010; Sato, 2022). Using\nthis model, based on the spatial location of the locally generated wave sources and their relative\nfrequencies, a wide range of stable, regular (non-random) wave patterns can be generated (Kopell &\nErmentrout, 1990).”\nWe have now also specifically discussed the different generative mechanisms of the specific patterns of\ntraveling waves on pages 11-12: “For example, a linear spatial gradient in frequencies, such as\nrelatively higher frequency on the left side of the ECoG array compared to the right side, can lead to a\nplanar traveling wave going from the left to right on the ECoG array (Ermentrout & Kopell,\n1984). Therefore, a presence of a planar traveling wave for a specific English letter would imply this\nfrequency heterogeneity in the left of the ECoG array compared to the right side represents that specific\nEnglish letter. Similarly, bullseyes or target waves can also be generated by frequency gradients using\nnearest-neighbor coupled oscillators (Ren & Ermentrout, 1998). For example, a source wave in an\nECoG array can be generated when the frequency at the center of the array is relatively higher\ncompared to all other places of the ECoG array, constituting a wave diverging away from the center of\nthe array. Therefore, a presence of a source wave for specific English letter indicates relatively higher\nfrequency at the center of the ECoG array responsible for cognitive representation for that specific\nletter. Similarly, a sink wave in an ECoG array can be generated when the frequency at the center of the\narray is lower compared to all other places of the ECoG array, constituting a wave converging to the\ncenter of the array. Therefore, a relatively lower frequency at the center of the array for a specific letter\nwould indicate that such frequency lowering is causally responsible for cognitively representing that\nspecific letter. Recently, in a new work, we have found the effect of the frequency heterogeneities on\ngenerating these planar and concentric traveling waves and when these waves are combined with\nrotational traveling waves, can produce a wide range of complex traveling waves (Das et al., 2025).\nMore broadly, our results of distinctive spatial wave patterns for separate memory items suggest that\nthere are spatial patterns of excitation across the cortex for separate memory items during cognition,\nsimilar to those seen in the hippocampus (Mickiewicz et al., 2025).”\nOur revised paper also describes the functional importance of spirals by showing that the existence of a\nspiral corresponds to a small group of electrodes at the core of the spiral that have large phase advances\nthat is functionally relevant for representing an English alphabet as we note on page 12: “Similarly, the\npresence of spiral waves is also informative about the spatial arrangement of the cortical pattern that\nsupports distinct behavioral states. Certain arrangements of spatial heterogeneities in an excitable\nmedia can produce spiral waves. For example, spatial heterogeneity in a 2×2 block of electrodes at\nthe middle of the ECoG grid where the phase advances by approximately an even ± step along each\nedge, such that the four edges together complete an approximate ±2𝜋 winding around the block, can\nproduce a spiral wave, with the 2×2 block of electrodes at the center constituting the spiral core\n(Greenberg & Hastings, 1978). Therefore, a presence of a spiral wave for a specific English letter\nindicates a spatial heterogeneity in a 2×2 block of electrodes at the middle of the ECoG array which\nconstitutes a spiral wave causally responsible for cognitively representing that specific letter.”\nThe reviewer asked about the Bhattacharya et al. paper, however the goal of our manuscript was not to\nreproduce the results of the Bhattacharya et al. paper experimentally. Rather, we mentioned this work to\ndiscuss the interesting question of how three-dimensional waves could be generated in the human\ncortex. We now describe this issue more clearly on page 13: “An additional way that traveling waves\ncould effectively propagate through sulci is by generation in subcortical areas, such as thalamus, as\ndemonstrated by recent computational modeling work (Bhattacharya et al., 2021) and experimental data\nin rodents (Ye et al., 2023). Subcortical generation of traveling waves also hints at the possibility of\nthree-dimensional traveling waves coordinating neural activity between the cortical and subcortical\nbrain areas, rather than waves propagating only along the cortical surface.” It is interesting to\nspeculate about the role of the thalamus in traveling wave generation. However, our dataset did not have\nthalamic recordings because typically, neurosurgeons do not put intracranial electrodes in the thalamus.\nTo directly show the potential impact on the cortical traveling waves from the thalamus, additional data\nneeds to collected from new epilepsy patients by running the same spatial episodic and verbal memory\ntasks, however, we think that such an endeavor is beyond the scope of the current manuscript and we\nreserve such work for future exploration."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have reviewed the revised manuscript and confirm that my previous concerns have been satisfactorily addressed—\nspecifically regarding the need to soften claims about robust decoding and to expand figure legends for improved clarity.\nI have one minor point requiring clarification in the final version (no further review needed): In the \"Electrophysiological recordings and preprocessing\" section, you note that different patients had different electrode configurations based on their\nclinical needs. However, when electrode locations are computed in standardized Talairach coordinates, the visualizations\nconsistently appear to show square grid patterns. Please clarify how this standardization produces apparently uniform grid\narrangements despite the variable electrode configurations across patients.", "reviewer_2": "The authors have addressed my comments.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.\nThe authors have responded to all of my comments and I suggest accepting the paper for publication.", "reviewer_4": "I agree with the other reviewers that this is a very interesting paper worthy of publication. However, it makes claims that,\nas far as I can see, it doesn’t support. Consider the abstract:\n“A fundamental challenge in neuroscience is explaining how widespread brain regions flexibly interact to support behaviors.\nWe hypothesize that {{traveling waves of oscillations are a key mechanism of neural coordination}}, such that they propagate\nacross the cortex in {{ distinctive spatial patterns that control how different regions interact. }} To test this hypothesis, we used\ndirect brain recordings from humans performing multiple memory experiments and an analytical framework that flexibly\nmeasures the propagation patterns of traveling waves. [[We found that traveling waves propagated along the cortex in not\nonly plane waves, but also spirals, sources and sinks, and more complex patterns. The propagation patterns of traveling\nwaves correlated with novel aspects of behavior, with specific spatial patterns reflecting particular cortical processes and\neven individual remembered items. Our findings show that large-scale cortical patterns of traveling waves reveal a diverse\nrange of spatial patterns in the brain]] that {{underlie neural computation}} and [[are relevant for neural decoding.]]”\nThe text in [[ ]] is what they did, and what they did is an excellent contribution. The text in {{ }} is what they claim, and I still\ndon’t see any evidence for that. They did give some helpful hints of how the waves might help with coordination, but they\ncertainly did not establish that their ideas are correct (or in some cases even plausible to me). I would be happy to see the\npaper published if the claims that go beyond the decoding and correlation are removed or at least that they remain\nhypotheses. That is, the paper needs to be clearer about the distinction between what they found and what they think may\nbe the consequences."}, "author_response": "Reviewer #1:\nI have one minor point requiring clarification in the final version (no further review needed): In\nthe \"Electrophysiological recordings and preprocessing\" section, you note that different patients\nhad different electrode configurations based on their clinical needs. However, when electrode\nlocations are computed in standardized Talairach coordinates, the visualizations consistently\nappear to show square grid patterns. Please clarify how this standardization produces apparently\nuniform grid arrangements despite the variable electrode configurations across patients.\nThank you for this question. Please note that in our analysis, we only included the subdural\nelectrocorticographic (ECoG) grid electrodes which are manufactured in a uniform grid arrangement.\nThis has been included in our revised manuscript in page 16: “Patients were implanted with different\nconfiguration of electrodes based on their clinical needs, which included both electrocorticographic\n(ECoG) surface grid and strips as well as depth electrodes. In this work, we only examined the ECoG\ngrid electrodes on the cortical surface.”. Also please note that square grid configurations such as 8*8,\n6*6, etc. are not the only type of grid configuration in our dataset, rather we had a very diverse range of\ngrid arrangements such as 4*8 grid, 4*6 grid, 4*5 grid, etc. (Supplementary Figure 1), capturing the\ndifferent clinical needs of the different patients.\nReviewer #4:\nI would be happy to see the paper published if the claims that go beyond the decoding and\ncorrelation are removed or at least that they remain hypotheses. That is, the paper needs to be\nclearer about the distinction between what they found and what they think may be the\nconsequences.\nAs per the reviewer’s suggestion, we have now removed the phrase “underlie neural computation” from\nthe Abstract. The Abstract has been revised accordingly as we note: “Our findings show that large-scale\ncortical patterns of traveling waves reveal a diverse range of spatial patterns in the brain that are\nrelevant for neural decoding.”"}]

Decision letter