Introduction

Traveling waves (TWs) have been described in human cortex since at least 1949 (1). They are a commonly described class of the more general phenomena of cortical phase dynamics, which we define as the set of all possible phase patterns. TWs may take the form of linear gradients (2–5), phase cones (6,7), spirals (8–10), and may also vary in velocity, i.e., ranging between standing and pure traveling waves (11,12). Large-scale TWs, specifically, are monotonic phase gradients measured over more than 8cm, up to the global scale of the cortical sheet.

Theoretical analyses have long predicted the existence of macroscopic wave dynamics (13– 15). Technically, the measured waves are group waves in the physical sense (16), that is, we measure the envelope of the heterogeneous activity at some temporal frequency (TF) as it moves through the cortex. The waves observed in neural mass models arise through a dynamic balance of excitation and inhibition and activated fields radiate outward from previously activated areas (13,14,17–19). Mathematical analyses indicate activity modes that are associated with damped TWs (14), and under appropriate boundary conditions, these correspond to large-scale TWs (see Figure 1). These waves are best understood as being coordinated by long-range myelinated axons (17). Being group waves, they do not directly indicate the velocity of information transmission (13,20), that is, the speed of transmission of axonal volleys (c.f. 18,19). They do, on the other hand, reflect the speed at which regions of activity can flip between different regimes, often as the result of small changes in coordinated synaptic inputs (20, c.f. 21).

Figure 1.

Alternative hypotheses have been proposed to explain large-scale TWs measured extra- cranially. We highlight three prior studies (21,22,25) here due their clear and explicit mathematical treatments, enabling an effective collision with the present research. Local TWs, driven by short-range intracortical connectivity, have been proposed as the origin of the extra- cranial measurements (21,22), as shown in Figure 1. Since these waves move slowly in the grey matter, but over a local region, when projected to the scalp they reproduce the fast velocity characteristics of the reported extra-cranial alpha waves (26–28). Another explanation for macroscopic TWs is in terms of localized oscillating cortical sources mixing at the sensor level (e.g., in electroencephalography, EEG). In this approach, two spatially offset and phase- lagged components are derived experimentally by source modelling, which in turn reproduces the measured extra-cranial waves (20,24; see Figure 1). Importantly for the present research, the dynamics associated with the alternative hypotheses can be re-represented as Fourier spectra of grey-matter activity, and these spectra give distinct and measurable predictions.

For the veridical TWs hypothesis, the expected peak in the spatial frequency spectrum of phase is at a low spatial frequency, while the alternative hypotheses for artefactual large-scale TWs predict peaks at higher spatial frequencies (SFs). In order to test these hypotheses, we therefore estimate the relative power of different bands of SFs comparing high (SF>16 cycles/m) to low (8<SF<16 cycles/m) SF contributions to the measurable phase; these SFs are approximately equivalent to wavelengths (𝜆 = 1 ⁄ SF) of λ<6cm vs. 6<λ<12 cm.

The alternative explanations that have been proposed for measurements of extra-cranial macroscopic TWs rely on two assumptions. First, for the localized activity to meaningfully reproduce the sensor level measurements, the activity proposed must dominate during the event; that is, the local TWs or oscillatory sources must arise within a relatively quiescent cortex. Otherwise, the local waves or oscillatory sources, as measured extra-cranially, will be washed out by the cortical activity that actually dominates. The second, and related assumption, is that measures aggregated over samples and/or trials accurately reflect relevant signals in single sample data (29). This aggregation assumption is related to the quiescence assumption because it allows local waves or oscillatory sources to be detected in otherwise ‘noisy’ environments. Typical sample aggregation procedures, such as averaging in event- related potentials, reverse correlation or computing the lead-field matrix in source localization, result in activity with high variation across aggregated quantities being made quiescent, thereby treating this activity as noise (29). By contrast, the present study estimates SF via low dimensional representations of the phase matrices using singular value decomposition (SVD), rather than for sample aggregates (e.g., trial means). SVD is a lossless decomposition of a participant’s phase data, up to the choice of cut-off rank for exclusion of noise-like phase vectors.

Previous research has consistently revealed that the spatial power of phase decreases monotonically with SF. Intracranial measurements show that the SF spectrum of phase declines with power from 12 cycles/m (8cm wavelengths) (2,7,8,29–31). The estimates, however, use cortical surface arrays with limited spatial coverage (≤8cm) and were therefore insensitive to longer wavelengths (see Figure 1). Likewise, animal studies also point to a SF spectrum for phase that monotonically decreases in power with increasing SF (6), but at measurement scales of 7mm or smaller (6,32–34). Estimates of the relative contribution of large-scale TWs to the SF spectrum of cortical phase dynamics have used EEG (7,30) (see Figure 1) or magnetoencephalogram (MEG) (11), and consistently describe a SF spectrum with maximum power at wavelengths near the array size, and decreasing in power with smaller wavelengths.

These studies, collectively, have not been able to unambiguously determine the relative importance of low SF components. Fourier analysis is sensitive to the sampling of the underlying data, both at the smallest range (the sampling frequency) and the maximum range (the size of the array) (35). Phase dynamics at the largest cortical scales, as broad-band signals, cannot be assessed using small intra-cortical arrays. Obversely, extra-cranial measurements are contaminated by volume conduction artefacts (EEG) or blurring by distance to the sensor array (MEG), both of which act as a low pass filter (18,29,36) (see Figure 1). A low pass filter will suppress higher SFs, thereby inflating the relative effect of the low SFs quantified from the phase dynamics and introducing an ambiguity. In the present research we use cortical depth electrodes to assess that part of the SF spectrum of phase between 8 and 50 cycles/m. Stereotactic EEG (sEEG) contacts have excellent spatial resolution and coverage (although sparse and irregular) which can extend to 25cm of the cortical sheet. By extending grey-matter measurements of the SF spectrum to longer wavelengths, we were able to accurately assess the relative contribution of large-scale phase dynamics to cortical activity.

The numerical methods described in the present research are novel. The spatial pattern of contacts in sEEG within the cortical grey matter has a high degree of irregularity and sparseness. In addition, the contacts are embedded in a thin, three-dimensional sheet. Viewed globally, the sheet covers the surface of a highly convoluted, three-dimensional, distorted spheroid. While methods exist for Fourier analysis on irregular grids (37–41), the combination of irregularity, sparseness and distortion in sEEG makes them unsuitable for the present research question. Methods also exist for computing spatial derivatives on irregular grids, but these generally assume denser sampling within a continuous volume (42). In the present research, we use a multi-scale, finite difference method to compute the SF of the phase, which is first estimated at each TF. To avoid spurious phase estimates, we first decompose the phase vectors into patterns of spatial covariance using SVD, retaining only spatial components that explain non-trivial amounts of variance.

In the present research we estimate the SF of the cortical phase dynamics. This means removing the amplitude information from the initial Fourier components estimated at each temporal frequency. This removal is for two reasons: (1) consistency with the research on TWs and other phase dynamics and (2) the assumption that the absolute amplitude measurable at each contact is a source of noise, since it depends on the exact details of contact placement which are not as exact as experimental studies in animals.

Despite the published research on SF spectra of cortical phase dynamics (6,7,30,43), decades of neuroscience research has favoured explanations of extra-cranial measurements in terms of mesoscopic TWs or focal cortical sources mixing at the sensor level (21,22,25). Here, we can provide a concrete test of the alternative hypotheses based on a simple criterion. We remove ambiguities in the measurement of purported macroscopic TWs by quantifying them within the grey matter. Consequently, the proposed mechanisms of localized TWs, oscillating sources or large-scale phase gradients can be tested against the SF spectrum of phase measurements for consistency with the proposed mechanism, as shown in Figure 1. Specifically, the characteristically small spatial scale of the proposed mesoscopic waves or oscillating sources will be revealed as maxima (or local maxima) in the SF spectrum. These high SF components have been hypothesized to combine, via blurring of signal, to form apparent low SF macroscopic TWs at the sensor level (21,22,25). This effect will therefore only be apparent if other cortical activity is relatively quiescent; the high SF components must really be a prominent part of the signal. Alternatively, a finding that low SF components dominate the spectrum of phase, at distances corresponding to large-scale TWs, cannot support the blurring or source mixing explanation for macroscopic TWs. This is because the sensor level measurements will in turn be dominated by the real low SF component of the cortical phase dynamics. Measurement of macroscopic TWs in EEG or MEG can only create a purely artefactual low SF gradient of phase, through mixing and blurring, if no such gradient dominates the underlying signal in the first place (compare Figure 1A & 1D).

