Traveling waves, measurements and origin.

A-C. Scales of measurement and degree of blurring in EEG, ECoG and stereotactic EEG. The maximum linear distances of unfolded cortex typically covered by each technique are shown in each panel. A. EEG measures the electrical activity of the cortex at the surface of the scalp (contacts shown in green). Due to volume conduction in the intervening tissues (blue), the cortical signal is highly blurred and preferentially samples gyri. The EEG array covers the entirety of the scalp in a quasi-regular grid, enabling a large extent of cortical coverage (green line, not to scale). B. Electrocorticogram (ECoG) measures the electrical activity of the cortex at the cortical surface. The volume conduction effects are less here, but the measurement array preferentially samples activity at the gyri. The distance between contacts is regular, forming a two-dimensional measurement array, but the size of the array is limited to 10cm or less. C. Stereotactic EEG measures electrical activity within the grey matter of the cortex, and therefore has a lower degree of blurring. The clinical placement of the linear arrays of contacts result in a highly irregular measurement array, though often with wide coverage of the cortex, though in general less than 30cm. D-F. Hypotheses for the origin of macroscopic TWs measured extra-cranially (adapted from 14). D. Extra-cranial TWs reflect real, coordinated, spatio-temporal activity across the entire cortex. The TWs are measurable at the scalp because the dominant spatial frequency of the activity has a long wavelength and is therefore detectable at the scalp. A sinusoidal wave corresponding to this hypothesis is shown in blue, with successive time samples in lighter shades. The fast Fourier transform of the model wave is shown in magenta. The expected spatial frequency spectrum, when measured in the grey- matter, has a peak at a low spatial frequency. E. A local wave (shown here in a sulcus) slowly traverses a small region of grey matter via intracortical connections. The wave appears at the scalp as a fast-moving global wave due to volume conduction effects. The wave’s (blue) expected spatial frequency spectrum (magenta), when measured in the grey-matter, has a peak at a high spatial frequency. F. Localized oscillatory sources that are out of phase appear as large-scale TWs at the scalp due to volume conduction effects. The oscillators’ (blue) expected spatial frequency spectrum (magenta), when measured in the grey-matter, has a peak at high spatial frequencies corresponding to the size of each oscillatory source, with a weaker elevation in low spatial frequency power due to the interaction between sources.

Schematic of data processing pipeline.

A. Cortical reconstruction of one participant, showing placement of sEEG contacts. Grey matter contacts are shown as black dots. White matter contacts are shown with colored circles (blue, or purple for contacts with 50% lowest amplitude signal). B. The phase values of a segment of data at filtered at 9.2Hz (real part only) is shown for all grey matter contacts in one participant. One contact’s phase values are shown in darker grey. A single vector of phase is indicated with the yellow line. C. SVD decomposes any matrix, A, into three separate matrices consisting of a rotation (U, the LSVs), a scaling (Σ) and another rotation (V, the right singular vectors). Here, the matrix of complex-valued phases for a participant at some centre frequency, 𝚽𝑠𝑡, is decomposed into LSVs that represent typical patterns of phase (red vertical bars). The values along the diagonal of Σ give the overall weighting of each LSV in the original matrix. The right singular vectors give the changing time- course for each of the LSVs in the original matrix (blue vertical lines). The bright and pale reds, greens and blues indicate how the least important components (pale) may be dropped, combining a small number of most important components to form a good approximation of the original matrix. D. Placement of grey matter contacts within inflated left hemisphere of cortex. Phase angles of a LSV indicated on the rainbow scale. Phase angles are shown on the colour wheel. E. Flattened representation of the approximate equilateral triangles of grey-matter contacts, into which the LSV values are resampled as triplets of values. F. Triplets of values from each LSV were used to estimate the rate of change of phase over regions of the cortex. The direction of phase change is indicated with grey arrow, and phases have been interpolated between vertices of the triangle for visualization purposes. G. The rate of phase of change of a triangle of contacts is converted to estimated spatial frequency over that region by incorporating the geodesic distances into the numerator. The inflated cortex is shown here, but the geodesic distances are calculated from the folded cortex i.e. shown in A.

Relation between SVD and Fourier analysis.

