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). Macroscopic TWs, specifically, are monotonic phase gradients measured over more than 15cm, 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 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,17,18). Mathematical analyses indicate activity modes that are associated with damped TWs (14), and under appropriate boundary conditions, these correspond to macroscopic 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,19), 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 macroscopic TWs measured extra-cranially. We highlight three prior studies (20,21,24) here due their 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 (20,21), 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 (25–27). 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 modeling, 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 measurable Fourier spectra, and these spectra give distinct predictions. While the expected peak in the spatial frequency spectrum for the veridical TWs hypothesis is at low spatial frequency, the alternative hypotheses for artefactual macroscopic TWs show peaks at high spatial frequencies. In order to test these hypotheses, we therefore estimate the relative power of different bands of spatial frequencies (f_s), comparing high (f_s>6 cycles/m) to low (3<f_s<6 cycles/m) spatial frequency contributions to the measurable phase; these spatial frequencies are approximately equivalent to wavelengths (λ = 1 / f_s) of λ<15 cm vs. 15<λ<35 cm.

The alternative explanations that have been proposed for measurements of extra-cranial macroscopic TWs rely on two assumptions. First, in order 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 (28). 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 result in activity with high variation across aggregated samples being made quiescent, thereby treating this activity as noise (28).

In the present study, spatial frequency estimates are made at the single sample level rather than for sample aggregates (e.g., trial means), up to the constraints imposed by estimating phase over a short temporal window and lossless decomposition of a participant’s phase data using singular value decomposition (SVD).

Research on TWs has consistently revealed that the spatial power of phase decreases monotonically with spatial frequency. Intracranial measurements show that the spatial frequency spectrum of phase declines with power from 5 cycles/m (10cm wavelengths) (2,7,8,28–30). The estimates, however, use cortical surface arrays with limited spatial coverage (<10cm) and were therefore insensitive to wavelengths longer than 10cm (see Figure 1). Likewise, animal studies also point to a spatial frequency spectrum for phase that monotonically decreases in power with increasing spatial frequency (6), but at measurement scales of 7mm or smaller (6,31–33). Estimates of the relative contribution of macroscopic TW to the spatial frequency spectrum of cortical phase dynamics have used EEG (7,29) (see Figure 1) or magnetoencephalogram (MEG) (11), and consistently describe a spatial frequency 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 spatial frequency 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) (34). Phase dynamics at the largest cortical scales 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,28,35) (see Figure 1). A low pass filter will suppress higher spatial frequencies, thereby inflating the relative effect of the low spatial frequencies quantified from the phase dynamics and introducing an ambiguity. In the present research we use cortical depth electrodes in order to assess that part of the spatial frequency spectrum of phase between 3 and 20 cycles/m. Stereotactic EEG contacts have excellent spatial resolution and coverage (although often-times sparse and irregular) which can extend to 25cm of the cortical sheet. By extending grey-matter measurements of the spatial frequency spectrum to long wavelengths, we are able to accurately assess the relative contribution macroscopic phase dynamics to cortical activity.

Results

Despite the published research on spatial frequency spectra of cortical phase dynamics (6,7,29,36), decades of neuroscience research has favored explanations of extra-cranial measurements in terms of mesoscopic TWs or cortical sources mixing at the sensor level (20,21,24). Here, we are able to 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 at the single trial level within the grey matter. As a consequence, the proposed mechanisms of localized TWs, oscillating sources or large-scale phase gradients must result in a spatial frequency spectrum of phase measurements that is consistent 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 or local inflection) in the spatial frequency spectrum. These high spatial frequency components have been hypothesized to combine, via blurring of signal, to form apparent low spatial frequency macroscopic TWs at the sensor level (20,21,24). This effect will therefore only be apparent if other cortical activity is relatively quiescent; the high spatial frequency components must really be a prominent part of the signal. Alternatively, single trial level measurements within the grey matter that reveal a dominance of low spatial frequency components in the spectrum of phase, at scales up to macroscopic TWs, cannot support the blurring or source mixing explanation for macroscopic TWs. This is because the sensor level measurements will be dominated by the real low spatial frequency component of the cortical phase dynamics. Measurement of macroscopic TWs in EEG or MEG can only create a purely artefactual low spatial frequency gradient of phase, through mixing and blurring, if no such gradient dominates the underlying signal in the first place (compare Figure 1A & 1D).

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 11 participants. Note that most subsequent procedures (construction of spatial frequency spectra, statistical testing of link to measurement array size and temporal frequency) were carried out at the single participant level, so each participant comprises a test of the hypothesis under scrutiny. The stereotactic EEG contacts in the grey matter were re-referenced using the average of low-amplitude white matter contacts (37). A typical participant’s contacts are shown in Figure 2D. From the stereotactic EEG time-series, complex-valued phase was analyzed at 36 frequencies, using short time-series methods (2 cycle Morlet wavelets). Centre frequencies ranged from 1.0 to 128.0Hz, on a logarithmic scale.

Figure 2.