Results

Overview

In the present research we analyze a publicly available data-set (RAM Public Data https://memory.psych.upenn.edu/RAM_Public_Data). We report the results from 23 participants. Note that most subsequent procedures (construction of spatial frequency spectra, statistical testing of link to measurement array size and TF) were carried out at the single participant level, so each participant comprises a test of the hypothesis under scrutiny. The sEEG contacts in the grey matter were re-referenced using the pooled average of the 50% lowest signal amplitude white matter contacts (44). A typical participant’s contacts are shown in Figure 2A. From the sEEG time-series, complex-valued phase was analyzed at 34 frequencies, using short time-series methods (2 cycle Morlet wavelets). Centre frequencies ranged from 1.0 to 97.0Hz, on a logarithmic scale. A segment of typical phase data with a centre frequency of 9.2Hz is shown in Figure 2B.

Figure 2.

The highly irregular sampling grid provided by the sEEG limits the statements that can be made about the cortical phase dynamics. Due to the inherent noise of the spatial array, it was not possible to assess the smoothness of the cortical phase dynamics in the present research. We therefore consistently refer to the SF of the cortical phase dynamics and separately discuss implications for TWs. Assessment of the ubiquity of smooth gradients of phase, TWs, has been made elsewhere (11,29).

The present research used linear algebra techniques to empirically construct a basis with which to decompose the phase signal, following and extending previous techniques (4,8,11). From this simpler set of empirical functions (45) we are able to estimate the spatial frequency spectra of the cortical phase dynamics. The general strategy for empirical construction of a basis using SVD is illustrated in Figure 2. In the present research, the phase data at each TF was entered into an SVD. The left singular vectors (LSVs) from the SVD comprise maps of spatial covariances in phase. There are important mathematical relations between SVD and Fourier analysis and for appropriate signals (see Methods) the SVD will empirically recover approximate sinusoidal basis functions.

In Figure 3, we present examples of the use of SVD to decompose signals into components approximating the Fourier basis. Figure 3A shows a one-dimensional random walk. Figure 3B shows how SVD decomposes this signal into pairs of sinusoids of increasing TF. These sinusoid pairs are quantitatively similar to the family of one-cycle Morlet wavelets (see 45 for a mathematical proof). Figure 3c shows an image of waves of water. We encoded the image as ‘time-series’ of phase at the dominant vertical frequency, then analysed the ‘spatial’ vectors of phase in the horizontal dimension using SVD. The LSVs are approximate sinusoids corresponding to the Fourier series in the spatial dimension. Figure 4A shows how the corresponding singular values of these LSVs (as amount of variance explained) can be read as a power spectrum, with each successively higher SF component explaining less of the variance in the phase encoded image. Figure 3F shows how these methods have been applied to MEG data (4,11,29) to extract global TW components of the signal, and components at higher SFs. Here, successive groups of LSVs from the MEG phase data show increasing wavenumber (approximately k=1, k=2 …). Figure 4B shows the SF spectrum computed for each of these wave maps using the method presented in the current research. The sum of these separate spectra represents the SF spectrum for the phase at this TF in the MEG.

Figure 3.

Figure 4.

The irregular sampling of cortical spatial coordinates in the sEEG poses a difficult problem for standard Fourier analysis, even for specialized techniques targeting non-uniform arrays (37– 39,42). Here, we make use of SVD to empirically decompose the phase data, before using multi-scale differencing on each LSV to compute the SF spectrum. The primary reasons for using SVD were (1) to prioritize patterns of high spatial covariance and so exclude noisy phase estimates, and (2) an assumption that the LSV vectors represent Fourier-like components of the signal, being a combination of the smooth sinusoidal bases recoverable from regular measurement arrays (45; and Figure 3) and low rank corrections to ‘noisy’ spatial sampling (40; see Methods). This assumption is supported via surrogate testing of our method using brain signals whose SFs is known (i.e. MEG).

The irregular sampling of cortical spatial coordinates inherent in sEEG was overcome by first decomposing the phase matrix (columns of contacts by rows of samples, at a given TF) into components that capture the largest portions of spatial covariance in the phase signal. The LSVs of the SVD (Figure 2C) were then further quantified using multi-scale finite difference methods to estimate the spatial frequency spectrum. The multi-scale difference methods consisted, first, of reorganizing the sEEG contacts into triplets corresponding to the vertices of approximately equilateral triangles within the cortical sheet. Between-contact distances were calculated as geodesic distances (47) in the cortical sheet to accurately reflect the underlying topology of the cortex. See Figure 5 and Supplementary Figure 1 for examples of geodesic paths on the cortical surface. These triangular areas were used to calculate the rate of phase change across each region of the LSV of phase. These techniques enabled the spatial frequency spectra of the cortical phase dynamics to be quantified for cortical regions of varying area, and this was done at each TF.

Figure 5.

Single participant tests

We first quantified the estimated spatial frequency spectrum for all single participants, across all temporal frequencies. Two participants’ results are shown in Figure 5. Lower temporal frequencies had more SF power than higher temporal frequencies (permutation test, p<0.05; Panel A, D). However, this is in large part to the better signal-to-noise ratio in the estimation of phase for the former (more time-samples in the 2-cycle Morlet wavelet) and similarly better signal-to-noise in signal measurement (higher signal amplitudes at low TFs). The effect is not necessarily a feature of higher TF signal per se. Panels B & E show the same spectra stacked over TF, and reveal the characteristic shape of the estimated spatial frequency spectrum. The spectra peak in power at a wavelength approximately equivalent to the size of the largest equilateral triangle defined on the measurement array, sometimes at a slightly lower SF. The spatial power of the phase dynamics decreased monotonically with increasing spatial frequency (permutation test over temporal frequencies, p<0.05). Within the range beginning at the maximum triangle size in the measurement array (equivalent to ∼8.5 and 11 cycles/m in these cases), the spatial power drops steeply, before flattening at higher spatial frequencies.

This effect was statistically significant in 19 out of 23 participants (permutation test over temporal frequencies, p<0.05; see Supplementary Figures 2 and 3). All participant’s mean SF slopes trended in the same direction, with participants in the top 50%, by numbers of triangles, all being significant at p<0.001. This suggests a practical cut-off value for number of triangles in the present method of n=1500 for being able to construct sufficiently low noise spectra to indicate a significant trend of SF versus SF power. This equates to around 50 grey- matter contacts. Only one of 782 distinct participant × TF regressions was significant with a slope in the opposite direction. We therefore conclude that SF of phase had the highest power at the lowest SF. This main result is consistent with previous studies of spatial frequency of cortical phase dynamics: the spatial power decreases with increasing spatial frequency (6,7,30,43). The present results extend the previously measured spectrum into scales equivalent to a large fraction of the cortical size.

Examples of the low SF phase dynamics are shown in Supplementary Videos 1, 2 and 3.

Surrogate tests

We used surrogate methods to perturb the sEEG phase signals with artificial signals whose ground truth SF was known. This perturbation was done over a wide range of relative weightings between real and surrogate signals. If the method accurately measures SF, addition of the surrogate signals will induce a peak at the SF of the surrogate signal, increasing in relative amplitude as the surrogate signal amplitude is increased. Examples of the surrogate signal are shown in Supplementary Figure 4. Figure 6A shows the results for one participant at 9.2Hz, perturbed by a surrogate anterior-posterior gradient at three different SFs. The shift in the peak of the SF spectrum is apparent at low levels of injection (1/3 ratio) of the surrogate signal, consistently toward the SF peak of the surrogate. As the perturbation increases in weight, the surrogate contribution to the spectrum becomes clearer. The perturbation is approximately at the correct SF when twice the strength of the empirical signal (2/3 ratio) and remains so until ratio 7/8.

Figure 6.

The surrogate SF signals were tested over 14 steps for a wide SF range (2c/m to 16c/m) and showed the same correct linear relation at each TF of empirical component of the phase. This was tested by regression of SF peak versus SF of surrogate (repeated for each surrogate direction and surrogate weighting) and was significant down to and including the 1/2 surrogate ratio (permutation tests over participants, p<0.001, ratio 1/3 not significant). Another example of surrogate tests on the same participant is shown in Supplementary Figure 4, here with a left-right gradient imposed over the same set of SFs. The Figure also shows cross-participant averages for these surrogate runs. Also shown are SF spectra with pure surrogate signal, giving a clear idea of the residual noise levels in the SF spectra estimation, since the only source of noise in this case is the irregularity of the sEEG array. Due to the highly irregular sampling of the cortex, the exact locations of the peaks for the veridical SF surrogate varied slightly from participant to participant for the same surrogate SF, but always trending within each participant in the correct ordering of SF peaks (peak SF vs. surrogate SF, permutation test over participants, p<0.001).

We also checked the accuracy of SF estimation by using the technique on neurophysiological data with known SF characteristics, here MEG. The method successfully recovered the distribution of SF power previously reported for phase signal (Figure 6B; compare Figure 2 in Alexander et al. 2016, and Figure 2 in Alexander et al., 2019). Additional checks were made on the method using MEG data. Figure 4 shows that each of the individual LSVs from SVD on MEG phase had different SF peaks using the method when computed separately, as expected from visual inspection of the LSVs (Figure 3F). This confirms that the peak SFs were being detected accurately, as the LSVs changed in SF characteristics. Figure 6B indicates how these separated curves sum to produce the completed SF spectrum in MEG.

However, we wanted to definitively rule out that the method was not just constructing distributions with left pointing tails, and that we were quantifying the left pointing tails as ∼1/f spectra. To this end we constructed an MEG time-series of phase in which the high SFs dominated. This was achieved by reversing the order, and so the weighting, of the first n singular values and reconstituting the phase matrix according to eq. 2 (Methods) and as described in Figure 2C. Figure 6C shows that this surrogate data had a SF spectrum that increased monotonically with SF over the relevant range, i.e., the method accurately characterized the new SF spectrum.

To additionally rule out whether the produced SF spectra were an artefact of the sparsity and irregularity of the sEEG array, we produced MEG time-series of phase in which the sensors were randomly removed until the phase vectors were only one third their normal length, i.e., 51 sensors from an original 151. This procedure was repeated for several random removal sequences. The results are shown in Supplementary Figure 5 and demonstrate that the effect of increasing measurement array sparsity and irregularity was to increase the noise of the SF estimate, as expected, but did not alter the overall shape of the spectra (wavelength vs. power, permutation test over runs, p<0.001). This procedure was repeated for the MEG surrogate in which high SFs dominated (i.e. SF power increasing with SF), with the same conclusions (wavelength vs. power, permutation test over runs, p<0.001).

We also tested for the effects of the average referencing used on the sEEG data (see Methods). Global synchrony appears as the DC component (SF=0) in spectra of phase. Computing an average reference increases the global synchrony of phase in an artefactual fashion. While this component is well separated from the minimum of the SF range of interest here, we needed to rule out the influence of the referencing method on our estimated spectra. To this end we produced spectra for MEG data, for which the case of a zero weight ‘average reference’ is also available (i.e., no reference). We compared spectra with an average reference added at increasing levels influence. Increasing the weight of average reference up to a factor of eight times the usual average reference did not disturb the conclusions of this report (see Supplementary Figure 6). This procedure was tested with the methodological step that removed negative SF components (see Methods) and without. The results showed that this step both lowers the DC component of the SF phase spectrum—and thereby ameliorates the effects of any artefactual global synchronous component—but also does not introduce a false maximum at the relevant low SFs due to lowering of the DC component.

Cross-participant analysis

We analysed the relation between wavelength and estimated power for each participant using linear regression. The regression lines for the logarithm of wavelength versus estimated power are shown in Figure 7, for all individual participants at 9.2Hz (permutation test over participants, p<0.001). The regression lines for all TFs are shown in Supplementary Figure 7 and all TFs showed a significant decrease in SF power with SF (permutation test over participants, p<0.001). The participants with the fewest triangles did not always show significant individual wavelength by SF power regressions at a given TF, but the regression lines almost universally trended (except one case) in the predicted direction. The present results extend the previous finding of ∼1/f spectra in cortex (43,48) to temporal frequencies up to 100Hz.

Figure 7.

The cross-participant analyses also confirmed that peak SF power declined linearly with increasing TF (logarithmic scale). This is show in Figure 7D for the participant-aggregate spectra (permutation test, p<0.001) and for the individual participant’s spectra in Figure 7E (permutation test over participants, p<0.001). This effect may be wholly or in part due to measurement and quantification noise at higher TFs, as previously mentioned. Importantly, the cross-participant analyses revealed no effect for TF on SF at peak SF power. This was true for participant aggregate spectra (Figure 7D) and for the permutation test over all individual participant’s regressions of (logarithm of) wavelength at peak SF power (Figure 7F). This means that the relative shapes of SF spectra were the same, regardless of TF. That is, the SF spectral power decreased monotonically over the relevant range from 8-12 to 50/cm, for TFs 1.0 to 100Hz.

The linear regressions we performed excluded SF estimates lower than the reciprocal of the maximum triangle size. This decision arises from considerations of Fourier analysis, in which the maximum wavelength detectable in broadband signal is a strict function of the maximum range of the (regular) grid size. In our results, the apparent SF power drops off for SFs lower than this, but this part of the spectrum is not meaningful. With the present method, maximum array size was weakly correlated with the SF at peak power for each participant (peak SF in TF aggregate spectrum vs. maximum inter-contact distance, permutation test, p<0.05, r=0.31). Maximum array size is defined by the two most (geodesically) distant contacts. It does not, therefore, always reflect the typical largest distances for the array. We describe in the Methods how peak SF is a function of the distribution of triangle sizes in irregular arrays and using the present SVD methods. For this reason, we use maximum triangle size as a bound for SF power regressions.

The array sizes used here did not reach the full extent of the cortex (∼45cm), however, we were at the outset interested in the peak spatial frequency for the case where the array size extended to the entire cortex. For visual convenience, the regression lines in Figure 5 and Supplementary Figure 3 are extended to larger cortical distances. The linear relationship suggests that triangle sizes extending to the full range of the cortical sheet will result in peak SFs estimates larger than the maximum reported for the present participants. This prediction is supported by previous analyses using EEG and MEG, showing that the peak spatial frequency coincides with the size of the whole head array (11,30), which in turn coincides with the entire length of the cortex accessible over the upper half of the skull. It would be informative to repeat the analyses presented here after careful selection of sEEG participants from all available databases for maximum cortical sampling extent, or as technical advances in cortical sampling further improve. We expect that participants with greater than 30cm coverage will extend the linear trend reported here.

Effects of TF power

The main analyses all used phase as the input variable, reflecting our interest in cortical TWs. However, the main conclusions are not disrupted when TF power is included in the calculations. Figure 8 shows the results when TF power is normalized on a contact-by-contact basis (see Methods). Individual participant’s spectra are largely unchanged, with a suggestion of a weak increase in power at higher SFs (>30 cycles/m) compared to lower SFs (<10 cycles/m). Supplementary Figure 8 shows the results of including (unnormalized) TF power. The clearest effect is the influence of decreasing TF power on SF power as TF increases (panels A and B) –TF power is small at higher TFs. We visually correct for this effect by plotting the relative values within the SF spectrum at each TF. Panels C and D show there is an apparent increase SF power between 12 and 20 cycles/m compared to the phase-only analyses, but not sufficient to overshadow the SF components 8-12 cycles/m. It is likely that much of this effect is due to noise in the measurement of signal amplitudes from contact to contact, since it largely disappears when power is normalized on a contact-by-contact basis. If a sensor array (e.g., a camera pixel array) is more sensitive at some sites than others, this will appear as an artefactual increase in higher SF components of the image.

Figure 8.

Discussion

In this research we aimed to quantify the low SF part of cortical phase dynamics spectrum, up to the range of large-scale TWs (wavelength>8cm), using cortical depth electrodes. We developed a novel method to fill this important gap in our understanding of macroscopic phase dynamics, addressing a current limitation in the literature. Previous estimates have either relied on extra-cranial measurements (3,4,11,27,28,49) or used smaller measurement arrays (≤8cm) (2,6,31,33). This means either the results were biased by a low pass filter (due to volume conduction in EEG or distance to measurement array in MEG), or, had a maximum measurable SF below the large-scale range, respectively. We show that the spectral power increases with wavelength, up to the limit imposed by the maximum distances in the recording array. The present results are therefore consistent with the existence of large-scale TWs in the cortex (Figure 1D) and are inconsistent with the alternative hypotheses for the measurement of these waves extra-cranially (Figure 1E & F). Low SF activity dominates the cortical phase dynamics; most of the spatial power is at the longest wavelengths. In other words, most of the power measured extra-cranially is therefore produced by long wavelength activity in the cortex.

Considered from the perspective of locally measured phase—i.e., local to a single grey matter contact—the measured SF spectra imply that most of the variance in phase is a function of low SF phase dynamics. Note that this statement does not include cross-frequency effects, outside the (relatively broad) range of the frequency bandwidth of the phase estimator used in this study. The present result partly answers the longstanding question of why there is so much variability in locally measured cortical responses, even to the same stimulus (45, c.f. 46). We suggest that low SF phase dynamics index a distinct scale of cortical activity and, judged from their dominance of that activity, comprise an important scale for understanding cortical function (52). The variability is a mystery if the response in a functionally relevant area is considered in the absence of the global phase dynamic context. Our results suggest the variability can be controlled for if the global context is also measured.

The shape of the SF spectrum in the range ∼8 to 50 cycles/m is therefore consistent with reports of the spectrum at smaller cortical scales (6,7,30). The cortical SF spectrum is generally described as 1⁄𝑓^𝛼, where the exponent 𝛼 is a function of the type of measurement (e.g., grey matter, cortical surface or extra-cranial) (6,30,43). Since peak SF also increases with the total size of a participants’ sEEG array (maximum inter-contact range), we assume that the 1⁄𝑓^𝛼 distribution of SF extends to the maximum size of the human cortex; this is approximately 45cm linear distance within the unfolded grey matter. This would be consistent with spectral peaks at lowest spatial frequencies measured using EEG (7,30) and MEG (11), since this is what is expected when a low pass filter is applied to a signal dominated by low spatial frequencies.

A surprising finding was that the shape of the SF spectrum did not vary much with TF. This finding implies that low SF cortical phase dynamics arise at temporal frequencies up to 100Hz, the maximum TF tested here. Some care should be taken in interpreting the SF spectra at temporal frequencies above 60Hz, since band-pass filtering alone is not sufficient to remove the low temporal frequency components of spiking activity (53). In addition, while volume conduction effects arise in local-field potentials (54), these are limited to 1cm range and so smaller than the minimum 2cm range (50c/m) reported here. These effects cannot explain the relatively uniform shapes of the SF spectra with changing TF.

Gamma (30-60Hz) and high gamma (>60Hz) synchrony have been functionally implicated across a range of scales in the cortex. Several studies have quantified the functional role of gamma band oscillations within cortical areas, notably in the visual and motor areas. Spike triggered analysis of local field potentials in V1 has shown that TWs have a greater range of propagation under low contrast visual stimulation, compared to high contrast (55). In multiple visual areas, Womelsdorf et al. (56) demonstrated that gamma phase synchrony between recording sites precedes high amplitude correlations between those same sites. V1 oscillations were analysed in the range 50-80Hz using the transfer entropy of the phase and show the activity to be compatible with directional wave propagation (57). Blakely et al., (58) analysed temporal covariance of gamma band activity recorded in ECoG arrays to map finger flexions in the motor cortex.

At a larger spatial scale, it has been demonstrated that the (putative) causal influence of frontal regions on more posterior zones using the phase slope index, across a wide range of temporal frequencies, including gamma (59). sEEG recording has been used to analyze long-range high frequency phase synchronization (60). Using zero- and nonzero-lag phase measures, these authors demonstrate that phase locking occurs in the high TF ranges 100 to 450Hz, across globally ranged cortical regions. None of these studies, however, specifically assessed the SF spectrum of the phase signals.

The existing literature on gamma-band TWs is patchy but consistent with our findings here. Macroscopic TWs are commonly reported in the delta to beta ranges (3,4,27,29,61). However, TWs have also been reported for low gamma frequency (40Hz) using MEG, in an event-locked paradigm (48). Similarly, EEG-based estimation of SF has been reported for temporal frequencies from 7 to 63Hz, with similar spectra (1⁄𝑓^𝛼) for spatial frequencies over this entire range (43).

We outlined in the introduction the theory of cortical TWs as group waves (13,14,17), which means the group wave velocity can exceed the speed of information propagation. Specifically, the rise and fall of coordinated cortical activity can outpace axonal conduction times, since it depends on patterns of spikes that have previously arrived (23). If TWs are group waves and do not directly reflect the speed of information propagation, this may appear to support the argument that the cortex’s electric fields are epiphenomenal, and that spiking activity of individual neurons is the ground truth. However, the literature on TWs shows that they can influence spiking patterns of individual neurons (33,53,54 cf. 35). To this we may add reports of the endogenous electro-magnetic field altering firing patterns of neurons (63,64). It seems likely that a major driver of cortical plasticity is precisely the dominant SFs found in the cortex, since post-synaptic membrane depolarization is the second major ingredient of long term potentiation (65). An evolutionary argument can also be applied, since dominant low SF phase dynamics will be more metabolically expensive than otherwise quiescent, mostly locally activated cortex. This implies critical function rather than epiphenomenon. Taken together, these arguments suggest that the pattern of phase organization in the cortex has a causal role and is not the epistemological equivalent of heat generated by a CPU.

The present results have important methodological implications. Considerations from Fourier analysis raised in the introduction indicate that the maximum spatial sampling range is critical to setting the lower bound of detectable SFs. This would not be the case if the SFs were pure and isolated, but for broadband, highly mixed signals it is. These considerations are supported empirically by the present results. The peak SF increases in wavelength with maximum contact distance for the participants (here, 0.12 to 0.24m), which in practice is related to the number of contacts, and maximum triangle size used in the present methods. In addition, we showed that there was a linear relation between peak SF and the triangle size used to estimate SF. The present findings also explain why detected TWs, regardless of cortical scale, generally have a wavelength ∼2π (29,66). If the spatial spectrum is monotonically decreasing in power over the entire large-scale to microscopic range, then the highest power waves detectable on the array will be those with a wavelength approximately the size of the array.

We conclude that reports of TWs should include estimates of SF to allow or disallow comparisons between studies with different extents of cortical coverage, since they may be accessing different SF ranges and therefore different functional connectivity ranges (c.f. 31,32). The present results, along with considerations from Fourier analysis, indicate that only methods that have large-scale coverage of cortical activity can assess the dominant scale of cortical phase dynamics. Unfortunately, source localization techniques, designed to improve upon sensor-level measurements, remove long range spatial correlations in the signal; both artefactual and real (18,68). This excludes these techniques from accurately estimating the SF spectra of phase dynamics since they filter out precisely the most dominant components, according to the present results. Similar considerations apply to methods to remove volume conduction effects from local field potential recording, such as taking the first- or second-order spatial derivatives.

For this reason, we argue, forward models of cortical activity are better suited to discovery of cortical dynamics from extra-cranial measurements (69–71). The present results suggest that the key methodological assumption of source localization techniques, namely that localized peaks in activity can be quantified to produce a mapping from extra-cranial patterns of activation to cortical locations, may not be sustainable. The contrary case, that high SF peaks exist in the cortical spectrum in real time (or that the low SF dominance is functionally irrelevant) should be demonstrated rather than implicitly assumed. Worse, confirmation of our present results would more broadly suggest that the endeavour itself is less meaningful than previously appreciated. If high SF activity forms only a small part of the total SF spectrum, then source localization is not just an ill-conditioned problem, but perhaps not the problem. At the least, estimates of relative amounts cortical activity explained, by sources thus localized, should be attempted (see 72,73 for some discussion of these issues).

TWs have been shown to allow prediction of future phase at single recording sites two wave cycles into the future. This effect was shown for both macroscopic MEG and ECoG (∼8cm range) measurements (8). The ability of TWs to predict future local activity shows that the long-range spatio-temporal correlations inherent in the structure of low SF activity are functionally significant. Future research could extend this approach using high temporal resolution sEEG to detect individual spikes and to characterize large-scale phase dynamics and to predict the former from the latter. More generally, most of functionally-related results concerning large-scale TWs have been obtained using extra-cranial methods (27,28,49,61,74). The reported linkage between late event-related potentials and TWs (4,29,75,76) means the body of research into late event-related potentials can guide experimental verification of the functional significance of macroscopic TWs at the single trial level.

Our conclusions regarding the dominance of large-scale cortical activity are consistent with other recent findings in neuroscience. There has been a resurgence of interest into the global structure of the cortical functional map, driven by developments in fMRI, large-scale open-access databases, and numerical analysis techniques (15,77,78). An intrinsic functional coordinate system has been proposed: along with hemispheric lateralization, a global coordinate axis has been described in terms of a visual to somatosensory/motor dimension, beginning at the occipital pole and calcarine sulcus, extending around the cortical hemispheres through the temporal lobe, frontal areas and to the central sulcus. Another dimension, at approximately twice the SF, is a function of the cortical distance between the primary areas and transmodal areas (77). Eigenmode analysis of fMRI data shows that macroscopic wavelengths dominate the spatio-temporal activity at time-scales characteristic of fMRI (15). This is consistent with previous findings in the EEG and MEG showing that long wavelengths are ubiquitous features of the spatio-temporal activity in the ranges delta to beta (27,29,31,67,79). Here, we have shown that low SF phase dynamics dominate the spatio- temporal organization of cortical activity measurable in the grey-matter, and that these global dynamics extend to fast cortical time scales into the high gamma range, at 10’s of milliseconds.

Methods

Data acquisition

We summarize the relevant details of the RAM data set here (35; https://memory.psych.upenn.edu/RAM_Public_Data). Full details of data acquisition are reported elsewhere (80). 35 human participants were part of a multi-centre sEEG monitoring study. The participants were drawn from experiments ‘FR5’ and ‘catFR5’, and we included those that had complete data sets for present purposes (time-series, contact locations, cortical surface meshes). From the initial pool of participants, 23 had sufficient numbers of grey-matter contacts to enable further analysis (see Numerical methods for criterion). Mean age was 34 years, nine female, four left-handed.

The participants’ task was designed to assess the effects of electrical stimulation on memory- related brain function. Each participant engaged in a delayed free recall task in which they studied lists of words for a later memory test. Trials with stimulation were not used in the present study, and our interest was in the general characteristics of cortical activity, not especially dependent upon this task. The research protocol was approved by the institutional review board at each hospital and informed consent was obtained from each participant. Electrophysiological data were collected from electrodes implanted within the brain parenchyma. Participants had temporal lobe depth electrodes, with frontal and parietal regions targeted as needed to best localize epileptogenic regions. Contacts on each electrode had a 5mm inter-contact spacing. The sampling frequency was 1000Hz, except for 8 participants, for which it was 500Hz. The number of grey-matter contacts ranged from 31 to 108 with a mean of 58.

The data for each participant was denoised by first calculating the standard-deviation of signal over samples, for each trial and sensor combination (44). If the mean standard-deviation for a trial was greater than ±3 the mean standard-deviations from all trials, then that trial was removed from further analysis. Then, if the mean standard-deviation (after recalculation) for a sensor was greater than ±1.5 mean standard-deviations from all sensors, it was removed from further analysis. Finally, if the standard-deviation (after recalculation) for a trial by sensor combination was greater than ±3 standard deviations from all trial by sensor combinations, it was replaced with the average signal (average calculated after excluding all over threshold trial by sensor combinations).

Using the automated contact labeling supplied within the dataset (checked and corrected manually by clinicians), contacts were partitioned into three groups: grey-matter, white-matter low-amplitude and other white-matter. We defined white-matter low-amplitude contacts as those whose standard deviation of signal amplitudes were in the bottom 50% of white matter contacts. The white-matter low-amplitude contacts were pooled to create an average reference for the grey-matter contacts (44). All numerical procedures were trialed using notch filters to remove the 60Hz line signal and its harmonics, and repeated without notch filters simply by avoiding these bands. No numerical differences were obtained for the common frequencies in the two versions of the data analyzed, so the former analyses are reported here.

Numerical methods

Phase estimation

Phase is represented as a unit-length complex number, ϕ.

We used complex-valued phase (8,81,81) in conjunction with the complex SVD to avoid the (error) noise introduced by spatial unwrapping phase when using real-valued techniques (c.f. 2,4,26,36,37). Phase is calculated from the sEEG time-series, for a particular TF, using short time-window Morlet wavelets (2 cycles, (84)). Short time series methods were used because previous research has shown TWs to often last only one or two cycles, and windows of 5 or 7 cycles can wash out the detail of fast temporal changes (4). Use of short-time window Morlet wavelets increases the temporal resolution of the phase estimate, at the cost of a decrease in TF resolution. For example, the wavelet with centre TF at 4.5Hz has a half-power bandwidth spanning 2.5 to 6.8Hz. This lower TF resolution has positive consequences, since it makes the analyses here-in more robust to small variations in TF across areas and regions.

We estimated phase for 34 centre-TFs, logarithmically spaced from 1.0 to 97.0Hz. Due to the low TF resolution of the 2 cycle Morlet wavelets and number of frequencies extracted, this amounted to an oversampling of the frequency range. We use indices of s and t to refer to measurement sites and sample times, respectively, giving a matrix of complex-valued phases for each TF, 𝚽_𝑠𝑡.

Because Morlet wavelets have poor frequency resolution, the very short time-window Morlet basis is only approximately sinusoid and robust to relatively wide variations in TF (85) and violations of sinusoidal shape. In practice, any surge of activity up and down with a time-course near the centre TF will have a high correlation with the applied Morlet basis. We are therefore confident that the present results are not confabulated by ‘ringing’ effects (86).

A potential criticism of the approach used here is that we estimate spatial frequency spectra for phase, while discarding amplitude information. The primary reason is our interest in TWs, which are generally defined in terms of phase. It is worth noting that single trial phase is a stronger determinant of mean signal amplitude than is single trial amplitude (29). This has been shown for EEG, MEG and ECoG. In other words, aggregate measures such as event- related potentials, coherence, and even values for the lead field matrix in beamformer techniques, are predominantly a function of phase interference over aggregated samples rather than true sample-by-sample variations in source amplitude (29,79). Conversely, if the amplitude recorded at contacts is sufficiently small, the estimated phase values will not affect the estimates of spatial frequency produced here; they will be noise-like and not contribute to the dominant components of the SVD. We introduce the effects of TF amplitude as an addendum to our main analyses.

SVD for empirical Fourier decomposition

Fourier methods are commonly used in scientific research to analyze signals that extend in time and space. The methods use sinusoids as a basis to decompose arbitrary n-dimensional series into weighted sums of sinusoids. Fourier methods have been used to analyze TWs in the EEG (3,87). The non-uniform spacing of the measurement array in the standardized montage (e.g., 10-20 system) is overcome by interpolating a regular grid between measurement sites. In sEEG, however, the contact spacing is irregular to a much larger degree, making interpolation not an option. Fourier analysis with multi-resolution grids is still an active area of research (42,88), yet these approaches are limited to a few scales of grid and do not account for irregular sampling. The analyses provided by nonuniform discrete Fourier transform methods (37–39,41) are also not suited to the extremes of sparsity and irregularity of the grey-matter sEEG contacts analyzed here.

SVD is a standard method to reorganize data into a more interpretable form. SVD has recently become a more common tool in neuroscience (15,89–94); for an alternative view, see (46)). Here we used SVD primarily to discount the influence of phase estimates that do not contribute strongly to the spatial covariance of phase. Also, the LSVs can be interpreted as Fourier components with corrections for the irregularity of the measurement array.

We have previously shown that SVD can be used to decompose phase data obtained via EEG, MEG or ECoG (4,29). The measurement coordinates of the scalp EEG or MEG sensor array define the (partial) surface of a three-dimensional spheroid. The SVD decomposes the phase data into families of phase gradients projected onto the surface of the spheroid. As illustrated in Figure 3F, the first three LSVs extracted using SVD are approximately equal in wavelength to the size of the measurement array. Since the LSVs were smooth gradients of phase, this meant these first three components were large-scale TWs. The first three components were (for EEG) anterior-posterior gradient, inferior-superior gradient and left-right gradient (4). In addition, macroscopic TWs with arbitrary directions could be described using weighted sums of the components, e.g., a mixture of the anterior-posterior and left-right gradients gives a diagonally propagating macroscopic TW. The next three LSVs were approximately wavenumber equal to two (i.e., twice the spatial frequency), and each successive three LSVs increasing in wavenumber by one (29).

We can draw out the relation between SVD and Fourier analysis more explicitly. Consider phase data with a smooth spatial organization, where the peak in the spatial frequency spectrum is at long wavelengths (i.e., λ≈size of measurement array), and the spectral power monotonically decreases with increasing spatial frequency. A matrix of a phase signals, A, then decomposes the many samples of spatial vectors of phase into a series of LSVs, ur in descending order of singular values, σr, considering only the r highest singular values (95).

The effectiveness of the SVD is that it re-organizes the data into components that are uncorrelated with each other (the LSVs), but that explain known amounts of variance in the data (the singular values). Because the correlation is zero, the components can be treated as independent observations of a system. Weighted sums of these components reproduce the original data in a lossless manner. Alternatively, only the most important components can be kept to more concisely describe the original data, while the ignored components can be treated as noise. One of the main reasons to use SVD in the present research is the reorganization of data into components that explain most of the variance, and others which are noise-like. This means the SVD discounts phase measurements that do not create consistent patterns of spatial covariation. This is a desirable effect, because being able to calculate phase at each time-point and contact does not imply the phase is meaningful. By choosing a cut-off in the reduced-rank SVD model that excludes negligible sources of explanatory variance, the meaningless measurements of phase are excluded. This would not be the case if spatial frequency estimates were made directly on the phase estimates.

SVD has been described as a generalization of Fourier analysis (96), in that it is used to empirically recover a set of bases that allow the original data to be represented as weighted sums of that basis. We illustrate in Figure 3 that analysis of suitable signals using SVD enables recovery of a set of bases that are approximations of sinusoidal basis functions, (see also (46) for a critical view). ‘Suitable’ here means (a) monotonically decreasing in power with frequency or calculated at a clear peak in the power spectrum, (b) containing no other strongly competing explanatory features (such as gradients) and, importantly, (c) smoothness of values in the underlying function. Suitable signals include digitized classical music, noise generated by random walk processes (proof in Shinn 2023) and intra- and extra-cranial cortical measurements (8,79). This specific ability of SVD to recover approximately sinusoidal bases means that, though empirically driven, in the present research there is a close theoretical relation between SVD and Fourier analysis.

The reasons that SVD finds a sinusoidal basis are treated in Shinn (2023), for the case of a uniform measurement array. Essentially, spatial smoothness in the signal means neighbouring values quantify the curvature in the form a central difference approximation of the second derivative. In the present data, the spatial smoothness is a function of real correlations in cortical activity across space, plus short-range volume conduction (<1cm; 91). The continuous function corresponding to this difference approximation is the harmonic oscillator, which is satisfied by linear combinations of sines and cosines. Depending on the boundary conditions of the signal, the SVD components may also involve the first order forward finite difference approximation of the first derivative. Here also, the corresponding function in the continuous case is also a sinusoid. Also depending upon the boundary conditions, the LSVs may not be pure Fourier components but weighted sums of the terms in the Fourier series (46).

In the present research, the vectors of phase (over the array) were decomposed into uncorrelated components using SVD. The matrix of complex-valued phases, 𝚽_𝑠𝑡, is submitted to the SVD. Two outputs from the SVD are used to compute the SF spectrum. These are the LSVs, which are each complex-valued vectors with indexes s. Each value within an LSV, 𝛇_𝑠, consisted of a phase, 𝛚_𝑠 = ∠𝛇_𝑠, and a magnitude, 𝒂_𝑠 = |𝛇_𝑠|, which could vary over the LSV, e.g. 𝒂_𝑠 tend to be low at singularities (rapid rates of change) in 𝛚_𝑠. The second output, the LSV scaling value, 𝛓_𝑟, is obtained from the main-diagonal values of the scaling matrix 𝚺_𝑟, and comprises the ordered contributions, or weightings, of each LSV’s contribution to the original data.

Subsequent to the SVD, it is crucial to analyze the nature of the components, rather than just assume they are physically meaningful decompositions of the measurements (98). This reality check arises the same sense for Fourier analysis, which is guaranteed to decompose any continuous signal into sums of sinusoids, even if the underlying physical process is not sinusoidal, such as in heat diffusion (99). Here we use the SVD as a preliminary step to then estimate the spatial frequency spectrum in sEEG, rather than make specific claims about the empirical meaningfulness of the components, although related methods and claims have been used in a wide range of neuroscience research (15,89–93).

The threshold for exclusion of ‘noise’ components is justified by the wide range of dimension reductions (between 10 and 40 components) that were trialed. The method does not depend on choosing a particular number of signal components. The noisier the component (from the point of view of spatially coherent patterns of phase), the lower its contribution to variance explained by the SVD (and hence to the SF calculation). In the results reported here we used the first 14 LSVs and singular values.

Here, we describe in more detail the relation between Fourier analysis and SVD methods to introduce the logic of the method. We limit our description to the one-dimensional case, without loss of generality. For suitable signals, as defined earlier in this section, each of the first few LSVs, 𝑣_𝑘, is an approximate complex sinusoid with wavenumber equal to k.

This means the SVD has decomposed the matrix of complex-valued phase into bases that approximate a Fourier decomposition in the spatial dimension

where Σn is the sum over the spatial array of length N samples. Note that this means the maximum spatial frequency that can be measured is equivalent in size to a one cycle sinusoid (i.e., when k=1) over the sampling array.

Just as the LSVs play the same role as the last term in equation 4, 𝑒 so the singular values describe the weights of each sinusoidal term, corresponding to the sums over 𝜓_𝑛 in Fourier analysis. The spatial frequency spectrum can therefore be estimated after the SVD by summing over the singular values assigned to unique (or by binning over a limited range of) spatial frequencies. This procedure is illustrated in Figures 3 and 4.

This correspondence, between Fourier analysis and SVD of spatially organized phase, allows spatial frequency spectra to be estimated on semi-regular arrays, such as EEG, MEG as well as ECoG arrays on the surface of the cortex (8,11,29). The first dozen or so LSVs can be assigned a spatial frequency by calculation of the circular range of phase across the array. This assignment is possible because the maps of phase so-produced are smooth and regular, even if the measurement array is semi-regular (11,29). Note that the smooth maps of phase arise from the regularities (covariances) in the phase data itself, and do not require information about the spatial locations of the measurement sites. The measurement locations are only required, post hoc, to estimate the spatial frequency of the maps, in cycles/m.

We hypothesized from our previous observation regarding LSV of phase (4,11,29) that the methods would generalize sufficiently well to estimate SF spectra in the case of highly irregular measurement arrays. Ruiz-Antolin & Townsend (41) show that Fourier analysis can be performed on non-uniform arrays, provided each (irregularly placed) sampling site can be re- mapped to a sampling site on a uniform array. The key observation to their procedure is that the effects of the re-mapping can be well approximated as a low rank matrix, and these effects can be computed by a limited number of additional, corrective, Fourier analyses. It seems reasonable then to assume that, just as SVD captures sinusoidal components on (quasi-) regular arrays, SVD can also capture components on highly irregular arrays that are combinations of sinusoidal and the (low-rank) effects of sampling site-displacement. Such a mathematical proof is beyond the scope of the present research. Instead, we rely on surrogate testing to show that the present method accuracy assesses the shape of the SF spectra. We discuss this issue further at the end of the Methods section.

Removal of DC component of phase spectrum

Alexander et al. (2016) introduced a quantitative method for splitting spatio-temporal waves into positive and negative components. This enables the pure TW component to be measured, while removing the effects of standing wave activity, including global phase-locked synchrony. We apply this method to the LSVs of phase to account for the effects of average reference, which introduce an artifactual global synchrony, thereby inflating the lowest SF’s power. Since we are considering the spatial spectra of phase, the removal of globally synchronous components has the effect of removing the DC component of the spectrum. In general, this component is outside the range of SFs of interest here (8 to 50 c/m), and inclusion or exclusion of this methodological step did not influence the present findings. It allowed us never-the-less to discount the DC component in the phase spectrum.

The positive SF component is calculated as follows. The procedure takes vector of Fourier components (phase-only or including the amplitude) over the array. Specifically, the LSVs are used as input to this algorithm. The LSV, which has complex magnitude and phase, 𝑟𝜙_𝑠, is first extended to one full cycle of phase

where × is the Cartesian product, 𝜙_𝑐 is a vector of unit length phases over n samples, where

and 𝑓_𝑠 is the sampling frequency and f is the centre frequency of the phase estimate. We then perform an SVD on the real part of the temporally extended phase vector. The pure TW can then be estimated combining the first and second LSVs into the real and imaginary parts of a complex-valued phase vector

where u are the LSVs of ℝ(𝑟𝜙_𝑠𝑐). Examples of the LSVs of 𝚽_𝑠𝑡 for the MEG data, before and after the removal of standing wave effects, are shown in Supplementary Figure 9.

Multi-scale differencing of phase

For the sEEG, the extreme irregularity of the array means the smoothness of the phase in the LSVs of 𝚽_𝑠𝑡 breaks down. To overcome this limit, we introduce further numerical analyses based on multi-scale differencing. We used multi-scale difference methods since we were interested in estimating the entire spatial frequency spectrum, including global ranges. The accuracy of the multi-scale differencing method was compared to SF spectra obtained using an alternative measurement modality (MEG) and using surrogate data whose ground truth was known by construction. Surrogate testing indicated that use of nearest neighbour differences alone was not sensitive to lower SFs known to be present in the data. Further, we demonstrate an important effect of the distortions introduced by the irregularity of sEEG contacts to smooth phase maps calculated via SVD. These distortions meant the largest inter- contact distances were most sensitive to low SF parts of the signal and the smallest contact distances were most sensitive to high SFs. This observation differs from smooth wave maps (4,29), such as those obtained from MEG (Figure 3F). Smooth wave maps provide the same estimates of rate of change of phase per metre at all spatial scales of differencing, because of their uniform rates of change in phase with distance.

The SF contributions of each LSV were quantified by considering of triplets of phase values, in the form of approximately equilateral triangles mapped onto the folded cortical surface. The rate of change of phase per metre, across triangles of contact sin each LSV, was used for each estimate of spatial frequency. The results therefore did not directly require the inference of spatial oscillatory components (98).

We group contacts into approximately equilateral triangles so the direction of phase flow could be unambiguously determined. Comparing phase between pairs of contacts, such as in phase coherence (18,100), along a line of contacts (3,26,27), or for a highly elongated triangle, only measures the phase flow along the axis in which the contacts are aligned. This procedure can greatly overestimate the rate of change of phase, in the case where the flow has a strong orthogonal component, i.e., is directed across the axis on which the contacts are arranged. We therefore re-organize the irregular field of sensors into triplets of approximately equilateral triangles.

We collated the unique triangles, i.e., contacts could appear in multiple triangles. First, the contact’s position was projected to the nearest point on the surface of the cortex. We then computed the geodesic distance (47) between every pair of contacts using the surface mesh derived from the participant’s MRI scan. Triangles were counted as approximately equilateral if the most acute angle was greater than π/4. Calculating geodesic distances enabled the true distances along the surface of the cortex to be accounted for and, since the TWs move through the cortical sheet, this is presumed more accurate than taking the Cartesian distances from the folded cortex. Use of the geodesic distance measure—which is only defined within a hemisphere—also meant the triangles were only defined within hemispheres, i.e., a triplet of contacts could not span both hemispheres. Contacts in the sEEG participants often exist in both hemispheres. However, the SVD step allowed patterns of phase to be detected from the regularities (covariances) across hemispheres. Examples of geodesically defined triangles are shown in Figure 5 and Supplementary Figure 1, from both large and small size triangle bins.

The number of triangles varied from 11051 down to 66 across participants, resampled from a range of grey-matter contacts 108 down to 31. Triangle spatial frequency, 𝜉_𝑇, is given by one over the triangle linear size (square-root of the area of the triangle):

𝐀̅̅𝐁̅̅ and 𝐂̅̅𝐀̅̅ are edge vectors of the triangle with vertices A, B, C and × is the cross-product. A is nominally at coordinate (0,0) within the triangle, and B and C are computed via the relevant geodesic distances across the cortical sheet. The spatial frequency distribution of triangle sizes formed a peaked distribution with maximum triangle counts occurring near the large end of triangle size range (see Supplementary Figure 1). An effective uniform distribution of 𝜉_𝑇 was achieved by dividing the range of 𝜉_𝑇 into 32 equally spaced spatial frequencies, corresponding to triangle sizes from 32cm to 1cm. An additional bin with the smallest triangle sizes was allowed to include the rare triangles of size <1cm and, depending on the participant, the lowest 𝜉_𝑇 bins could be empty if the contacts spanned a more limited range than 16+cm. A normalization term, 𝑁_𝑇 was then applied, such that the weight of each triangle’s contribution to the SF calculation was inversely proportional to the number of triangles in that triangle’s bin. This meant the SF estimates of phase were undertaken using an underlying set of distances with an effective uniform distribution in 𝜉_𝑇.

For each LSV, and for each triangle of three contacts, we computed the rate of change in phase per metre on the triangle in the abstracted plane implied by the three inter-contact geodesic distances. Hence the spatial frequency, 𝜉, is given by:

where ∇_𝐻 is the horizontal (i.e., two dimensional) del-operator, so

and 𝛚 = ∠𝛇. Each vector in eqs. 9 and 10 has an implicit triangle index, T, which indexes a triple of contacts drawn from s, and LSV index, r, so 𝜉 can be rewritten as 𝜉_𝑟𝑇.

The weighting of each triangle’s contribution to the spatial frequency spectrum is computed as a function of the scaling value of the r^th LSV (i.e. the amount of variance explained), the magnitude of the LSV’s complex-value at the triple of relevant contacts and the normalization term for each triangle, and is given by:

where 〈𝒂〉_𝑇 is the mean element-wise magnitude of 𝛇_𝑠 for the subset of three vertices in s for LSV r. Combining the spatial frequency estimate for each geodesic triangle in each LSV (eq. 9), with the weighting of each triangle’s contribution to the power (eq. 11), allows the construction of the estimated spatial frequency spectrum. For K SF bins, we have a bin width, w:

with each i^th weight 𝐖_𝑟𝑇 being assigned to the bin corresponding to its the computed i^th SF, 𝜉_𝑟𝑇. The estimated power in arbitrary units (A.U.) for the k^th bin, 𝑘 ∈ [1, 𝐾] is then,

We used 32 SF bins to plot the SF spectra, after choosing min(𝜉) and max(𝜉) on a reasonable range (e.g. SFs equivalent to twice the maximum triangle size through to 50c/m for plotting the sEEG data; SFs equivalent to the maximum triangle range through to 50c/m for computing regression lines). For some additional analyses, e.g. the effects of average referencing on the spectrum, the lowest SF plotted was 0.5c/m. The estimated power is in arbitrary units because it is a function of the number of contacts submitted to the SVD and number of components retained. Since these are both constant for the calculation of a given spectra, the relative power and different SFs is the relevant quantity for each participant.

Surrogate testing

The primary method of surrogate testing introduced an artificial SF component to the sEEG signal, whose ground truth was known. The artificial SF components were created by inflating the cortical mesh of each participant to a sphere, then projecting a phase gradient onto the sphere, in either the anterior-posterior or left-right directions. The phase gradients were varied in wavelength from 2 to 16 c/m, and the positions of the actual contacts projected onto the inflated sphere used to infer surrogate basis vectors. Surrogate right singular vectors (i.e. time- series of weightings for the LSVs) were constructed by reversing the time-order of the empirical right singular vectors for that participant at the relevant TF, then taking the complex conjugate to re-establish the correct direction of phase change over time. A surrogate time series in complex-valued phase was constructed using each surrogate basis vector as LSV, a scaling value of 1.0 and a surrogate right singular vector, according to eq. 2. The surrogate time-series in phase was added to the empirical phase signal at each TF, with the surrogate signal weighted, as 𝑤_𝑠𝑢𝑟, across a range of values to both (1) perturb the empirical signal (1/3 surrogate to 2/3 signal ratio), and (2) to see the effects of perturbation by realistic signals on the veridical surrogate phases (7/8 surrogate to 1/8 signal ratio).

where Φ_𝜉 is the phase surrogate signal at some SF and direction and Φ_𝑠𝑢𝑟is the combined signal.

We additionally analyzed two participants from a different, MEG, data-set (101), as a source of neurophysiologically realistic time-series with known SF spectra (Alexander et al., 2016). Surrogate testing of the method using MEG was included because the large number of quasi- regularly spaced contacts enables recovery of smooth wave maps using SVD, with visually apparent SF values for each map. We then used the present methods estimate the distinct frequency spectra for each individual LSV. In addition, we constructed MEG surrogates using the first 28 LSVs of phase, but reversing the relative weighting of the singular values, i.e. 𝜍_𝑖, in eq. 2. This produced an MEG surrogate signal, in which small spatial scales of activity dominated the signal, called hereafter the f spectrum (as opposed to the usual 1/f spectrum). We additionally made use of the large number of contacts in the MEG array to test the effects of sparsity of measurement contacts. The phase data for the f spectrum surrogate, along with the unmodified MEG phase data, were further evaluated by randomly removing subsets of contacts until only 1/3 remained (51 of 151). This analysis took advantage of the large number of sensors available with MEG, so that highly irregular arrays are constructed by increasing the sparsity. This analysis was to rule-out that the estimated SF spectra were artefactually created by the sparsity of the sEEG arrays, when analyzed using our method. We also used these MEG time-series to test the effects of an average reference on our method, since with MEG the case of no reference is also available to analyze.

Effects of TF power

The main analyses here were undertaken using complex-valued phase but generalize naturally to complex values that are not unit length e.g. by not excluding TF power. Each value computed by the Morlet wavelet is denoted 𝑟𝜙 = 𝑟𝑒^−𝑖𝜃 where r is the power of the signal. The matrix 𝐑𝚽_𝑠𝑡 was submitted to the SVD, and the SF spectra estimated. There is inherent noise in the TF power, mostly due to the accuracy and residual uncertainty of contact placement relative to the mm scales of the cortical sheet. This results in large variation in the mean amplitude signal across contacts. For this reason, we also computed the normalized power at each contact, 𝑟̅, such that the normalized power at a given contact varied between zero and one over the recording, and 𝐑̅𝚽_𝑠𝑡 was submitted to the SVD.

Statistical analysis

Statistical testing of linear regressions reported in the results was undertaken as correlation between the empirical values (𝑦) and the regression model values (𝑦^). For regressions of SF versus estimated power, linear fits of the model, 𝝆, were undertaken using a log-log scaling of wavelength (1.0/SF) with power. This combination produced the best linear fits. The statistical significance of these correlations, over multiple comparisons, was tested by permutation methods. Each permutation of the data comprised a correlation between 𝑦s in the regression, and random shuffling of 𝑦^s (10,000 permutations). The probability of arriving at a mean correlation (ignoring the sign of the mean), 〈𝝆〉_𝑛, over n frequencies, participants or other analysis dimension was compared to the actual (unsigned magnitude of the) mean correlation (i.e., a two-tailed test). The tests for the significance of a family of regressions are reported as “variable A versus variable B (permutation test over n-condition, p-level)”. The p-level was set to p<0.05 for planned analyses (individual participant regressions of wavelength versus power, over all temporal frequencies) and p<0.001 for others. When n was one, i.e., only one correlation rather than the mean of a group of correlations, this is indicated by omission of the n-condition term.

Effects of triangle size on SF estimates

We have assumed that the LSV derived from irregular arrays can be understood as Fourier components along with low rank ‘corrections’ for contact displacements (41). As part of our methodological approach to apply multi-scale difference estimates of SF, we examined peak SF as a function of triangle size. Examination of SF spectra revealed that triangle bin size and SF peaks derived from each bin had a strong linear relationship. That is, small triangles were more sensitive to high SFs, large triangles were more sensitive to low SFs. Statistical analysis of spectra produced from each separate triangle size bin showed that triangle size (as mean triangle size over a bin) linearly increased with the peak in the SF spectra for that bin size (as logarithm of wavelength; permutation test over participants, p<0.001). Examples of this effect are shown in Supplementary Figure 10. This effect is not expected under standard Fourier analysis on regular grids, for which the SF estimates will not vary if calculated by difference methods at alternative values for delta-distance. Likewise, the SF estimates from standard Fourier analysis will not substantively differ at most (except highest) SFs if the data are sampled with different underlying but uniform contact spacings. This effect arises in the present sEEG data due the irregular spacing of the contacts, which become encoded in the LSVs of the SVD.

We analyzed this effect with surrogate methods, manipulating the inputs to the SVD of the random walk surrogate shown in Figure 2A. We performed the same SVD while systematically changing the distribution of intercontact distances in the measurement array. The results are shown in Supplementary Figure 11. A sparse, uniform array, covering the same spatial range as our first analysis, retrieved the same SF components (panel F, column 1). A dense array with a narrow range of contacts retrieved only high SF components (panel F, column 5), relatively speaking, but exactly the components defined by the maximum array size. These results are expected from considerations from standard Fourier analysis. Intermediate cases of grid spacing were produced via modulo rearrangements of the contacts (102). The maximum contact range and the Nyquist frequency was varied systematically between the two extremes. Careful visual inspection indicates that SVD produced LSVs with high SF components that were a function of the minimum intercontact distance; this is shown by the rate of change of values in the LSVs in runs of nearest neighbour contacts (panel F, columns 2-4). However, lower SFs representations were also present in the extracted components; this is demonstrated by interpolating between equivalent points in each run of nearest neighbour contacts (panel G).

Consider the application of a simplified, one-dimensional version of the present methods i.e. estimating SF spectra from the LSVs via multi-scale contact distances rather than triangle sizes. When using a multi-scale contact array, the highest SFs can be extracted from smallest neighbouring intercontact distances. The lowest SFs can be extracted from the largest intercontact distances. This effect arises from the way the SVD reorganizes the Fourier information under conditions of non-uniform sampling (cf. Ruiz-Antolin & Townsend, 2018). This is the same effect we see for triangle size in the sEEG data, i.e. small triangles are sensitive to high SFs and large triangles to low. We therefore used multi-scale differencing to estimate the SF spectra, in order to capture the full range of SF information present in each LSV. The accuracy of this approach was then confirmed by the surrogate testing of the methods.

Code availability

The key Python routines used in this research can be found here: https://github.com/DugueLab/Traveling-wave-analysis/blob/main/estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb

Supplementary figures

Supplementary Figure 1.

Supplementary Figure 2.

Supplementary Figure 3.

Supplementary Figure 4.

Supplementary Figure 5.

Supplementary Figure 6.

Supplementary Figure 7.

Supplementary Figure 8.

Supplementary Figure 9.

Supplementary Figure 10.

Supplementary Figure 11.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 852139 – Laura Dugué).

