The dominance of large-scale phase dynamics in human cortex, from delta to gamma
Figures
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 10 cm or less.(C) Stereotactic EEG measures electrical activity within the gray matter of the cortex, and therefore has a lower degree of blurring. The clinical placement of the linear arrays of contacts results in a highly irregular measurement array, though often with wide coverage of the cortex, though in general less than 30 cm. (D–F) Hypotheses for the origin of macroscopic TWs measured extracranially (adapted from Hindriks et al., 2014). (D) Extracranial TWs reflect real, coordinated, spatio-temporal activity across the entire cortex. The TWs are measurable at the scalp because the dominant SF 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 SF spectrum, when measured in the gray matter, has a peak at a low SF. (E) A local wave (shown here in a sulcus) slowly traverses a small region of gray 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 SF spectrum (magenta), when measured in the gray matter, has a peak at a high SF. (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 SF spectrum (magenta), when measured in the gray matter, has a peak at high spatial frequencies corresponding to the size of each oscillatory source, with a weaker elevation in low SF power due to the interaction between sources.
Schematic of data processing pipeline.
(A) Cortical reconstruction of one participant, showing placement of sEEG contacts. Gray 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 filtered at 9.2 Hz (real part only) are shown for all gray matter contacts in one participant. One contact’s phase values are shown in darker gray. 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 center 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 gray matter contacts within inflated left hemisphere of cortex. Phase angles of a LSV indicated on the rainbow scale. Phase angles are shown on the color wheel. (E) Flattened representation of the approximate equilateral triangles of gray 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 a gray arrow, and phases have been interpolated between vertices of the triangle for visualization purposes. (G) The rate of phase change of a triangle of contacts is converted to estimated SF 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, that is 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 SF 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 15 Hz 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). The 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 colored horizontal line, with phase values shown on the color wheel. The first plot shows the dominant ‘spatial’ (horizontal) frequency, captured as a single cycle wave. Each successive component has an increment in the 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 center TF of 9.2 Hz. One sensor’s data is highlighted in darker gray. (F) Wave maps (Alexander et al., 2006) extracted from the participant’s phase data at 9.2 Hz. The MEG sensors are mapped to a sphere, and phase values are interpolated (values indicated by the color 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 (Alexander et al., 2006; Alexander et al., 2013). These global patterns of phase are used to extract the dominant SF of TWs in the data (Alexander et al., 2019; Alexander et al., 2016). 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. (Colored lines below each map correspond to the spectra lines in Figure 4B).
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 50 cycles/m only) and associated spectra are shown with colored 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, aggregated by triangle bin size. (F) Smoothed spectra for (E).
Surrogate analysis of multi-scale contact arrays using random walk data.
(A) A random walk (Nt=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 FFT. Since the data are one-dimensional they can be regarded as temporal (frequencies in Hz) or as sampled in one spatial dimension (frequencies in cycles/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 spaced 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 are strongly a function of the characteristic intercontact distances in the multi-scale arrays. The LSVs from the original SVD (Figure 3A) are shown as gray 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 neighborhood (purple line) also reveals that low SF components are represented in the LSVs. Plots are the same as (F), row 1, panels 2 and 3 (indicated by purple squares).
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 with discrete numbers of cycles, the Fourier spectrum can simply be estimated by reading off the singular values (shown here as % variance explained). The dominant SF 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 colors correspond to the line labels in Figure 3F, and 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 participants.
(A). Cortical surface with phase values of the first five LSVs from a single participant at center TF 3.5 Hz. Gray matter contact locations are shown with spheres, and the phase value of the LSV at each site is shown on color 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–97.0 Hz (color 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 TFs (colored lines), for the participant shown in (D) (only regression lines significant at p<0.05 shown). Spectrum values shown as dots, colours indicating TF. The mean regression line (mean of slopes and offsets of significant regressions) is shown with a black line and extended to 0.32 m 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 (colored lines), for the participant shown in (G). Conventions the same as (F).
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 (colored 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 and D) Same data as (A and 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 TFs 1.0–97.0 Hz (color scale), for all participants.
Participant 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 TFs 1.0–97.0 Hz (colored lines), for all participants.
Participant 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 Figure 5—figure supplement 2 (only individual regression lines significant at p<0.05 shown). Spectrum values shown as dots, colours indicating TF. 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.32 m wavelength for visualization purposes.
Surrogate tests on individual participants.
(A) Effects of injection of veridical SF surrogate signal on the empirically derived spectrum. The left panel shows the spectrum of a pure surrogate signal (S) with an anterior-posterior gradient and a SF at 4 cycles/m (top row within plot). The empirical spectrum for the phase of this participant at a center frequency of 9.2 Hz is also shown (bottom row within plot). The effects of different weightings of empirical and surrogate signals vary vertically over the plot. The middle panel follows the same conventions, except a surrogate signal at 9 cycles/m is shown. The right panel shows the effect of a surrogate signal at 14 cycles/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.2 Hz 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.
More surrogate testing examples.
(A) Participant’s inflated cortex with surrogate phase values, here in the form of a left-right gradient at 4 cycles/m. Gray matter contacts are shown as spheres, with phase shown on the color scale. Colored lines below each cortex refer to the color spectra in (D,E and F,B). Same conventions as (A), except 9 cycles/m. (C) Same conventions as (A), except 14 cycles/m. (D) Spectra for pure surrogate signals injected as a left-right gradient to the participant in (A), over a range of SFs 2–16 cycles/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. The left panel shows the spectrum of a pure surrogate signal (S) with the left-right gradient and a SF at 4 cycles/m (top row within plot, same surrogate signal as shown in A). The empirical spectrum for phase of this participant at a center frequency of 9.2 Hz is also shown (bottom row within plot). The effects of different weightings of empirical and surrogate signal vary vertically over the plot. The middle panel follows the same conventions, except a surrogate signal at 9 cycles/m is shown (same surrogate signal as shown in B). The right panel shows the effect of a surrogate signal at 14 cycles/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.21 m (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 SF spectrum of a single participant for TFs 1.0 (blue) to 32.0 Hz (red), estimated from the sensory array (and corresponding time-series) shown in (A). (E) Stacked SF spectrum of a single participant for TFs 1.0 (blue) to 32.0 Hz (red), estimated from the sensory array (and corresponding time-series) shown in (C). (F) Wavelength by SF power regression lines for two MEG participants at 9.2 Hz, 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 SF spectrum from the surrogate MEG data with reversed SF spectrum, conventions same as (D) (all sensors). (H) Stacked SF spectrum from the surrogate MEG data with reversed SF spectrum, conventions same as (E) (two-thirds of sensors removed). (I) Wavelength by SF power regression lines for two MEG participants with reversed SF spectrum at 9.2 Hz, conventions same as (F).
Surrogate testing for the effect of average reference and removal of DC component on the phase spatial spectrum.
(A) SF spectrum of phase for a single MEG participant, at 1 Hz with increasing average reference (color scale). DC component has not been removed. Black dot shows the SF corresponding to the maximum inter-sensor distance. (B) Same conventions as A, except 9.2 Hz. (C) Same conventions as A, except 32.0 Hz. (D) SF spectrum of phase at 1 Hz with increasing average reference. DC component has been removed. (E) Same as (D), except 9.2 Hz. (F) Same as (D), except 32.0 Hz.
Illustration of the removal of negative SF components.
(A) Three LSVs from an MEG subject at 6.0 Hz. 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 (color 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.
Cross-participant analyses.
(A) Regression lines, for all participants, for wavelength versus SF power at center TF 9.2 Hz. 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 colored dots (data from Figure 5—figure supplement 2). Only individual regressions at 9.2 Hz 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 colored scale. Individual lines are the mean spectra over 23 participants at that TF. The maximum triangle size (root area, 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 (gray dots) and SF at peak power versus TF (black crosses; data from B). Significant regression line for peak power is shown in gray. (E) Single participant regression lines for peak power versus TF. Color conventions the same as (A), and the mean slope and offset of the significant regressions are 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. Color conventions the same as (A).
Cross-participant analyses by TF.
Regression lines for (logarithm of) wavelength versus SF power at each center TF, for all participants. In each plot, participants are ordered from most number of triangles (n=11,051; light blue) to least number of triangles (n=56; lilac). Spectra values are shown as colored dots (data from Figure 5—figure supplement 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 TF are shown as a black line.
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 color 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 colors, 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 color 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.
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 color 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 colors. In order to visualize the SFs at higher TFs, the SF power has been normalized between 0 (cold) and 1 (hot) at each TF, 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 color 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.
Videos
Dorsal view of one participant’s dominant spatial patterns of phase at 3.5 Hz.
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 6 s of cortical activity. Color values are cosine of the phase (yellow=-1, red = 1). Values between contacts are interpolated for visualization purposes. The anterior region of the cortex is to the left.
Anterior view of one participant’s dominant spatial patterns of phase at 8.0 Hz.
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 6 s of cortical activity. Color values are cosine of the phase (yellow=-1, red = 1). Values between contacts are interpolated for visualization purposes. The left hemisphere of the cortex is to the left.
Left view of one participant’s dominant spatial patterns of phase at 8.0 Hz.
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 6 s of cortical activity. Color values are cosine of the phase (yellow = -1, red = 1). Values between contacts are interpolated for visualization purposes. The anterior region of the cortex is to the left.