The highly irregular sampling grid provided by the stereotactic EEG limits the statements that can be made about the cortical phase dynamics. Selecting model components peaked in either the macroscopic or mesoscopic range still results in broad spectral curves with tails well outside the selected range. This reflects the inherent noise of the spatial array. For this reason, it was not possible to assess the smoothness of the cortical phase dynamics in the present research. We therefore consistently refer to the spatial frequency 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,28).

The numerical methods described in the present research are novel. They rely on the use of complex-valued phase (8) to avoid the (error) noise introduced by unwrapping phase when using real-valued techniques (c.f. 2,4,26,36,37). The irregular sampling of cortical spatial coordinates via stereotactic EEG was partly overcome by the resampling of the phase data into triplets corresponding to the vertices of approximately equilateral triangles within the cortical sheet. Between-contact distances were calculated as geodesic distances (40) in the cortical sheet to accurately reflect the underlying cortical topology. We used higher-order SVD, which is a generalization of SVD to greater than two dimensions, i.e., for use with tensors. This enabled the tensor of samples × triangles × vertices to be distilled into the dominant phase patterns and then characterize flows of phase across a limited set of uncorrelated model triangles. These techniques enabled the spatial frequency spectra of the cortical phase dynamics to be quantified for each temporal frequency.

The strategy of this research was to use linear algebra techniques to empirically construct a Fourier-like basis with which to decompose the phase signal, following and extending previous techniques (4,8,11).

From this simpler set of empirical functions (41) we are able to estimate the spatial frequency spectra of the cortical phase dynamics. The general strategy for empirical construction of a Fourier-like basis is illustrated in Figure 2. In this simplified case, a one-dimensional measurement array of real-valued waves (Figure 2A) is first converted to a time × space matrix of unit-length, complex-valued phase using Morlet wavelets tuned to the dominant temporal frequency. The matrix of phase is decomposed into singular vectors, each consisting of an approximately sinusoidal function in the spatial dimension (Figure 2B). The singular vectors enumerate a complex basis analogous to Fourier analysis. Summing these basis vectors with appropriate weights reproduces the original wave data. Along the way, the SVD quantifies the amount of variance attributable to each singular vector. These singular values can be used as proxy for spectral power (Figure 2C).

We have previously shown that SVD can be used to decompose phase data obtained via EEG, MEG or ECoG (4,28). 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. The first three singular vectors extracted using SVD were approximately equal in wavelength to the size of the measurement array. Since the singular vectors were smooth gradients of phase, this meant these first three components were macroscopic 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 singular vectors were approximately wavenumber equal to two (i.e., twice the spatial frequency), and each successive three singular vectors increasing in wavenumber by one (28).

We applied the SVD methods for spatial frequency spectrum estimation to the stereotactic EEG data. The methods and results for this analysis are described in Supplementary Section A. The results (see Supplementary Figure 2) support the conclusions presented here, i.e., low spatial frequency cortical phase dynamics dominate over higher spatial frequency components. However, the singular vectors did not produce smooth maps of phase, since they were strongly affected by the noise of the irregular stereotactic EEG arrays. The maps were therefore difficult to tease apart for further verification. In particular, we required a version of the SVD method that allowed different wavelengths of signal to be separated and verified, and way to assess the relation between spatial power and spatial frequency in a fine-grained manner.

Here we present a slightly different entry-point to the problem, more amenable to verification of the elements of the decomposition. In short, we define sets of triangles on the cortical sheet, with the stereotactic EEG contacts forming the vertices of each triangle. An example participant’s triangles are shown in Figure 2F. All subsequent single participant examples reported here-in are from this same participant. The triangles are used to create a larger data structure from the phase values at each temporal frequency, of dimensions samples × triangles × vertices. The new data structure is submitted to a higher-order SVD (42), so as to characterize the phase flow across triangles of different scales.

At each temporal frequency, we used higher-order SVD to reduce the dimensionality of the phase along both the samples dimension and the triangles dimension. The results reported here were not dependent on the number of reduced dimensions. We trialed 5 reduced-samples, 5 reduced-triangles and successively doubled the numbers until 40 reduced-samples, 40 reduced-triangles. Choice of reduced dimensions, within this range, did not alter the substantive results. The remainder of the results reported here are for the case of 20 reduced-samples, 20 reduced-triangles, giving 400 model components in total (see Equation 12). The numerical pipeline for the analysis is illustrated schematically in Figure 2G-L.

Each model component captured a characteristic triplet of phase relations as phase values on the vertices of the reduced-triangles, giving a direction (relative to the three vertices) and magnitude of phase change. Examples of model components from a single participant are shown in Figure 3, which also shows the distributions of phase and phase-differences for the contact triangle which had the highest weighting on the given model component.

Figure 3.