A. A random walk of t time-steps is first filtered to remove very low frequency components, then successive segments (w=1024 samples) from the time-series are windowed with a Gaussian to produce a Hankel matrix 𝐀𝑤𝑇, where T = t - w/2. Here, the one unwindowed sequence of the filtered random walk is shown. B. The LSVs of the matrix 𝐀𝑤𝑇 come in pairs, with each pair approximating the real and imaginary parts of a one cycle Morlet wavelet. Here, the first four such pairs are shown, showing increasing frequency of the empirically derived Fourier components. The frequency of the first pair is closely tied to the size of the windowed sampling array. C. Image of waves. The waves have a peak spatial frequency at 15 cycles per 1000 pixels vertically, and this dimension is treated as if it were time by scanning down the image (nominally) at 1ms per pixel. The waves in the image are aligned diagonally, moving towards the bottom right. The waves complete approximately one cycle every 512 pixels in the horizontal direction. This is the width chosen for the size of the measurement array, scanning the image vertically three times (white boundaries). The complex-valued phase at 15Hz along the ‘temporal’ (vertical) axis is first calculated, and the phase dynamics along ‘spatial’ (horizontal) dimension are extracted using SVD. D. First four components of the SVD of the vector of phases from image in C. Real part of the component is shown in blue, imaginary part in orange. The phase angle corresponding to these real and imaginary values is shown by the coloured horizontal line, with phase values shown on the colour wheel. The first plot shows the dominant ‘spatial’ (horizontal) frequency, captured as a single cycle wave. Each successive component has an increment in number of cycles represented, from the previous. E. Time-series of (the real part of) phase extracted from MEG data recorded with 151 magnetometers, here filtered at a centre TF of 9.2Hz. One sensor’s data is highlighted in darker grey. F. Wave maps (4) extracted from the participant’s phase data at 9.2Hz. The MEG sensors are mapped to a sphere, and phase values interpolated (values indicated by the colour wheel). The head is viewed from the top, with the nose pointing left. The first three maps show global patterns of phase, with directions of flow characteristic of MEG data, namely left-right, inferior-superior and anterior-posterior (4,29). These global patterns of phase are used to extract the dominant SF of TWs in the data (8,11). The components in this type of analysis are inherently spatio-temporal, since a spatial vector of smoothly changing phase implies a traveling wave over time. Successive groups of LSVs indicate components of the signal with increasing SF. (coloured lines below each map correspond to the spectra lines in Figure 4B)

SVD and estimating spectral power.

A. Singular values obtained from SVD of the phase vectors from the image shown in Figure 3C. Because the first four components from the SVD are approximate sinusoids, the Fourier spectrum can simply be estimated by reading off the singular values (shown here as % variance explained). The dominant spatial frequency is at wavenumber equals one cycle per 512 pixels (the size of the array), and the next highest SF components contribute successively less to the signal. B. Estimated phase spectra for the wave maps of MEG data shown in Figure 3F, computed by the present method. The line colours correspond to the line labels in Figure 3F, the black dot indicates the SF equivalent to the maximum distance between sensors. The spectra peak at a range of different SFs, with the global maps in Figure 3F showing the lowest SFs and highest power, and maps 6-14 showing the highest SF and least power.

LSVs of sEEG participants, equilateral triangles of LSV phase and estimated SF spectra for two partcipants.

A. Cortical surface with phase values of the first five LSVs from a single participant at centre TF 3.5Hz. Grey-matter contact locations are shown with spheres, and the phase value of the LSV at each site is shown on colour scale. Phase values on the cortical surface are interpolated between contacts for visualization purposes. B. Examples of approximately equilateral triangles of contacts, used to calculate the SF spectrum. The shortest geodesic distance between contacts is shown with the black line. Phase on the section of cortical surface is interpolated for visualization purposes. The scale of each equilateral triangle (square-root of triangle area) is indicated with the black reference bar. C. Same examples shown in B, except projected onto the inflated cortex to show the unfolded triangle. D. SF spectra for one participant, over TFs 1.0 to 97.0Hz (colour scale). The largest triangle size is indicated with the black dot. E. Same as D, except shown as a stacked spectrum. F. Regression lines for (logarithm of) wavelength versus SF power at all temporal frequencies (coloured lines), for the participant shown in D (only regression lines significant at p<0.05 shown). Spectrum values shown with coloured dots. The mean regression line (mean of slopes and offsets of significant regressions) is shown with a black line and extended to 0.32m wavelength for visualization purposes. G. SF spectra for a second participant. H. Same as G, except shown as a stacked spectrum. I. Regression lines for wavelength versus SF power (coloured lines), for the participant shown in G. Conventions the same as F.