We want to thank to Dr. R. Austin Benn for helpful advice on measuring geodesic distances in the cortex, and to Dr. Kirsten Petras for reading versions of the manuscript and providing useful feedback.

Additional files

Supplementary video 1. Dorsal view of one participant’s dominant spatial patterns of phase at 3.5Hz. The video shows the time-series of a low rank model of phase dynamics, constructed from the first three LSVs. Together, these LSVs account for 43.6% of the variance in phase. This low SF activity dominates the dynamics. The entire video comprises 6s of cortical activity. Colour values are cosine of the phase (yellow=-1, red=1). Values between contacts are interpolated for visualization purposes. Anterior region of cortex is to the left.

Supplementary video 2. Anterior view of one participant’s dominant spatial patterns of phase at 8.0Hz. The video shows the time-series of a low rank model of phase dynamics, constructed from the first four LSVs. Together, these LSVs account for 39.3% of the variance in phase. This low SF activity dominates the dynamics. The entire video comprises 6s of cortical activity. Colour values are cosine of the phase (yellow=-1, red=1). Values between contacts are interpolated for visualization purposes. Left hemisphere of cortex is to the left.

Supplementary video 3. Left view of one participant’s dominant spatial patterns of phase at 8.0Hz. The video shows the time-series of a low rank model of phase dynamics, constructed from the first three LSVs. Together, these LSVs account for 40.8% of the variance in phase. This low SF activity dominates the dynamics. The entire video comprises 6s of cortical activity. Colour values are cosine of the phase (yellow=-1, red=1). Values between contacts are interpolated for visualization purposes. Anterior region of cortex is to the left.