We first quantified the estimated spatial frequency spectrum for a single participant, across all temporal frequencies (Figure 4). Lower temporal frequencies had slightly more spatial power than higher temporal frequencies (Panel A). However, this is likely due 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 therefore not a feature of higher temporal frequency signal per se. Panel B shows the characteristic shape of the estimated spatial frequency spectrum, peaking in power at a wavelength slightly larger than the largest inter-contact distance on the measurement array. The spatial power of the phase dynamics decreased monotonically with increasing spatial frequency. Within the range beginning at the maximum inter-contact distance in the measurement array (equivalent to ∼4.8 cycles/m in this case), the spatial power drops steeply, before flattening at higher spatial frequencies. This was true of all participants (see Supplementary Figure 4). The cross-participant aggregate spectrum is shown in Panel C, confirming the observed distribution. 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,29,36). The present result extends the previously measured spectrum into scales equivalent to a large fraction of the cortical size, as well as to temporal frequencies up to 100Hz. In order to check the accuracy of the spatial frequency estimates, we selected components that had peaks in power spectrum in the macroscopic range (2-3 cycles/m), or alternatively, in the mesoscopic range (6-50 cycles/m). Figure 5 shows single sample phase estimates that had a good fit (i.e., high ρ, see methods) to either the macroscopic or mesoscopic model components at a temporal frequency of 8Hz. Supplementary Figure 5 shows samples taken from the same participant, but across a wide range of temporal frequencies. Samples with good fits to the macroscopic components had approximately three times the wavelength compared with samples with good fits to the mesoscopic components. These samples also illustrate the limits imposed by the noise of the uneven contact distances in the sampling of the cortex. The plots do not show the smooth phase gradients obtained from more regular measurement arrays, e.g., ECoG and EEG (8,28).

Figure 4.

Figure 5.

A second analysis was carried out using the decomposition method on subsets of triangles to restrict the mean distances and thus sample the total measurable range. The results of this analysis are shown in Supplementary Figure 6. The peak spatial power for these spatial frequency estimates increases nonlinearly with measurement array size. Linear regressions of log₂ of power at peak wavelength versus mean contact distance were significant for all participants (Wald test, p<0.05), and showed that power increased greater than linearly with increasing peak wavelength. Given the noise inherent in the contact placements, and the complex interactions of various sources of noise with the higher order SVD, we did not here attempt to estimate the exact function corresponding to each participant’s spectra. This result, along with the shape of the complete spatial frequency spectrum for each participant, is consistent with previous studies showing that the spectrum can be approximated by the function 1/f_s^α at all cortical distance ranges (6,7,11,29).

We quantified the relation between the power at peak spatial frequency and temporal frequency, using the participant-aggregate spatial frequency spectra derived from all contact triangles. There was a strong inverse relation between temporal frequency and power at peak spatial frequency (Wald test, p<0.05; Figure 4D). However, this effect, while neatly linear at the participant aggregate level, was for a truncated range. In other words, the changes in spatial power at peak with temporal frequency were small compared to the changes in power over the entire spatial frequency range. Even at the highest temporal frequencies, the spatial frequency spectrum followed a monotonically decreasing curve from the wavelength equivalent to the size of the array down to the smallest measured wavelengths. It might also be that the lower spatial frequency at peak spatial power for high temporal frequencies was at least in part a function of higher noise in the estimation of phase at these frequencies, i.e., fewer numbers of samples in the 2-cycle Morlet wavelet.

Additionally, low temporal frequencies peaked in spatial power at slightly higher spatial frequencies than did high temporal frequencies (Wald test, p<0.05; Figure 4E, participant-aggregate spatial frequency spectrum). However, these effects were small compared to the nonlinear decline in spatial power with spatial frequency. The weakness of the former relation illustrates the second major finding of this research: the shape of the spatial frequency spectra was similar across all temporal frequencies tested. The participant’s spectra shown in Figure 4A is typical in this respect. This result was surprising, as it implies that fast activity well into the gamma range is strongly coordinated in the cortex at the largest spatial scales, since long wavelength activity also dominates these temporal frequencies.

Considerations from Fourier analysis, outlined in the Introduction, show that the maximum detectable spatial frequency depends strongly on the size of the measurement array. Since the array sizes used here did not reach the full extent of the cortex, we were at the outset interested in estimating the peak spatial frequency for the case where the array size extended to the entire cortex. The results show that the peak spatial frequency across participants (regardless of temporal frequency) was dependent on the maximum size of the contact array. We first regressed the peak wavelength, for each temporal frequency, against all participants’ maximum array size (Figure 4F). The slopes of the regression lines ranged from -0.23 to 1.5, and all but one frequency in the delta range were positive. For visual convenience, the regression lines are extended to larger cortical distances. These slopes are strongly affected by the variation in participant’ ss contact placements, which here acts as a source of noise. In particular, the maximum inter-contact distance for a participant only weakly reflects their distribution of inter-contact distances.