Surrogate tests on individual participants.

A. Effects of injection of veridical SF surrogate signal on the empirically derived spectrum. Left panel shows the spectrum of a pure surrogate signal (S) with an anterior-posterior gradient and a SF at 4c/m (top row). The empirical spectrum for phase of this participant at a centre frequency of 9.2Hz is also shown (bottom row). The effects of different weightings of empirical and surrogate signal vary vertically over the plot. Middle panel follows the same conventions, except a surrogate signal at 9c/m is shown. Right panel shows the effect of a surrogate signal at 14c/m. B. SF spectra of phase for a participant with MEG measurements. Each spectrum shows the characteristic monotonic decrease in power with increasing SF. This data is from the same participant as Figure 4B, where the unsummed spectra for each LSV at 9.2Hz are shown. C. SF spectra for surrogate MEG measurements. Here, the singular values of the first 28 components are reversed in order before reconstructing the phase time-series, creating phase signal in which high SFs dominate. The plot shows that the present methods successfully characterize the surrogate signal as SF spectra with slope reversed from the empirical signal.

Cross-participant analyses.

A. Regression lines, for all participants, for wavelength versus SF power at centre TF 9.2Hz. Participants are ordered from most number of triangles (n=11051; light blue) to least number of triangles (n=56; lilac). Spectra values are shown as coloured dots (data from Supplementary Figure 2). Only individual regressions at 9.2Hz that were significant at p<0.05 are shown as regression lines. The mean slope and offset (of only significant regressions) is shown as a black line. B. Participant aggregate spectra for all TFs. TF is shown on the coloured scale. Individual lines are the mean spectra over 23 participants at that TF. The maximum triangle size (expressed as equivalent SF, i.e. 1/m) is shown for all 23 participants, as black dots. C. The same spectra as B, except shown as a stacked histogram, giving the overall shape of the spectra. D. Regression of peak power versus TF (grey dots) and SF at peak power versus TF (black crosses; data from B). Significant regression line for peak power is shown in grey. E. Single participant regression lines for peak power versus TF. Colour conventions the same as A, and the mean slope and offset of the significant regressions is shown as a black line. Peak power is shown on a logarithmic scale for visualization purposes only. F. Single participant regression lines for SF at peak power versus TF. Colour conventions the same as A.

Effects of including normalized TF power on SF spectra.

A. Participant’s stacked SF spectra with normalized TF power included in the estimate. TFs are shown on the colour scale. Compare with the same participant, phase only spectra, Figure 5H. B. Aggregate stacked SF spectra with normalized TF power included, over 23 participants. Compare with the aggregate stacked spectra shown in Figure 7C. C. Participant’s SF spectra as SF by TF plot. The left plot shows the SF spectra with normalized TF power included (data taken from A). Power values are shown as colours, normalized between 0 (cold) and 1 (hot). The middle plot shows the phase only SF spectra (data taken from Figure 5H), normalized between 0 (cold) and 1 (hot). The right plot shows the difference between the two SF spectra, with differences given on the colour scale. D. Aggregate spectra for 23 participants, conventions the same as C. Data are taken from B and Figure 7C for the left and middle plots.

Examples of the equilateral triangles between contacts used in the multi-scale differencing.

A. Cortical surface for one participant, with contacts (black dots) and geodesic distances (coloured lines) corresponding to edges of triangles in the largest triangle bin used in the multi-scale differencing. B. Contacts and geodesic distances corresponding to the smallest triangle bin for this participant. C & D. Same data as A & B but shown on the inflated cortex. While the geodesic distances are calculated on the folded cortex, the inflated cortex gives a clearer view of the equilateral triangles connected by shortest distances across the cortex. E. Distribution of triangle sizes (as equivalent SFs) for this participant. Left green wedge and right red wedge indicate the bounds of each triangle bin. F. Normalization factor applied to each triangle of contacts in the calculation of the SF spectrum. This effectively flattens the distribution to give equal weighting to triangles across the size range.

Stacked SF spectra, across temporal frequencies 1.0 to 97.0Hz (colour scale), for all participants.

Participant’s plots are ordered by number of equilateral triangles, highest (n=11051) to lowest (n=66). The largest triangle size for each participant is indicated with the black dot.