We therefore turned to the analysis using higher order SVD on subsets of triangles, in order to partition the triangle size ranges and so manipulate mean triangle size (see Supplementary Figure 6F). The regressions for peak wavelength (calculated after aggregating spectra over all temporal frequencies) versus mean contact distance were significant for 9 of 10 participants (Wald test, p<0.05), the exception being one participant with only 5 triangle size bins, and thus statistically underpowered. Using this artificial constriction of triangle sizes, we show that the peak wavelength estimate obtained for each subset is strongly constrained by the typical inter-contact distances. It increases linearly, and for visual convenience, the figure shows the regression lines extended beyond the maximum sized subset of triangles. In this regression analysis of subsets of triangles, the slopes of the regression lines were in the range 1.9 to 4.1 across participants. These slopes (>1) reflect two effects. First, the mean distance is given as the distance measure, so the triangles included in each analysis contained contact distances up to approximately 1.5x this mean distance. Second, while the maximum wavelength detectable is strongly constrained by inter-contact distances, the model phase triplets were optimally tuned to wavelengths slightly larger than the triangle’s linear size. For this reason, the peak spatial frequency for each participant’s subset analysis was larger than the maximum inter-contact distance. This effect can also be seen for the decomposition of 1-dimensional model data provided in the methods section, and in each participant’s (all triangle) spatial spectra (Supplementary Figure 4 and Figure 4B).

Finally, we did not consider the method of extending regression lines to large cortical distances to be accurate enough to give a concrete answer to the question of the expected peak spatial frequency when measured over the entire cortex, again due to the inherent limitations imposed by the spatial noise of participant’s contact placements. Even so, the good fits of the linear regressions suggest that the trend increases linearly with further increases in array size. This conclusion 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,29), 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 stereotactic EEG participants from all available databases. We expect that participants with greater than 30cm coverage will extend the linear trend reported here.

Previous studies of cortical phase dynamics, and specifically, TWs at a variety of spatial scales, have illustrated their functional significance (2,3,8,30,43–50). The role of the different temporal frequency bands is better understood (51) than spatial frequency, at least to the extent that different temporal frequencies are known to subserve different functions. We therefore explored how phase dynamics might differ at different spatial frequencies by examining their relationship to temporal frequency. We extracted model components with peaks in either the macroscopic or mesoscopic range, and calculated the fit of these model components to the original phase measurements, over all time samples. This procedure was repeated over the range of temporal frequencies. We then correlated the fits between each temporal frequency and all other temporal frequencies (see Figure 6). In general, the (temporal) cross-frequency correlations showed that macroscopic phase dynamics of a given frequency had lower correlations to nearby temporal frequencies. This effect is illustrated by the example participant in Figure 6A. Here, we see that the peak in (temporal) cross-frequency correlation drops off rapidly from the diagonal i.e. equal frequencies. By contrast, the off-diagonal correlations of temporal frequencies for the mesoscopic spatial frequency band were more strongly correlated with each other. It should be noted that this participant is not the most convincing example of this effect, since there are regions near the temporal frequency diagonal with negative values, but is shown here for consistency of exemplar participant. There was some variation across participants in this pattern from temporal frequency to temporal frequency likely due to the highly uneven sampling of the space of spatial frequency × cortical regions imposed by the stereotactic EEG contact placements. However, there was a general trend across participants for macroscopic scale phase dynamics to be narrower in their (temporal) cross-frequency correlations than their mesoscopic counterparts. This is shown in the participant averages of the cross-frequency correlations for macroscopic and mesoscopic dynamics and the difference in correlations between them (see Figure 6B; permutation test, p<0.001 for each frequency-by-frequency cell).

Figure 6.

The narrow peaks in cross-frequency correlations for macroscopic fit occurred despite the low frequency resolution in the estimation of phase (due to the high time resolution of the Morlet wavelets, using 2-cycle window at all temporal frequencies). It suggests that when low spatial frequency dynamics are active, they tend to be active at specific, narrow temporal frequency bands. Higher spatial frequencies, by contrast, tended be more broadband in the temporal frequency domain. This effect was strongest at the highest temporal frequencies tested (gamma), due to stronger (neigbouring) cross-frequency correlations for mesoscopic phase dynamics in the gamma range. This is true even though temporal frequency was sampled on a logarithmic scale, so that the range 32-100Hz had the same number of centre frequencies as 2-14Hz. Broadband temporal interactions for mesoscopic dynamics and narrowband interactions for macroscopic dynamics tended to occur simultaneously, and not as separate effects. This can be seen in the participant average plot of macroscopic versus mesoscopic cross-frequency correlations of fit. Regions of positive correlation form vertical bands, showing that narrowband interactions in temporal frequency for macroscopic dynamics coincided with broader temporal frequency bands for mesoscopic dynamics.

Discussion

In this research we aimed to quantify the low spatial frequency part of cortical phase dynamics spectrum, in the range of macroscopic TWs (wavelength>15cm), 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,26,27,47) or used smaller measurement arrays (<10cm) (2,6,30,32). This means either the results were biased by a low pass filter, or had a maximum measurable spatial frequency below the macroscopic 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 macroscopic TWs in the cortex (Figure 1D) and not with the alternative hypotheses for the measurement of these waves extra-cranially (Figure 1E & F). Low spatial frequency 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—local to a single grey matter contact—this implies that the majority of the variance in phase is a function of low spatial frequency phase dynamics. 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 spatial frequency phase dynamics index a distinct scale of cortical activity and, judged from the dominance of that activity, comprise an important scale for understanding cortical function (54).

The effect of spatial frequency on spatial power is greater than linear, so that power is substantially higher at low spatial frequencies than at high spatial frequencies. This is shown directly by the steep drop-off in power from the spatial frequency equivalent to the maximum array size (Figure 3A-C, Supplementary Figure 1). Indirectly, the nonlinearity appears as the logarithmic scaling of power at peak spatial frequency with mean contact triangle size. This effect of triangle size was demonstrated by grouping contact triangles into subsets of approximately equal size, sampled from the entire triangle size range (Supplementary Figure 6E & G). The shape of the spatial frequency spectrum in the range ∼3 to 20 cycles/m is therefore consistent with reports of the spectrum at smaller cortical scales (6,7,29). The cortical spatial frequency spectrum is generally described as 1/f^α, where the exponent α is a function of the type of measurement (e.g., grey matter, cortical surface or extra-cranial) (6,29,36). Since peak spatial frequency also increases with the total size of a participants’ stereotactic EEG array (maximum inter-contact range), we assume that the 1/f^α distribution of spatial frequency extends to the maximum size of the human cortex; approximately 45cm linear distance within unfolded grey matter. This would be consistent with spectral peaks at lowest spatial frequencies measured using EEG (7,29) 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 spatial frequency spectrum did not vary much with temporal frequency. This finding implies that low spatial frequency cortical phase dynamics arise at temporal frequencies up to 100Hz, the maximum temporal frequency tested here. Macroscopic TWs have previously been reported in the delta to beta ranges (3,4,26,28,45). However, TWs have also been reported for low gamma frequency (40Hz) using MEG, in an event-locked paradigm (55). Similarly, EEG-based estimation of spatial frequency has been reported for temporal frequencies from 7 to 63Hz, with similar spectra (1/f^α) for spatial frequencies over this entire range (36). We outlined in the introduction the analysis 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 (22).

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,54, cf. 35). To this we may add reports of the endogenous electro-magnetic field altering firing patterns of neurons (57,58). It seems likely that a major driver of cortical plasticity is precisely the dominant spatial frequencies found in the cortex, since post-synaptic membrane depolarization is the second major ingredient of long term potentiation (59). An evolutionary argument can also be applied, since dominant low spatial frequency 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 spatial frequencies. These considerations are supported empirically by the present results. The peak spatial frequency increases in wavelength with maximum contact distance for the participants (here, 0.16 to 0.24m). In addition, we showed that there was a tight linear relation between peak spatial frequency and average triangle size, if subsets of triangles were selected from the triangle size range for the higher-order SVD analysis. The present findings also explain why detected TWs, regardless of cortical scale, generally have a wavelength ∼2π (28,60). If the spatial spectrum is monotonically decreasing in power over the entire macro-to microscopic range, then the strongest detectable waves will be those with a wavelength approximately the size of the array.

We conclude that reports of TWs should include estimates of spatial frequency to allow or disallow comparisons between studies with different extents of cortical coverage, since they may be accessing different spatial frequency 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 are able to 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,62). This excludes these techniques from accurately estimating the spatial frequency spectra of phase dynamics since they filter out precisely the most dominant components, according to the present results.

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 mesoscopic 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 macroscopic and mesoscopic TWs are functionally significant. An interesting finding of the present research is that macroscopic phase dynamics tend to be more narrowband in the temporal frequency domain than their mesoscopic counterparts. It would be interesting to confirm this effect in more regular measurement arrays (e.g., EEG, MEG and ECoG). An intriguing question is whether the relation between cortical scale and temporal frequency specificity is related to successful prediction of future activity, i.e., large-scale phase dynamics are good at predicting future activity because of their temporal specificity. More generally, most of functionally-related results concerning macroscopic TWs have been obtained using extra-cranial methods (26,27,45,47,49). The reported linkage between late event-related potentials and TWs (4,28,43,44) 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,63,64). 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 spatial frequency, is a function of the cortical distance between the primary areas and transmodal areas (63). 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 (26,28,30,61,65). Here, we have shown that low spatial frequency 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 (66). 35 human participants were part of a multi-centre stereotactic EEG monitoring study. From the initial pool of participants, 11 had sufficient numbers of grey-matter contacts to enable further analysis (see Numerical methods for criterion). Mean age was 33.2 years, five female, three 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 2 participants, for which it was 500Hz. The number of contact ranged from 51 to 108 with a mean of 76.

The data for each participant was denoised by first calculating the standard-deviation of signal over samples, for each trial and sensor combination (37). 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 (calculated after excluding all such above threshold trial by sensor combinations).

Using the supplied automated contact labeling (checked and corrected manually by clinicians), contacts were partitioned into three groups: grey-matter, white-matter low-amplitude and other white-matter. White-matter low-amplitude were those contacts whose standard deviation of signal amplitudes were in the bottom 50% of white matter contacts. The white-matter low-amplitude contacts were used to create an average reference for the grey-matter contacts (37). All numerical procedures were trialed using notch filters to remove the 50Hz 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 is represented as a unit-length complex number, ϕ.