Regression lines for wavelength versus SF power at temporal frequencies 1.0 to 97.0Hz (coloured lines), for all participants.

Participant’s plots are ordered by number of triangles, from highest (n=11051) to lowest (n=66). Regressions are from the values of the spectra shown in Supplementary Figure 2 (only individual regression lines significant at p<0.05 shown). Spectrum values are shown with coloured dots. The mean regression line for each participant (mean of slopes and offsets of significant regressions) is shown with a black line and extended to 0.32m wavelength for visualization purposes.

More surrogate testing examples.

A. Participant’s inflated cortex with surrogate phase values, here in the form of a left-right gradient at 4c/m. Grey-matter contacts are shown as spheres, with phase shown on the colour scale. Coloured lines below each cortex refer to the colour spectra in D, E & F.B. Same conventions as A, except 9c/m. C. Same conventions as A, except 14c/m. D. Spectra for pure surrogate signals injected as a left-right gradient to the participant in A, over a range of SFs 2 to 16 c/m. Black dot shows maximum array distance for this participant. E. Same as D, except for evenly weighted surrogate and empirical signals. F. Same as E, except for 1:2 surrogate to empirical signal. G. Effects of injection of veridical SF surrogate signal on the spectrum derived from the empirical signal. Left panel shows the spectrum of a pure surrogate signal (S) with the left-right gradient and a SF at 4c/m (top row, same surrogate signal as shown in A). The empirical spectrum for phase of this participant at a centre frequency of 9.2Hz is also shown (bottom row). The effects of different weightings of empirical and surrogate signal vary vertically over the plot. Middle panel follows the same conventions, except a surrogate signal at 9c/m is shown (same surrogate signal as shown in B). Right panel shows the effect of a surrogate signal at 14c/m (shown in C). H. Participant averaged spectra (n=23), conventions otherwise the same as G.

Effects of randomly removing sensors from a semi-regular array.

A. MEG sensor array, with 151 sensors, showing triangles of sensors in the size bin with mean size 0.21m (smallest and largest triangle in this bin without transparency). B. MEG sensor array, with one third of sensors removed, conventions same as A. C. MEG sensor array, with two thirds of sensors removed, conventions same as A. D. Stacked spatial frequency spectrum of a single participant for temporal frequencies 1.0 (blue) to 32.0Hz (red), estimated from the sensory array (and corresponding time-series) shown in A. E. Stacked spatial frequency spectrum of a single participant for temporal frequencies 1.0 (blue) to 32.0Hz (red), estimated from the sensory array (and corresponding time-series) shown in C. F. Wavelength by spatial frequency power regression lines for two MEG participants at 9.2Hz, shown with closed and open circles. Light blue shows the regressions for full sensor arrays, violet shows the regression lines for the sensor arrays with two thirds of sensors removed. Contacts were removed pseudo-randomly in five steps between these ranges (inclusive), in three runs. Only regressions that were significant (p<0.001) are shown as lines. Mean slope over all regressions shown with black line (p<0.001, permutation testing). G. Stacked spatial frequency spectrum from the surrogate MEG data with reversed spatial frequency spectrum, conventions same as D (all sensors). H. Stacked spatial frequency spectrum from the surrogate MEG data with reversed spatial frequency spectrum, conventions same as E (two thirds of sensors removed). I. Wavelength by spatial frequency power regression lines for two MEG participants with reversed spatial frequency spectrum at 9.2Hz, conventions same as F.

Surrogate testing for the effect of average reference and removal of DC component on the phase spatial spectrum.

A. Spatial frequency spectrum of phase for a single MEG participant, at 1Hz with increasing average reference (colour scale). DC component has not been removed. Black dot shows the spatial frequency corresponding to the maximum inter-sensor distance. B. Same conventions as A, except 9.2Hz. C. Same conventions as A, except 32.0Hz. D. Spatial frequency spectrum of phase at 1Hz with increasing average reference. DC component has been removed. E. Same as D, except 9.2Hz. F. Same as D, except 32.0Hz.

Cross-participant analyses by temporal frequency.

Regression lines for (logarithm of) wavelength versus SF power at each centre TF, for all participants. In each plot, participants are ordered from most number of triangles (n=11051; light blue) to least number of triangles (n=56; lilac). Spectra values are shown as coloured dots (data from Supplementary Figure 2). Only individual participant (at a given frequency) regressions that were significant at p<0.05 are shown as regression lines in each plot. The mean slope and offset for all significant regressions at the relevant temporal frequency are shown as a black line.