Phase is calculated from the stereotactic EEG time-series, for a particular temporal frequency, using short time-window Morlet wavelets (2 cycles, (67)). 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). We estimated phase for 36 centre-frequencies, logarithmically space from 1.0 to 128.0Hz. Due to the low frequency 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 c to refer to samples and contacts, respectively, giving a matrix of complex-valued phases for each temporal frequency, Φ_sc. For some participants, a lower sampling frequency at the time of acquisition meant the maximum frequency examined was 97Hz.

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 (28). 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 amplitude (28,65). 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.

SVD is a standard method to reorganize data into a more interpretable form. 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. An SVD of such a phase signal, M, then decomposes the many samples of spatial vectors of phase into a series of singular vectors, v_r in descending order of singular values, σ_r, for the r highest singular values (68).

The power of the SVD is that it re-organizes the data into components that are uncorrelated with each other (the singular vectors), but that explain known amounts of variance in the data (the singular values). The components are special vectors of space defined by the data in that they form an orthonormal basis. Weighted sums of these uncorrelated 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. 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 (69). In the present research, the components consisted of triplets of phase values, whose spatial frequency was estimated. The results therefore did not directly require the inference of spatial oscillatory components (69). In addition, the values of the estimated spatial frequency spectrum were checked against the spatial pattern of single samples of phase, in order to verify that the spatial frequency estimation procedure was correct. The exclusion of ‘noise’ components is justified by the wide range of dimension reductions (between 5 and 40 components, for both samples and triangles) that were trialed. The method does not depend on choosing a set number of signal components.

Here, we describe in more detail the relation between Fourier analysis and SVD methods to empirically derive Fourier-like basis functions. We limit our description to the one-dimensional case, without loss of generality. For appropriate signals, i.e., those with power that decreases monotonically with frequency, each of the first few singular vectors, v_k, 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 S_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 closely circumscribed by the size a one cycle sinusoid (i.e., when k=1) over the sampling array.

Just as the singular vectors play the same role as the last term in equation 4, , so the singular values describe the weights of each sinusoidal term, analogous to the sums over ψ_n. The spatial frequency spectrum can therefore be estimated using SVD by summing over the singular values assigned to each set of singular vectors with unique (or by binning over a limited range of) spatial frequencies. This procedure is illustrated in Figure 1A-C.

This analogy, 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. The first dozen or so singular vectors 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,28). 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.

The reason for grouping contacts into equilateral triangles in the present study was so the direction of phase flow could be unambiguously determined. Comparing phase between pairs of contacts, such as in phase coherence (18,70), along a line of contacts (3,25,26), 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 (40) between every pair of contacts. Triangles were counted as approximately equilateral if the geodesic distance between the closest two contacts in the triple was at least 75% of the geodesic distance between the furthest two contacts. Calculating geodesic distances enabled the distances along the surface of the cortex to be accounted for and this is more accurate than taking the Cartesian distances of the folded cortex, since the TWs move through the cortical sheet. 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. However, contact triangles could exist in both hemispheres, so the SVD procedure allowed patterns of phase to be detected from the regularities (covariances) of triangles across hemispheres.

The resulting data structure had a high degree of redundancy, requiring the use of higher-order SVD (42) in order, first, to reduce the number of cases, second, to arrive at uncorrelated estimates of the phase flow at different scales. Care was therefore also taken to equalize the size distribution of triangles such that triangle spatial frequency formed an approximately uniform distribution, where triangle spatial frequency, ξ_T, is given by one over the 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) and B and C are computed via the relevant geodesic distances. The uniform distribution of ξ_T was achieved by dividing the range of ξ_T into 10 equally spaced bin ranges containing approximately equal numbers of triangles each, depending on the total number of candidate triangles for each participant. The number of triangles required per bin was set so that the least filled bin was at least 75% full. Those bins that were not full were given an additional weighting in the calculation of spatial power of up to 1.3 (i.e., 1/0.75), depending on how full. For size bins that were over-represented, triangles were excluded from the bin that were furthest from equilateral until the bin was no longer overfull. This procedure of collation of triangles and flattening of the ξ_T distribution resulted a maximum of 9100 triangles for one participant, with participants being excluded from further analysis if they accrued less than 1000 triangles, i.e., 100 triangles in each of 10 bins.

This procedure of forming triplets of approximately equilateral triangles maps the phase matrix, Φ_sc, of dimensions samples × contacts, onto a tensor with dimensions samples × triangles × vertices, Φ_STV. 11 participants had sufficient contact triangles to continue analyses. Higher-order SVD decomposes a tensor, 𝒜, into a smaller core tensor, 𝒮, and a series of mode matrices, U_m, which play the same role as singular vectors in SVD.

where x is matrix-tensor multiplication. Higher-order SVD shares the desirable properties of SVD: it enables construction of a version of the phase data with only the most important components (reduced rank) as well as producing components that are uncorrelated (orthonormal). We decompose the phase/triangle tensor, Φ_STV, into a core tensor with reduced samples and reduced triangles, 𝒮_stV, along with its associated mode matrices:

where the subscripts indicate the dimensionality of the tensor: samples, S; triangles, T; and vertices, V. Indices s and t were the reduced-sample and reduced-triangle dimensions. The second dimension of the mode matrices is left implicit in this indexing. No further reduction or orthogonalization was required for the vertices dimension.

The higher-order SVD was computed by flattening the tensor along the non-samples dimensions, so that the tensor Φ_STV becomes the matrix Φ_SM where the size of M is the product of the other two dimension sizes, T and V. The SVD is performed upon Φ_SM:

and the first core tensor, 𝒮_sTV, is computed as the flattened matrix S_sM from the first s components of V^∗:

Likewise the second core tensor is computed by flattening 𝒮_sTV along the non-triangles dimensions, as S_TM, then:

where the size of M is the product of the other two dimension sizes, s and V. We store the first t components of U as U_t.The second core tensor 𝒮_stV, is computed via the flattened matrix S_tM:

Then for each reduced triangle dimension, τ ∈ {1,…t}, we project the second core tensor back onto the space of triangles as:

where T_τ is the same size as the original set of triangles, but refers to the model components produced from the τ^th reduced triangle dimension. This resulted in t model components, each of which consisted s sets of model phase triplets.

The pair of (1) a model phase triplet from the model components with (2) the three inter-contact distances from the original set of triangles is called a model triangle. Model triangles therefore consist of a directed flow of phase over known distances in metres. In total, there were t*s model triangles. We compute the rate of change in phase per metre on the triangle in the plane implied from the three inter-contact geodesic distances. Hence the spatial frequency, ξ, is given by:

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

The weighting of each model triangle on each model component is given by:

where |X| is the complex magnitude of X. Combining the spatial frequency estimate for each model triangle with the weighting of each model triangle on each model component allows the construction of the estimated spatial frequency spectrum. For K spatial frequency bins, we have a bin width, w:

and the estimated power (A.U.) for the k^th bin, k ∈ [1, K],

We also compute the fit of the model components to the phase data. This first required reconstructing a reduced rank version of the phase/triangle tensor from the first s and t components of U_S and U_T:

In order to illustrate the ability of the higher order SVD to successively identify components at different spatial frequencies, we generated artificial phase data on a one-dimensional cortex. The choice of a one-dimensional cortex was for didactic reasons; it is easier to illustrate and convey the method as a prelude to the introduction of the case for real cortex. The individual samples in the signal comprised sums of sinusoids with wavenumbers k ∈ {0.5,1.0,2.0,3.0,4.0} plus some random phase at each site. The spectral power of the samples was varied so that some wavenumber had twice the power of the wavenumber with next highest power. Analyses were identical to that described (Equations 6-19), except that triangles of contacts were replaced with pairs of contacts for the one-dimensional case. The size of the cortex was 1000 locations nominally over 50cm, with 500 pairs of recording sites with distances drawn from a uniform distribution of 1/f_s. 10,000 time samples were generated. Figure 7 shows the spatial distribution of phase for four typical samples, along with the phase and relative phase distributions over all samples for one pair of sites. The methods were successfully able to extract the low spatial frequency component when that dominated the signal (Figure 7 A-B). The same is true when medium spatial frequency components were the dominant signal Figure 7 C-D).

Figure 7:

For some analyses, we calculated the fit of the model components to the phase data, ρ, which was then computed as:

Where ℜ(X) is the real part of X, |𝒳| is the complex magnitude of 𝒳, and ⟨ 𝒳 ⟩ _Y is the mean over the dimension Y. Where Φ_STV and (unit length) agree, their ratio will be a complex unit vector pointing towards zero radians on the unit circle, giving real part close to one. Where Φ_STV and (unit length) have random relationships in the their argument (angle), their ratio averaged over vertices and triangles will have a real part approximately equal to zero.

For some analyses we computed the fit from a limited spectral range of model triangles. The macroscopic model triangles were defined as those that had a peak spatial power in the range 2.0 to 2.9 cycles/m, the mesoscopic 5.7 to 53 cycles/m (λ=0.35 to 0.5m vs. 0.18 to 0.019m). This resulted in spectra peaked near the middle of these respective ranges (see Figure 5A, for an example). For one participant with a small maximum inter-contact distance, there were no macroscopic model components (as defined here) at some temporal frequencies. Another participant likewise had no mesoscopic components at some temporal frequencies, resulting in 9 participants for this analysis.