Effects of including TF power on SF spectra.

A. Participant’s stacked SF spectra with TF power included in the estimate. TFs are shown on the colour scale. Compare with the same participant, phase only spectra, Figure 5H. Higher TFs have smaller contributions to SF power due to the smaller amplitudes of higher TF signal. B. Aggregate stacked SF spectra with TF power included, over 23 participants. Compare with the aggregate stacked spectra shown in Figure 7C. C. Participant’s SF spectra as SF by TF plot. Power values are shown as colours. In order to visualize the SFs at higher TFs, the SF power has been normalized between 0 (cold) and 1 (hot) at each temporal frequency, prior to plotting. The left plot shows the SF spectra with TF power included (data taken from A). The middle plot shows the phase only SF spectra (data taken from Figure 5H). The right plot shows the difference between the two SF spectra, with differences given on the colour scale. D. Aggregate spectra for 23 participants, conventions the same as C. Data are taken from B and Figure 7C for the left and middle plots.

Illustration of the removal of negative SF components.

A. Three LSVs from an MEG subject at 6.0Hz. These LSVs showed stronger than usual standing wave patterns in components 3 and 4. The standing wave patterns can arise from SVD, even when the underlying samples are all pure TWs, if the balance of the two directions is almost equal. As discussed in Alexander et al. (2016), standing wave patterns give rise to misleading estimates of phase velocity (here, to SF), since the interference pattern of the two wave directions creates both slowly changing gradients (at the peaks and troughs of the standing wave) and fast changing gradients (at the inflection points of the standing wave). Panels show left, top and back views of MEG array, with phase values (colour scale) interpolated for visualization purposes. B. Same data as A, except the negative SF components have been removed, leaving only the positive SFs. The gradients of phase are smoother, reflecting the SF of the underlying Fourier-like component of the phase. Conventions the same as A.

Effects of bin size and calculation of spectral peak.

A. Spectra for an individual participant, aggregated by triangle bin size. Triangle bin sizes are indicated with black dots (as 1/m, shown out to 50c/m only) and associated spectra are shown with coloured lines. B. Same data as A, except spectra have been smoothed with a rolling average over 10 steps. Peaks of spectra are read off from the maxima of these curves. C. Spectra for another individual participant, conventions same as A. D. Smoothed spectra for participant shown in C. E. Mean spectra across 23 participants, aggregate by triangle bin size. F. Smoothed spectra for E.

Surrogate analysis of multi-scale contact arrays using random walk data.

A. A random walk (t = 1,000,000) is constructed by taking cosine of a randomly generated number between -π and π and adding it to the cumulative total. The series is high- pass filtered to remove the DC offset, trend and large scale fluctuations. B. The frequency spectrum of the filtered random walk, generated using FTT. Since the data are one- dimensional they can be regarded as temporal (frequencies in Hz) or as sampled in one spatial dimension (frequencies in c/m). C. A Hankel matrix was constructed from successive segments of w = 1024 samples, windowed by a Gaussian function (shown). The Hankel matrix 𝐀𝑤𝑇, where T = t - w/2, is submitted to the SVD. D. Multi-scale contact grids (red dots) were constructed by applying modulo arithmetic to an initial sparse but regularly space grid (first panel). The next three grids (panels 2-4) are comprised of separated individual contacts, but the spatial scale of the Figure is not sufficient to reveal them as separated dots. E. Histograms of intercontact distances for the five contact arrays shown in D. F. LSVs from five SVD analyses using each of the measurement arrays to sample the same random walk. Successive pairs of components are shown down each column, showing an increase in SF tuning down the column. Each row shows a pair of components (e.g. 1st and 2nd components, with lowest SF, in the first row), as extracted from the five SVD analyses. The SF properties of the component is strongly a function of the characteristic intercontact distances in the multi-scale arrays. The LSVs from the original SVD (Figure 3A) are shown as grey and black lines for comparison purposes. G. High SF representations are visually apparent in each of the LSVs produced from multi-scale SVDs. Interpolation of a sinusoid from the leading contact in each local neighbourhood (purple line) also reveals that low SF components are represented in the LSVs. Plots are the same as F, row 1, panels 2 & 3 (indicated by purple squares).