We aimed to test the effects of array size by manipulating the size of triangles available to the higher order SVD. To this end, participants’ triangles were partitioned into 5, 10 or 20 equally sized bins along the linear distance range (square root of triangle area). The number of triangles per bin were chosen to be in a similar range across participants, resulting in bins containing 250 to 580 triangles each, depending on the total number of candidate triangles for each participant and the number of bins. This procedure of partitioning the triangles into equally size bins resulted a maximum of 9100 triangles over 20 bins for one participant, with participants being excluded from further analysis if they accrued less than 1000 triangles, i.e., 200 triangles in each of 5 bins. The triangle collation procedure resulted in 10 participants available for this analysis. For each participant, and for the triangles in each bin, a higher order SVD was used to analyze the tensor of phase, as described. The mean linear distance of the triangles for each run was compared to the spatial frequency spectra so produced.

Statistical testing of linear regressions reported in the results used the Wald test, for which the null hypothesis is that regression slope is zero. Fits, ρ, were correlated across time samples for the cross-frequency correlations. The significance of these correlations was tested by permutation methods, in which each permutation of the data comprised a correlation between the first time series, as is, and the second time series, offset by a random number of samples with circular boundary conditions. The reason for this choice of shuffling was that the time series were locally (temporally) correlated due to the way phase is calculated, and this local correlation is stronger for low temporal frequencies. The shuffling method preserves the degree of local temporal correlation, while having a null hypothesis of zero correlation.

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.

Supplemental information

Supplementary section A Methods

Unless otherwise detailed here, methods were the same as for the main analysis. The purpose of this supplemental analysis was to show that previously published methods of phase map decomposition (8,11,28) reproduced the substantive results reported in the main analysis. Phase time-series were estimated using short time window Morlet wavelets (2 cycles) over the frequency range 2 to 42Hz. The analysis used all 35 participants, since the alternative method reported here did not require as many grey-matter contacts as the method reported in the main analysis. As with the main analysis, the irregularity of the stereotactic EEG contact array makes standard methods of computing spatial frequency untenable. We therefore computed the singular vectors of each participant’s complex-valued, unit length phase vectors, at each frequency. Example code is given here: https://github.com/ScaleSymmetry/Traveling-wave-analysis/blob/main/example_morlet_wavemap.ipynb This meant we could characterize the phase in terms of spatial maps, with each map explaining a known amount of variance within the samples × contact phase array. The maps of phase derived by these methods show a range of wave-types, typically some combination of standing wave and traveling wave (11). In the present reported results, we discarded all but the first 12 singular vectors. Testing of the number of singular vectors to be included over a range (6-24) indicated the results were not dependent on this choice.

Standing waves with a small traveling component can result in erroneous estimations of the velocity and wavelength of the component wave, unless care is taken to separate out the pure traveling wave components comprising the wave. This is because, near the spatial point of transition between peak and trough, the apparent local rate of change of phase is inflated by the cancellation of the two component TWs (see Supplementary Figure 1). We therefore estimated the underlying TW components for each singular vector, using the following numerical procedure: Each phase map was extended from representing the TW as a single sample of values to encompass the full cycle of a wave (11). The cycle × contact matrix of complex-valued phases was then decomposed into pure TW components using SVD. Essentially this procedure forces the decomposition to construct the underlying wave components that sum to the phase map, i.e., the pure TW versions of the phase map. Example code is given here: https://github.com/ScaleSymmetry/Traveling-wave-analysis/blob/main/example_velocity_decomposition.ipynb

We then applied local difference methods (i.e., change in phase per change in distance) in order to estimate spatial frequency. Estimates of local change in phase were made over areas defined by a triangle of three nearest neighbour contacts. The reason for this was that estimating between pairs of contacts can overestimate the phase velocity and underestimate the spatial frequency, if the actual direction of the wave is orthogonal to the axis of the contact pair. Prior to estimation, therefore, the geodesic distances (40) between contacts were calculated, and these distances used to find triplets of nearest neighbour contacts. From these triplets we collated the set of local triangles, and triangles were excluded if they were too obtuse (defined here as having a smallest inter-contact distance less than 50% of largest inter-contact distance). For each triangle in the residual set, the local change in phase per metre was computed, using the phase values from the pure TW maps computed in the previous step and the geodesic distances. The weighting of each triangle in the original phase data was also computed (see main methods), allowing the construction of a histogram of spatial frequency versus estimated power.

Supplementary Figure 1.

Results

We find that stereotactic EEG spectral power decreases with increasing spatial frequency. This is true over a broad range of temporal frequencies. The aggregate spatial frequency spectrum for all 35 participants is shown in Supplementary Figure 2. The spatial frequency spectrum shows a similar shape to that found for the main results (see Figure 4C). The slopes of the spatial frequency spectra were estimated for the temporal frequency 9.2Hz, for each participant. On a log-log scale, the fits were linear (Supplementary Figure 3), but the slopes depended strongly on the sampling details of each participant’s stereotactic EEG contact placements. This dependence means the exponent of the 1/f^α distribution were not solely a function of the underlying cortical phase dynamics, so only the general form of the distribution is inferred.

Supplementary Figure 2.

Supplementary Figure 3.

Additional Supplementary Figures

Supplementary Figure 4.

Supplementary Figure 5.

Supplementary Figure 6.