EEG state-trajectory instability and speed reveal global rules of intrinsic spatiotemporal neural dynamics

Spatiotemporal dynamics of EEG/MEG (electro-/magneto-encephalogram) have typically been investigated by applying time-frequency decomposition and examining amplitude-amplitude, phase-phase, or phase-amplitude associations between combinations of frequency bands and scalp sites, primarily to identify neural correlates of behaviors and traits. Instead, we directly extracted global EEG spatiotemporal dynamics as trajectories of k-dimensional state vectors (k = the number of estimated current sources) to investigate potential global rules governing neural dynamics. We chose timescale-dependent measures of trajectory instability (approximately the 2nd temporal derivative) and speed (approximately the 1st temporal derivative) as state variables, that succinctly characterized trajectory forms. We compared trajectories across posterior, central, anterior, and lateral scalp regions as the current sources under those regions may serve distinct functions. We recorded EEG while participants rested with their eyes closed (likely engaged in spontaneous thoughts) to investigate intrinsic neural dynamics. Some potential global rules emerged. Time-averaged trajectory instability from all five regions tightly converged (with their variability minimized) at the level of generating nearly unconstrained but slightly conservative turns (∼100° on average) on the timescale of ∼25 ms, suggesting that spectral-amplitude profiles are globally adjusted to maintain this convergence. Further, within-frequency and cross-frequency phase relations appear to be independently coordinated to reduce average trajectory speed and increase the variability in trajectory speed and instability in a relatively timescale invariant manner, and to make trajectories less oscillatory. Future research may investigate the functional relevance of these intrinsic global-dynamics rules by examining how they adjust to various sensory environments and task demands or remain invariant. The current results also provide macroscopic constraints for quantitative modeling of neural dynamics as the timescale dependencies of trajectory instability and speed are relatable to oscillatory dynamics.


Introduction
EEG/MEG offer non-invasive methods to study large-scale neural dynamics with millisecond temporal precision (see Varela et al., 2001 for a review). For the purpose of investigating neural dynamics, EEG/MEG from each scalp site is typically decomposed into a time series of frequency components (e.g., by using a Morlet wavelet convolution). Temporal associations among spectral amplitudes and/or phases are then examined within each site or between pairs of sites within the same frequency bands (cross-site only) or across different frequency bands (within-and cross-site). Temporal associations have also been examined between the phase of a lower-frequency band (typically q-a) and the amplitude of a higher-frequency band (typically g) either within or across sites. These analyses using time-frequency decomposition have elucidated how amplitude-amplitude, phase-phase, and/or phase-amplitude couplings within and across sites in different brain regions may contribute to various perceptual, attentional, and cognitive functions (e.g., Fries, 2005;Palva & Palva, 2007;Duzel, Penny, & Burgess, 2010;Engel & Fries, 2010;Dahaene & Changeux, 2011;Donner & Siegel, 2011;Palva & Palva, 2012;Klimesch, 2012;Jensen, Gips, Bergmann, & Bonnefond, 2014;Michalareas, Vezoli, van Pelt, Schoffelen, Kennedy, & Fries, 2016).
The current study was motivated by a complementary goal. Instead of trying to identify various types of oscillatory coupling as neural correlates of mental activities, we investigated general rules that may govern large-scale spatiotemporal dynamics. We focused on "intrinsic" dynamics reflected in EEG signals recorded while people rested with their eyes closed (while asked to allow spontaneous thoughts). Many EEG/MEG studies (also fMRI studies based on inter-region temporal correlations of BOLD signals on a much slower timescale) have examined what is called resting-state activity for diagnostic classification (e.g., Adebimpe, Arabi, Bourel-Ponchel, Mahmoudzadeh, & Wallois, 2016;Antonakakis et al., 2016;Gonzalez et al., 2018) as well as for identifying intrinsic spatial networks of phase-phase, amplitude-amplitude, and/or phase-amplitude couplings within or between specific frequency bands (e.g., Osipoval, Hermes, & Jensen, 2008;Hipp, Hawellek, Corbetta, Siegel, & Engel, 2012;Müller, Perdikis, Oertzen, Sleimen-Malkoun, Jirsa, & Lindenberger, 2016;Olejarczyk, Marzetti, Pizzella, & Zappasodi, 2017;Stankovski, Ticcinelli, McClintock, & Stefanovska, 2017;Scally, Burke, Bunce, & Delvenne, 2018). While these studies have revealed much in terms of spatial distributions and graph-theoretic metrics of amplitude and phase relations, their primary goals were to target specific frequency bands to identify neural correlates of individual differences, dysfunctions, or tasks, or to distinguish spatial networks mediated by different frequency bands. The goal of the current study was complementary in that we aimed to elucidate general rules that may govern large-scale intrinsic neural dynamics by directly examining EEG spatiotemporal dynamics as multidimensional state vectors (without applying time-frequency decomposition).
We coded EEG from k sites as a k-dimensional state vector ! " ####⃗ (&) ( Figure 1A). Global EEG dynamics are entirely captured by the k-dimensional trajectories of a state vector. Large-scale neural processes are likely reflected in how these trajectories evolve. For example, a persistent trajectory movement in the same direction would indicate a continuation of the same process. The degree of directional instability of a trajectory (as a function of time) for a specific timescale Dt can be defined as the angular change q in trajectory direction  The gray curves illustrate state trajectories with black dots indicating temporal sampling (512 per second). The red arrows illustrate Δ! " # − ∆#/2 and Δ! " # + ∆#/2 , each linearly interpolated over Dt, used to compute trajectory instability at time t on timescale Dt in terms of the angle between them, q (roughly curvature on timescale Dt). Left. A smoothly curved trajectory is characterized by smaller q values for a shorter Dt (e.g., Dt = 2 timepoints) and larger q values for a longer Dt (e.g., Dt = 4 timepoints). Right. A rough but overall straight trajectory is characterized by larger q values for a shorter Dt (e.g., Dt = 2 timepoints) and smaller q values for a long Dt (e.g., Dt = 4 timepoints).
Our goal was to identify general rules governing intrinsic neural dynamics. Our strategy was to identify dynamical properties of EEG-state trajectories that are globally preserved. Specifically, we divided the "whole brain" (those aspects reflected in the current sources estimated from scalp EEG) into several regions that likely perform distinct computations in order to identify dynamical properties that remained invariant across those functionally distinct regions. Further, region specific differences in dynamical properties may elucidate how global rules may change to accommodate different functional demands.
We thus divided the 64 sites into five regions, posterior, central, anterior, left-lateral and right-lateral regions (see Figure 2) because prior research suggests that these regions may subserve broadly distinct processes that may require different dynamics (e.g., Laufs et al., 2003;Jann et al., 2010;Luck & Kappenman, 2012;Pratt, 2011;Smulders, Miller, & Luck, 2012;Swaab, Ledoux, Camblin, & Boudewyn, 2012;Singh, Bauer, Chowanski, Sui, Atkinson, Baurley, Fry, Evans, & Bianchi-Berthouze, 2014). Overall, prior research suggests that EEG activity in the posterior region primarily reflects visual processes (including attention and working memory), activity in the central region primarily reflects auditory, tactile, linguistic, mnemonic, and motor processes (also including attention and working memory), and activity in the anterior region primarily reflects cognitive, regulatory, attentional, and decisional processes, while activity in any of these regions may be lateralized. Further, our recent study has demonstrated that these regions engage in distinct patterns of crossfrequency amplitude-amplitude and phase-amplitude interactions (Menceloglu, Grabowecky, & Suzuki, 2020).
We thus applied the trajectory-instability and trajectory-speed analyses to the resting-state EEG recorded from these regions to identify the dynamical properties of EEG-state trajectories that remained invariant across all regions and those that were regionally altered to accommodate specific functional demands.

Participants
Forty-eight Northwestern University students (33 women, 14 men, and 1 non-binary; ages 18 to 29 years, M = 20.83, SD = 2.65) gave informed consent to participate for monetary compensation ($10 per hour). All were right-handed, had normal hearing and normal or corrected-to-normal vision, and had no history of concussion. They were tested individually in a dimly lit room. The study protocol was approved by the Northwestern University Institutional Review Board.

EEG recording procedures and preprocessing
EEG signals were recorded from 64 scalp electrodes at a sampling rate of 512 Hz using a BioSemi ActiveTwo system (using the default filter settings; attenuation to -3 dB at 1/5 of the sampling rate or 102.4 Hz; see www.biosemi.com) for about 5 minutes. The participants were instructed to rest with their eyes closed while minimizing movements and freely engage in whatever thoughts that spontaneously arose. We expect that the degree of mind wandering varied from participant to participant. This may have contributed some variability in our results given that resting fMRI-BOLD activity has been shown to vary when mind wandering is encouraged versus discouraged (Kawagoe, Onoda, & Yamaguchi, 2018).
Two additional electrodes were placed on the left and right mastoid area. The EEG data were preprocessed with MATLAB using the EEGLAB and ERPLAB toolboxes (Delorme & Makeig, 2004;Lopez-Calderon & Luck, 2014). The data were re-referenced offline to the average of the two mastoid electrodes, bandpass-filtered at 0.01-80 Hz, and notch-filtered at 60 Hz (to remove power-line noise that affected the EEG data from some of the participants). Based on our inspection of the EEG data, we did not identify any reasons to remove EEG signals over any time segments from any electrodes for any participants.
Because we compared EEG-state trajectory dynamics across regions, it was important to localize the underlying neural current sources from scalp-recorded EEG, especially to minimize any influences of volume conduction. To this end, we used a surface-Laplacian transform which has been shown to appropriately approximate the underlying current sources (especially the superficial sources near the skull) in a data-driven manner with minimal assumptions about the underlying sources (e.g., Hjorth, 1980;Kayser and Tenke, 2006;Tenke and Kayser, 2012). This transform also essentially removes effects of reference electrode choices and substantially reduces volume conduction effects largely confining them within adjacent electrodes for a 64electrode montage like that used in the current study (e.g., Cohen, 2014). We surface-Laplacian transformed the EEG data using Perrin and colleagues' method (e.g., Perrin, Pernier, Bertrand, Giard, & Echallier, 1987;Perrin, Pernier, Bertrand, & Echallier, 1989a;1989b) with a typical set of parameter values (Cohen, 2014).
Here we refer to surface-Laplacian transformed EEG simply as EEG.

EEG data analysis
The continuous ~5 min EEG data (excluding the first 10 sec) were divided into consecutive 5-sec time segments and the first 58 segments from all participants were used for the analyses (some participants yielded 59 segments because their recording sessions were slightly longer). We chose to analyze the data in 5-sec segments for several reasons. First, as we focused on neural dynamics on the sub-second timescale, the segmenting allowed us to generate phase-scrambled control data (see below) for each 5-sec segment to control for any effects of gradual (slower than ~0.2 Hz) changes in spectral amplitudes. The choice of 5 sec was partly motivated by our recent finding that a half cycle of several seconds is a characteristic rate at which EEG spectral amplitudes fluctuate (Menceloglu et al., 2020). Second, we wanted to reliably characterize trajectory instability and speed up to the timescale of about a quarter of a second which corresponded to about a half-second cycle for the relevant underlying oscillatory activity (see below); thus, the use of 5-sec segments allowed us to include ~10 cycles across all timescales examined. We note in passing that the dynamical properties identified in the current study were relatively stable over a ~5 min period (see below), implying that the use of a longer analysis segment would have yielded comparable results.
While phase-scrambling can be performed using several different methods, we chose discrete cosine transform, DCT (e.g., Kiya et al., 2010). In short, we transformed each 5-sec EEG waveform with type-2 DCT, randomly shuffled the signs of the coefficients, and then inverse-transformed it with type-3 DCT, which yielded a phase-scrambled version. DCT-based phase-scrambling is similar to DFT (discrete Fourier transform) -based phase-scrambling except that it is less susceptible to edge effects especially when phase-scrambling needs to be performed on relatively short waveforms.
The spectral amplitudes from the five regions showed typical profiles. Fast Fourier transform (FFT) was applied to each 5-sec EEG segment and the results were averaged across all segments per site per participant; the values were then averaged across sites within each region and averaged across participants. The resultant regional time-averaged spectral-amplitude profiles are shown in Figure 2. All regional profiles included the typical 1/f b component that partly reflects an Ornstein-Uhlenbeck process which is a component of the standard integrate-and-fire model of a neuron, accounting for the stochastic membrane-potential dynamics in the absence of input currents (e.g., Koch, 1999;Mazzoni et al., 2008;see He, 2014, for a review of the various contributing factors). As expected from EEG recordings with the eyes closed, the posterior spectralamplitude profiles included a pronounced a (~10 Hz) peak (the black curve in Figure 2) that was reduced in the central and left/right lateral profiles (the blue, orange, and green curves) and virtually absent in the anterior profiles (the red curve).
We represented the global spatiotemporal dynamics of each region as a 17-dimensional (k = 17) state ∆+^ -Eq. 4. Both trajectory instability and speed may depend on phase-independent spectral-amplitude profiles (e.g., reduced high-frequency amplitudes would generate smoother trajectories regardless of phase relations) as well as on phase relations. Comparisons with phase-scrambled controls (see above) allowed us to isolate the phase-dependent aspects of trajectory instability and speed indicative of active phase coordination. In particular, by using two types of phase scrambling (see Results and Discussion), we distinguished the contributions of cross-site within-frequency and cross-frequency phase relations to EEG-state trajectories. We examined the temporal average and standard deviation of trajectory instability and speed, as well as the temporal correlation between trajectory instability and speed, computed for both real data and the corresponding phase-scrambled controls as a function of Dt. These statistics were first computed for each 5sec segment and then averaged across all segments per region per participant. We demonstrate statistical reliability with small error regions (indicating ±1 standard error of the mean with participants as the random effect) relative to mean differences, virtually identical patterns obtained from the odd and even numbered participants, and clearly separated distributions of values.

Effects of spectral-amplitude profiles on trajectory instability and speed
The timescale dependence of EEG state-trajectory instability was characteristically different across the posterior, central, anterior, left-lateral, and right-lateral regions ( Figure 3A). On the shorter timescales (Dt < ~25  Figure 3B) that minimized the inter-region standard deviation in trajectory instability, (2) by the narrow distribution of the trajectoryinstability levels (mostly within 90°-110°; Figure 3C) at the critical timescale, and (3) by the high precision of convergence, with the inter-region standard deviations of trajectory instability dipping to only 1°-2° at the critical timescale for most participants ( Figure 3D). The virtually identical patterns of timescale dependence of trajectory instability for the odd and even numbered participants ( Figure 3E) demonstrate statistical reliability.
We note that the tight global (inter-region) convergence of state-trajectory instability on the timescale of ~25 ms is primarily a characteristic of baseline spectral-amplitude profiles independent of phase relations (see below).
State-trajectory speed overall diminished with increasing timescale Dt in all five regions ( Figure 3F). In addition, the posterior region had the overall fastest trajectories, followed by the left/right lateral regions and the anterior region with the central region having the slowest trajectories. Further, the trajectory speeds in the posterior, left/right lateral, and anterior regions converged on the short (Dt < ~10 ms) and long (Dt > ~130 ms) timescales, while the trajectory speeds in the central region remained low across all examined timescales.
These patterns are reliable as those obtained from the odd and even numbered participants were virtually identical ( Figure 3G). We note that these patterns of state-trajectory speed are also primarily a characteristic of baseline spectral-amplitude profiles independent of phase relations (see below). Notably, the five curves tightly converge at Dt ~ 25 ms where the level of trajectory instability is ~100°(vertical and horizontal dashed lines). To illustrate the consistency of this convergence across participants, the cross-participant distributions (histograms) of the Dt values and trajectory-instability levels at the convergence are shown in B and C, respectively. To illustrate the precision of convergence, the cross-participant distribution of the inter-region standard deviation (sd) of trajectory instability at the convergence is shown in D (mostly 1°-2°and all less than 5°). E. The time-averaged trajectory instability pattern in A shown separately for the odd and even numbered participants; the virtually identical patterns demonstrate statistical reliability. F. Time-averaged trajectory speed for the posterior (black), central (blue), anterior (red), left-lateral (orange), and right-lateral (green) regions as a function of Dt. G. The time-averaged trajectory speed in F shown separately for the odd and even numbered participants; the virtually identical patterns demonstrate statistical reliability. The error bars represent ±1 standard error of the mean.

Effects of phase relations on trajectory instability and speed
Comparing the real trajectory instability and speed with those computed using the corresponding phasescrambled controls revealed trajectory characteristics that depended on spectral-phase relations over and above baseline spectral amplitudes (our phase-scrambled controls closely tracked slow variations in spectral amplitudes as we used 5-sec segments as the unit of analyses; see Material and Methods). The phasedependent variability in time-averaged trajectory instability was small for all regions, limited to less than ±2° ( Figure 4A), indicating that the characteristic regional profiles of time-averaged trajectory instability including the tight convergence at Dt ~ 25 ms ( Figure 3A) primarily reflect baseline spectral-amplitude profiles.
Nevertheless, phase relations produced some consistent effects on trajectory instability. The posterior region yielded increased instability (relative to the phase-scrambled control) on the short < 15 ms and longer ~100-140 ms timescales (the black curve in Figure 4A) with the latter peak accompanied by the left and right lateral regions (the orange and green curves). Further, all regions except the anterior region (red) yielded reduced instability on the mid-range ~30-80 ms timescales, all regions except the posterior region (black) yielded reduced instability on the short < 8 ms timescales, and all regions yielded reduced instability on the long > ~180 ms timescales. These patterns were obtained from both the odd and even numbered participants (the right panels in Figure 4A).
While phase relations had overall small effects on time-averaged trajectory instability ( Figure 4A), they substantially reduced time-averaged trajectory speed in all regions. Further, the levels of speed reductions in dB were largely invariant across all examined timescales in all regions ( Figure 4B). This suggests that similar mechanisms of phase coordination divisively reduce trajectory speed in all regions. This also indicates that the characteristic regional profiles of time-averaged trajectory speed ( Figure 3F) primarily reflect baseline spectralamplitude profiles.
Active control of phase relations, which includes synchronization and desynchronization, likely influences temporal variability in trajectory instability and speed. Accordingly, the temporal variability (measured as temporal standard deviation) in both trajectory instability ( Figure 4C) and speed ( Figure 4D) were substantially increased in all regions (relative to the phase-scrambled controls). Further, the substantial increases in variability (in degrees for the instability variability and in dB for the speed variability) were observed across all examined timescales in all regions. Nevertheless, the variability in both trajectory instability and speed was most elevated in the posterior region, followed by the anterior and left/right lateral regions and least elevated in the central region ( Figure 4C-D). These patterns were virtually identical for the odd and even numbered participants (the right panels in Figure 4C-D).
Given this evidence of coordination of trajectory instability and speed through modulation of phase relations, they may be jointly coordinated. Specifically, a positive temporal correlation between trajectory instability and speed would indicate faster trajectories being less stable (and slower trajectories being more stable), a negative correlation would indicate faster trajectories being more stable (and slower trajectories being less stable), whereas a lack of correlation would suggest independent control of trajectory instability and speed. We estimated the temporal correlations between trajectory instability and speed that reflected active spectral-phase coordination as the Fisher-Z transformed Spearman's r values (robust against outliers) computed for the real data minus those computed for the phase-scrambled controls. Across all five regions trajectory instability and speed were overall either positively correlated or uncorrelated, implying that the dynamics in all regions are coordinated such that faster trajectories are less stable (and slower trajectories are more stable) ( Figure 4E). In addition, all regions except the anterior region (red) yielded positive correlations on the ~20-80 ms timescales, all regions except the posterior region (black) yielded positive correlations on the fast < 8 ms timescales, and all regions yielded positive correlations on the slow > ~180 ms timescales. The lack of correlation for the anterior region (red) on the timescales of ~15-100 ms is notable, suggesting that trajectory instability and speed are independently controlled in the anterior region on these timescales. All these patterns were virtually identical for the odd and even numbered participants (the right panels in Figure   4E). For each plot, the reliability of the data patterns is demonstrated by the accompanying plots on the right showing virtually identical patterns obtained from the odd and even numbered participants. A. Time-averaged trajectory instability as a function of Dt (real data minus the corresponding phase-scrambled controls). B. Time-averaged trajectory speed as a function of Dt (real data vs. the corresponding phase-scrambled controls in dB). C. The temporal variability (measured as standard deviation) in trajectory instability as a function of Dt (real data minus the corresponding phase-scrambled controls). D. The temporal variability (measured as standard deviation) in trajectory speed as a function of Dt (real data vs. the corresponding phase-scrambled controls in dB). E. The temporal correlation between trajectory instability and speed as a function of Dt (real data minus the corresponding phase-scrambled controls in z-transformed Spearman r). In all panels, the data for the posterior, central, anterior, left-lateral, and right-lateral regions are shown in black, blue, red, orange, and green curves, respectively (see the region map). The error bars represent 1 standard error of the mean.

Effects of within-frequency versus cross-frequency phase relations
Phase relations include both within-frequency phase relations-phase relations within the same frequency bands across sites-and cross-frequency phase relations-phase relations between different frequency bands within and across sites. The analyses presented in Figure 4 include the contributions of both types of phase relations because the phase-scrambling procedure randomized spectral phase independently for each frequency component at each site. We examined the unique contributions of cross-frequency phase relations using "cross-frequency phase-scrambled" controls, computed by independently randomizing phase for each frequency component but applying the same set of cross-frequency phase shifts to all sites. This procedure randomized cross-frequency phase relations within and across sites (e.g., X-Hz and Y-Hz components were phase-shifted by Dx and Dy degrees), but preserved within-frequency phase relations across sites (e.g., X-Hz components at all sites were phase-shifted by Dx degrees and Y-Hz components at all sites were phase-shifted by Dy degrees).
The effects of cross-frequency phase relations alone (measured relative to the cross-frequency phasescrambled controls) were generally smaller than the combined effects of cross-and within-frequency phase relations (measured relative to the full phase-scrambled controls) ( Figure 5 vs. Figure 4). Nevertheless, some results are informative. Note that comparisons with the full phase-scrambled controls reflect combined influences of both cross-site within-frequency and cross-frequency phase relations, whereas comparisons with the cross-frequency phase-scrambled controls reflect influences of only cross-frequency phase relations. Thus, (1) a unique contribution of cross-site within-frequency phase relations may be inferred based on a reliable difference from the full phase-scrambled control with no reliable difference from the cross-frequency phasescrambled control. (2) A unique contribution of cross-frequency phase relations may be inferred based on reliable but equivalent differences from both the full and cross-frequency phase-scrambled controls. (3) Joint contributions of cross-site within-frequency and cross-frequency phase relations may be inferred based on a large reliable difference from the full phase-scrambled control and a substantially smaller but still reliable difference from the cross-frequency phase-scrambled control. We used this logic to interpret the results relative to the full and cross-frequency phase-scrambled controls.
For time-averaged trajectory instability, cross-site within-frequency phase relations uniquely contributed to increasing trajectory instability in the posterior region on the fast < 15 ms timescales as it was elevated relative to the full phase-scrambled control (the black curve in Figure 4A) but equivalent to the cross-frequency phase-scrambled control (the black curve in Figure 5A). Cross-frequency phase relations uniquely contributed to reducing trajectory instability in all regions except the posterior region on the fast < 8 ms timescales as the reductions were equivalent relative to both the full and cross-frequency phase-scrambled controls (the nonblack curves in Figures 4A and 5A). Similarly, cross-frequency phase relations uniquely contributed to reducing trajectory instability in all regions on the slow > ~180 ms timescales as the reductions were equivalent relative to both the full and cross-frequency phase-scrambled controls (all curves in Figures 4A and 5A). All these patterns were virtually identical for the odd and even numbered participants (the right panels in Figures 4A and   5A).
For time-averaged trajectory speed, both cross-site within-frequency phase relations and cross-frequency phase relations contributed to its global (across all regions) and largely timescale invariant reductions as the reductions were consistently larger relative to the full phase-scrambled controls than relative to the crossfrequency phase-scrambled controls (compare Figures 4B and 5B), but the smaller reductions relative to the cross-frequency phase-scrambled controls were still consistent ( Figure 5B). These patterns were virtually identical for the odd and even numbered participants (the right panels in Figures 4B and 5B).
For the temporal variability in trajectory instability, cross-site within-frequency phase relations uniquely contributed to its global (across all regions) and largely timescale invariant increases as the increases were substantial relative to the full phase-scrambled controls ( Figure 4C) but virtually absent relative to the crossfrequency phase-scrambled controls ( Figure 5C); if anything, cross-frequency phase relations contributed to reducing the temporal variability in trajectory instability in the posterior region on the ~15-40 ms and ~60-90 ms timescales (the black curve in Figure 5C). These patterns were virtually identical for the odd and even numbered participants (the right panels in Figures 4C and 5C).
For the temporal variability in trajectory speed, both cross-site within-frequency phase relations and cross-frequency phase relations contributed to its global (across all regions) and largely timescale invariant increases as the increases were consistently larger relative to the full phase-scrambled controls than relative to the cross-frequency phase-scrambled controls (compare Figures 4D and 5D), but the smaller increases relative to the cross-frequency phase-scrambled controls were still consistent ( Figure 5D). These patterns were virtually identical for the odd and even numbered participants (the right panels in Figures 4D and 5D).
For the temporal correlation between trajectory instability and speed, cross-frequency phase relations uniquely contributed to the patterns of positive correlations in the left/right-lateral and anterior regions across all examined timescales and in the central and posterior regions on the timescales slower than ~20 ms as those patterns were equivalent relative to both the full and cross-frequency phase-scrambled controls ( Figures   4E and 5E); an exception may be that in the posterior region on the timescales of ~30-70 ms, the positive correlations were more elevated relative to the full phase-scrambled control (the black curve in Figure 4E) than relative to the cross-frequency phase-scrambled control (the black curve in Figure 5E) suggesting an additional contribution of cross-site within-frequency phase relations. In addition, in the posterior and central regions on the fast < ~15 ms timescales, cross-frequency phase relations increased but cross-site within-frequency phase relations reduced the positive correlations as they were more elevated relative to the cross-frequency phasescrambled controls (the black and blue curves in Figure 5E) than relative to the full phase-scrambled controls (the black and blue curves in Figure 4E). It is noteworthy that phase relations appear to be coordinated in such a way that trajectory instability and speed are uniquely independent in the anterior region on the timescales of ~15-120 ms as the correlations did not deviate from the full phase-scrambled control (the red curve in Figure   4E). All these patterns were virtually identical for the odd and even numbered participants (the right panels in Figures 4E and 5E). The cross-frequency-phase-relation dependent aspects of trajectory instability and speed and their temporal correlation as a function of timescale Dt relative to the corresponding cross-frequency phase-scrambled controls (that preserved cross-site withinfrequency phase relations). For each plot, the reliability of the data patterns is demonstrated by the accompanying plots on the right showing virtually identical patterns obtained from the odd and even numbered participants. A. Time-averaged trajectory instability as a function of Dt (real data minus the corresponding phase-scrambled controls). B. Time-averaged trajectory speed as a function of Dt (real data vs. the corresponding phase-scrambled controls in dB). C. The temporal variability (measured as standard deviation) in trajectory instability as a function of Dt (real data minus the corresponding phase-scrambled controls). D. The temporal variability (measured as standard deviation) in trajectory speed as a function of Dt (real data vs. the corresponding phase-scrambled controls in dB). E. The temporal correlation between trajectory instability and speed as a function of Dt (real data minus the corresponding phase-scrambled controls in z-transformed Spearman r). In all panels, the data for the posterior, central, anterior, left-lateral, and right-lateral regions are shown in black, blue, red, orange, and green curves, respectively (see the region map). The error bars represent 1 standard error of the mean.

Menceloglu et al. 20
Finally, we evaluated the temporal stability of the results presented in Figures 3-5 by computing the corresponding quantities separately in successive 50-sec intervals, dividing each ~5-min EEG recording period into 5 successive intervals. The data from the 5 intervals are shown as superimposed curves in Figure 6 with later intervals indicated with lighter colors. The time-averaged trajectory instability and speed as a function of Dt (presented in Figure 3A and F), which primarily depend on phase-independent spectral-amplitude profiles ( Figure 2), are replotted in Figure 6A-B. The full-phase-relation dependent aspects (relative to the full phasescrambled controls) of the time-averaged trajectory instability and speed, their temporal variability, and their temporal correlations as functions of Dt (presented in Figure 4A-E) are replotted in Figure 6C-G. The crossfrequency-phase-relation dependent aspects (relative to the cross-frequency phase-scrambled controls) of these variables (presented in Figure 5A-E) are replotted in Figure 6H-L. It is clear from Figure 6 that all results are reasonably stable over a ~5-min period. Further, all the results are remarkably consistent for the odd and even numbered participants.

General Discussion
We investigated potential global rules governing large-scale intrinsic neural dynamics by characterizing the spatiotemporal dynamics of resting-state EEG as trajectories of k-dimensional state vectors where k is the number of scalp sites (with each site corresponding to a current source estimated with the surface-Laplacian transform). We focused on the timescale dependence of state-trajectory instability and speed because these variables, approximately the 2 nd and 1 st temporal derivatives, are informative for characterizing general properties of state trajectories (e.g., Figure 1). We examined these trajectory properties in five broadly defined regions (posterior, central, frontal, left-lateral, and right-lateral regions) that likely perform distinct computations based on prior research including ours (e.g., Laufs et al., 2003;Jann et al., 2010;Luck & Kappenman, 2012;Pratt, 2011;Smulders et al., 2012;Swaab et al., 2012;Menceloglu et al., 2020a). Our strategy was to characterize the global rules governing intrinsic neural dynamics by identifying the trajectory properties that are preserved across all these regions. A secondary goal was to examine how those rules may be altered in separate regions to shed light on how global rules may change to support regional computational needs.
To interpret the current results in the form of timescale dependencies of EEG-state trajectory instability and speed, we consider how behaviors of these state variables relate to underlying neural spectral dynamics.
Neurophysiological signals including scalp-recorded EEG have been known to exhibit characteristic spectralamplitude profiles, generally including a 1/f b background profile with embedded peaks at specific frequency bands such as q, a, and b bands (e.g., see van Albada and Robinson, 2013, for a review; also see Figure 2 for examples from the current study). The 1/f b background is thought to reflect random-walk type neuronal membrane fluctuations (in the absence of input currents) described by an Ornstein-Uhlenbeck process (Eq. 5) included in the standard integrate-and-fire model (e.g., Koch, 1999;Mazzoni et al., 2008), where v(t) is the membrane potential, g is 1/tm, the reciprocal of the membrane time constant, D is the noise intensity multiplied by g 2 , and gn(t) is Gaussian noise. The membrane potential, v(t), has a monotonically decreasing spectral-power profile, reducing to half power (or 1 √2 ⁄ amplitude) at the frequency, h i/? = j ?k , and approximating a 1/f 2 power profile (or a 1/f amplitude profile) at higher frequencies (see He, 2014, for a review of contributions of this and other mechanisms to the 1/f b spectral-power profiles). Peaks at specific frequency bands are thought to reflect oscillatory neural activity, and their roles in coordinating processes that mediate perceptual, attentional, memory, and cognitive operations have been well studied (e.g., van Beijsterveldt and van Baal, 2002;Fries, 2005;Palva & Palva, 2007;Duzel et al., 2010;Engel & Fries, 2010;Dahaene & Changeux, 2011;Donner & Siegel, 2011;Palva & Palva, 2012;Klimesch, 2012;Jensen et al., 2014;Michalareas et al., 2016;Scally et al., 2018). However, it has been less clear how the diverse spectral-amplitude profiles observed in different brain regions may be coordinated to maintain the global integrity of neural dynamics. The current results may contribute to an understanding of dynamical parameters relevant to such global coordination.

Influences of spectral-amplitude profiles on EEG state-trajectory instability and speed
The time-averaged state-trajectory instability in the posterior, central, anterior, left-lateral, and right-lateral regions tightly converged on the timescale of ~25 ms at the instability level of ~100° (Figures 3A-D). Because phase relations made minimal contributions to the time-averaged trajectory instability (less than ±2° across all examined timescales in all regions; Figure 4A), this convergence likely reflects a global coordination of baseline spectral-amplitude profiles. That is, the baseline spectral-amplitude profiles appear to be globally (across brain regions) coordinated in such a way that, while state trajectories in different regions make different average degrees of turns on most timescales, on a universal timescale of ~25 ms, all state trajectories uniformly make ~100° turns on average, that are near orthogonal (i.e., near 90°) but slightly conservative. This implies that successive state-trajectory changes on the timescale of ~25 ms are globally coordinated to be near-maximally unrelated, that is, near-minimally constrained, but slightly restoring.
How might time-averaged trajectory instability be globally maintained at ~100° on the timescale of ~25 ms? The timescale dependence of time-averaged trajectory-instability for each region depends on the exact spectral-amplitude profiles at its constituent sites. Nevertheless, there are some factors that make systematic contributions. For example, a Gaussian-random-walk type of EEG activity occurring at all constituent sites would yield a constant level of time-averaged trajectory instability at 90° regardless of timescale (making memory free orthogonal turns at all timescales on average).
Oscillatory activities occurring at any of the constituent sites would also make characteristic contributions.
For example, if EEG at most sites oscillated at the same frequency, the state trajectory would move along an approximate elliptical path. The elliptical path would be planar only if EEG at all sites oscillated at the same frequency (with a non-singular distribution of phase). Presence of oscillations in other frequencies and the 1/f b spectral background would distort the elliptical path in a variety of ways. Nevertheless, if oscillations at a specific frequency were dominant across sites relative to other contributions, the path would locally approximate a curved helix. In that case, average trajectory instability would be near 180° on the timescale of a half period (i.e., approximately reversing directions from half period to half period) and near 0° on the timescale of a period (i.e., approximately moving in the same direction from period to period). Oscillatory directional changes would again optimally impact the state trajectory on the timescale of 1.5 periods, and minimally impact it on the timescale of 2 periods, and so on, peaking at the odd multiples and dipping at the even multiples of the oscillatory half period. For real EEG data that include a broad range of oscillatory frequencies, the impact of a specific oscillation frequency would diminish for the higher harmonics because the directional instability on longer timescales would be more strongly influenced by slower oscillatory components and slower-varying non-oscillatory components (including the random-walk type activity that partially underlies the 1/f b spectral profile).
An example of an oscillatory contribution can be seen in the current results for the posterior region whose constituent sites yielded a substantial spectral peak in the a-band (~10 Hz with 50 ms half-period) (the black curves in Figure 2). Based on our reasoning above, the a-band contribution to the time-averaged trajectory instability should generate diminishing peaks on the timescales of 50 ms (a half-period), 150 ms (3 halfperiods), etc., and diminishing troughs on the timescales of 100 ms (2 half-periods), 200 ms (4 half-periods), etc. This is what we observed (the black curve in Figure 3A). However, note that even the largest peak and trough on the half-period and one-period timescales only reached ~130° and ~75°, respectively (rather than  Figure 3A), consistent with the fact that the a-band peaks were less pronounced in the central and left/right-lateral regions (the blue, orange, and green curves in Figure 2) relative to the posterior region. The anterior instability profile was relatively flat on the timescales beyond ~15 ms, consistent with the absence of any substantial spectral peaks (the red curve in Figure 2). Nevertheless, this relatively flat portion reflects oscillatory activities in a broad range of frequencies rather than a dominance by the random-walk type spectral background because the plateaued instability levels over ~15-70 ms timescales were consistently above the 90° level expected by random walk ( Figure 3A). The trajectory-instability profiles from all regions gradually converged toward the level of 90° on the longer timescales > ~150 ms ( Figure 3A), indicating that state trajectories on average make relatively unconstrained directional changes (with successive trajectory changes being maximally unrelated) on the timescales beyond ~150 ms.
Given that trajectory instability profiles are sensitive to oscillatory activities and 1/f b spectral baselines at constituent sites, the fact that the profiles from all five regions tightly converged at the level of ~100° on the timescale of ~25 ms ( Figure 3A) may have functional relevance rather than being epiphenomenal. The convergence may suggest, for example, that intrinsic neural dynamics are globally calibrated in such a way that bottom-up sensory and top-down attentional/cognitive processes could exchange information in ~25 ms packets at a matched level of nearly unconstrained but slightly conservative (~100°) state-trajectory flexibility.
This matching may facilitate directionally balanced interfacing of bottom-up and top-down processes so that they contribute on an equal footing to construct internal representations that optimally incorporate both behavioral goals and environmental constraints. Because b-band activity (~20 Hz, ~25 ms half-period) would have the largest impact on state-trajectory instability on the convergent timescale of ~25 ms, an efficient way to maintain this calibration may be to network-wide tune b-band activity to compensate for any cross-network variations in multi-frequency oscillatory and non-oscillatory activity.
While these inferences are highly speculative, the data provide some indirect support for the interpretation that the inter-region convergence of time-averaged trajectory instability at ~100° at ~25 ms may be actively maintained. If the convergence were actively maintained, time-averaged trajectory instability should be particularly stable on the critical timescale of ~25 ms in all five regions, both across participants and over time. Indeed, the cross-participant standard deviation of trajectory instability was minimized at the critical ~25 ms timescale in all five regions ( Figure 7B) and the results were virtually identical for the odd and even numbered participants (the right panels in Figure 7B). To evaluate the temporal stability of time-averaged trajectory instability, we computed average trajectory instability over each consecutive 5-sec time segment, then evaluated its temporal variability as the standard deviation across all segments-inter-time-segment sd.
This measure of temporal variability also minimized at the critical ~25 ms timescale in all five regions (though shifted to ~20 ms in the anterior region) ( Figure 7C); the results were again virtually identical for the odd and even numbered participants (the right panels in Figure 7C). The histogram of the variability-minimizing timescale values (derived from the inter-time-segment sd profiles averaged across all regions) was sharply peaked at the critical ~25 ms timescale ( Figure 7D), indicating that time-averaged trajectory instability was most temporally stable on the critical ~25 ms timescale for most participants.
Note that perturbations in any dimensions orthogonal to the plane of a trajectory turn would shift the degree of the turn toward 90°, and the shift would be larger for trajectory turns at angles more deviated from 90° (with no shirts for 90° turns). Thus, random perturbations would generate greater variability when average turns (reflected in time-averaged trajectory instability) are at angles more deviated from 90°. One can confirm this prediction for timescales longer than the critical ~25 ms timescale. For all regions, when the time-averaged trajectory instability substantially deviated from 90° (e.g., on ~50 ms and ~100 ms timescales for the posterior, central, and left/right-lateral regions; Figure 7A), both the cross-participant and inter-time-segment sd's tended to be at local maxima ( Figure 7B-C), whereas when the time-averaged trajectory instability crossed (or approached) the 90° level (e.g., on ~75 ms and ~160 ms timescales for the posterior region; Figure 7A), both the cross-participant and inter-time-segment sd's tended to be at local minima ( Figure 7B-C). Remarkably, despite the fact that the time-averaged trajectory instability was substantially higher than 90° on the critical ~25 ms timescale ( Figure 7A  In contrast to the feature-rich profiles of time-averaged trajectory instability, time-averaged trajectory speed largely yielded monotonically decreasing profiles as a function of timescale for all five regions ( Figure   3F). Gaussian-random-walk type EEG activity occurring at all sites would yield a power-function profile of timeaveraged trajectory speed, which would be a negatively sloped line in the log-log plot, with the intercept proportional to the temporal standard deviation of EEG. An oscillatory activity would contribute maximal directional changes twice per oscillation period. Thus, on timescales sufficiently shorter than a half-period, the path segments linearly interpolating over Dt (used in the current analyses; Figure 1)  about a half-period, thereby making the profile upwardly bowed with a broad peak trailing the timescale of a half period. This reasoning is consistent with the fact that the posterior state trajectory with its prominent aband peak (the black curve in Figure 2) yielded a substantially bowed trajectory-speed profile with the curvature broadly peaking on the timescale a little longer than ~50 ms (the black curve in Figure 3F). The degree of bowing is less for the central and left/right-lateral regions (the blue, orange, and green curves in Figure 3F) with less prominent a-band peaks (the blue, orange, and green curves in Figure 2). The anterior region with no particularly prominent spectral peaks yielded a relatively linear (in log-log) time-averaged trajectory-speed profile (the red curve in Figure 3F).
We note that the notch filter at 60 Hz used to remove line noise (see Materials and Methods) attenuated oscillatory components that would have strongly influenced trajectory instability and speed on ~10 ms and shorter timescales (note that the half period of 60 Hz is ~8 ms). The filter thus caused the general flattening of the time-averaged trajectory-speed profiles ( Figure 3E) as well as the general falling of the time-averaged trajectory-instability profiles ( Figure 3A) on ~10 ms and shorter timescales.

Influences of phase relations on EEG state-trajectory instability and speed
While the discussion so far has focused on the influences of spectral-amplitude profiles on EEG-state trajectory instability and speed, our results suggest that they are also influenced by phase relations. We used two types of phase scrambling to infer the contributions of cross-site within-frequency and cross-frequency phase relations. Comparisons with the full-phase-scrambled controls (EEG from each site independently phase-scrambled) elucidated the joint contributions of cross-site within-frequency and cross-frequency phase relations to trajectory instability and speed, whereas comparisons with the cross-frequency phase-scrambled controls (EEG from all sites phase-scrambled with the same randomization seed thus preserving cross-site within-frequency phase relations) isolated the contributions of cross-frequency phase relations.

Influences of phase relations on trajectory instability
We first consider the potential impact of cross-site within-frequency phase relations on state-trajectory instability, and compare them with the results. At one extreme, if oscillatory activities in a specific frequency band were all phase aligned (i.e., 0° phase lagged) across sites, they would contribute an approximately oscillatory component to a state trajectory. At another extreme, if their phases were minimally aligned across sites, that is, if the phases were uniformly distributed between 0° and 180° across sites, the oscillatory activities would contribute an approximately circular component to a state trajectory. Strictly speaking, the contribution would be circular only if oscillation amplitudes were equal across sites. Uneven amplitudes would make the contribution elliptical. Given a set of uneven amplitudes, the contribution would have the smallest aspect ratio (i.e., closest to being circular) when the phases of large-amplitude oscillations were evenly distributed between 0° and 180°. An intermediate level of cross-site phase alignment would contribute an approximately elliptical component with a higher (narrower) aspect ratio when the phases are more tightly aligned across sites. Note that phase scrambling would reduce the aspect ratio of these elliptical contributions but would not make them circular because a randomized phase distribution would necessarily deviate from a uniform distribution. As expected, when we phase scrambled artificial data consisting of single-frequency oscillatory sources with evenly distributed phases across sites, the originally circular trajectory became elliptical.
Because trajectory instability is approximately a timescale dependent measure of trajectory curvature, a narrow elliptical component would generate large temporal variability within each oscillatory cycle because directional changes would be large over the two high-curvature portions and small over the two low-curvature portions. Because each of the high-and low-curvature portions spans about a half cycle, an elliptical component would particularly elevate trajectory-instability variability measured on the timescale of a quarter period (sharply changing directions over the half period containing the high curvature and largely moving in the same direction over the half period containing the low curvature). This cyclic variation in trajectory instability would become less pronounced for elliptical components with lower aspect ratios generated by less tight crosssite phase alignment because the high and low curvatures would be less different. Thus, a tighter cross-site within-frequency phase alignment of oscillatory activities should increase trajectory-instability variability especially on the quarter-period timescale. In contrast, time-averaged trajectory instability is not expected to be much affected by cross-site within-frequency phase alignment because the average curvature of an elliptical path would be relatively independent of its aspect ratio. We confirmed this reasoning by simulating the waveform at each site as a sinewave oscillating at f Hz and varying the tightness of the phase distribution across 17 sites (because our EEG analyses included 17 sites per region). As predicted, tighter cross-site phase alignment yielded greater trajectory-instability variability, especially on the timescale of the quarter oscillatory period (and also on the timescales of its harmonics if no 1/f noise was included), with little effect on time-averaged trajectory instability. Our simulation further demonstrated that cross-site phase synchronization over a broad range of frequencies elevated trajectory-instability variability across a broad range of timescales with little effect on time-averaged trajectory instability. Specifically, we simulated the waveform at each site as either a Gaussian noise (with a flat spectral-amplitude profile) or a Gaussian random-walk noise (with a 1/f spectral-amplitude profile typically observed in human EEG; e.g., Figure 2), thus simulating phasesynchronized broadband activity across sites. In each case, the trajectory-instability variability for the broadband cross-site phase-synchronized noise was substantially elevated in a largely timescale invariant manner relative to the full phase-scrambled control that randomized cross-site within-frequency phase relations, whereas the time-averaged trajectory instability minimally deviated from the control.
The resting EEG data across all regions mirrored these simulation results (largely timescale invariant elevations in trajectory-instability variability, Figure 4C, with minimal, less than ±2°, effects on time-averaged trajectory instability, Figure 4A), suggesting that intrinsic neural dynamics equivalently facilitate phasesynchronization across a broad range of frequencies. The regional differences in the elevation of trajectoryinstability variability ( Figure 4C) suggest that cross-site broadband synchronization was most facilitated in the posterior region, slightly less facilitated in the anterior and left/right-lateral regions, and least facilitated in the central region. These differences may reflect local computational needs that may require different degrees of broadband phase synchronization.
Interestingly, the variability in the trajectory instability was not elevated relative to the cross-frequency phase-scrambled controls ( Figure 5C). This indicates that any coordination of within-site cross-frequency interactions is such that it does not substantially influence trajectory-instability variability. The potential crossfrequency mechanisms underlying the small but consistent reductions in trajectory-instability variability in specific ranges of timescales in the posterior region (the black curve in Figure 5B) are unclear.

Influences of phase relations on trajectory speed
A theoretical consideration suggests that the coordination of cross-site within-frequency phase synchronization and the coordination of within-site cross-frequency interactions may both influence time- in this example) and increase its temporal variability (from 0 to ∆∆ n × √v in this example). We can also reason that these synchronization effects would be enhanced when EEG changes at individual sites were sharper (i.e., larger ∆∆ n ) and more intermittent (enhancing the impact of synchronization; c.f. if EEG was constant at each site, cross-site synchronization would have no impact).
Given these considerations, increases in cross-site synchronization and the coordination of within-site cross-frequency phase relations that generates sharper and intermittent EEG changes may independently contribute to reductions in time-averaged trajectory speed and increases in its variability relative to the full phase-scrambled controls ( Figure 4B and 4D) and the smaller effects relative to the cross-frequency phasescrambled controls ( Figure 5B and 5D). The fact that the large effects relative to the full phase-scrambled controls were largely timescale invariant implies the coordination of broadband cross-site phase synchronization. In partial support of this possibility, in the simulation experiments discussed above, the phasesynchronized broadband noise with random cross-frequency phase relations generated largely timescaleinvariant reductions in time-averaged trajectory speed and increases in its variability relative to only the full phase-scrambled controls (with little effect of cross-frequency phase scrambling).
The fact that the smaller effects relative to the cross-frequency phase-scrambled controls were more tightly timescale and region invariant ( Figure 5B and 5D) suggests that cross-frequency phase relations are globally coordinated in such a way that EEG waveforms at each site include sharp and intermittent changes across a broad range of timescales. This may be analogous to the way in which different spatial-frequency components in natural visual scenes are phase aligned to generate sharp, high-contrast, and sparse object borders across a broad range of spatial scales (partly due to changes in viewing distance).
If cross-site within-frequency synchronization and within-site cross-frequency phase coordination caused both reductions in time-averaged trajectory speed and increases in its variability, the reductions and increases should be closely associated over time. For example, during a period when time-averaged trajectory speed is strongly reduced its variability should be strongly increased. To support the role of cross-site within-frequency synchronization, we should see a strong temporal association between average-speed reductions and speed-variability increases measured relative to the full phase-scrambled controls beyond the cross-frequency phasescrambled controls. To support the role of within-site cross-frequency coordination, we should see a strong temporal association between average-speed reductions and speed-variability increases measured relative to the cross-frequency phase-scrambled controls. Because the average-speed reductions and speed-variability increases were largely timescale invariant ( Figures 4B, 4D, 5B, and 5D), we averaged the values across timescales (though we verified that the temporal associations were equivalent across timescales). As we computed all statistics (i.e., temporal averages and standard deviations of trajectory instability and speed as well as their phase-scrambled controls) within each consecutive 5-sec time segment (see Materials and Methods), we computed the temporal associations between average-speed reductions and speed-variability increases across time segments per region per participant. We used Spearman r to minimize the effects of potential outliers.
In support of our predictions, we obtained strong negative temporal associations between average speed and speed variability in all five regions from all participants (r < -0.5 for nearly all participants) both relative to full phase-scrambled controls (beyond cross-frequency phase-scrambled controls) ( Figure 8A) and relative to cross-frequency phase-scrambled controls ( Figure 8D). To confirm the specificity of these associations, we obtained no consistent temporal associations when the types of phase-scrambled controls were crossed between the measures of average-speed reductions and speed-variability increases ( Figure 8B-C). Thus, both cross-site within-frequency synchronization and within-site cross-frequency phase coordination (e.g., potentially generating sharp and intermittent changes) globally reduce time-averaged trajectory-speed and increase its variability, but these mechanisms operate independently.
Earlier we reasoned (based on the EEG and simulation results) that cross-site within-frequency synchronization also increases the variability in trajectory instability. If so, increases in trajectory-instability variability should be temporally associated with the aspects of average-speed reductions and speed-variability increases attributable to the coordination of cross-site within-frequency phase synchronization. Indeed, when these variables were measured relative to the full phase-scrambled controls (beyond the cross-frequency phase-scrambled controls), increases in trajectory-instability variability were negatively associated with average-speed reductions (r < 0 for nearly all participants; Figure 9A) and positively associated with speedvariability increases (r > 0 for nearly all participants; Figure 9B) in all five regions.
Thus, the temporal association results confirm that the coordination of cross-site within-frequency synchronization and the coordination of within-site cross-frequency phase relations independently contribute to reducing time-averaged trajectory speed and increasing its variability, with the coordination of cross-site withinfrequency synchronization additionally contributing to increasing the variability in trajectory instability.

Fig 8. Temporal associations (across 5-sec time segments) between reductions in time-averaged trajectory speed and increases in its standard deviation (sd).
Each variable was measured relative to the full phase-scrambled controls (beyond the cross-frequency phase-scrambled controls) to reflect the effects of cross-site within-frequency phase relations, and relative to the cross-frequency phase-scrambled controls to reflect the effects of cross-frequency phase relations. Temporal associations are organized in a 2-by-2 design where the two hypothetical underlying mechanisms, the coordination of cross-site withinfrequency phase relations and the coordination of cross-frequency phase relations, are matched on the diagonal (A and D) and crossed on the off diagonal (B and C). Each panel shows the cross-participant distributions (as interpolated probability density functions) of Spearman r values for the five regions (posterior: black, central: blue, anterior: red, left-lateral: orange, and right-lateral: green). The insets show the corresponding distributions for EEG-state trajectories for the global non-adjacent site configuration (see Figure 10).

Fig 9. Temporal associations (across 5-sec time segments) between increases in the standard deviation (sd) of trajectory instability and reductions in time-averaged
trajectory speed and increases in the sd of trajectory speed. All variables were measured relative to the full phase-scrambled controls (beyond the cross-frequency phase-scrambled controls) to reflect the effects of cross-site within-frequency phase relations. A. Temporal associations between increases in trajectory-instability sd and reductions in time-averaged trajectory speed. B. Temporal associations between increases in trajectory-instability sd and increases in trajectory-speed sd. Each panel shows the cross-participant distributions (as interpolated probability density functions) of Spearman r values for the five regions (posterior: black, central: blue, anterior: red, left-lateral: orange, and right-lateral: green). The insets show the corresponding distributions for EEG-state trajectories for the global non-adjacent site configuration (see Figure 10).

Influences of phase relations on the temporal correlation between trajectory instability and speed
Overall, we obtained positive temporal correlations between trajectory instability and speed relative to both the full phase-scrambled and cross-frequency phase-scrambled controls, suggesting that cross-site within-frequency and cross-frequency phase relations are both adjusted in such a way that faster trajectories are less directionally stable and slower trajectories are more directionally stable. Note that an oscillatory component would generally contribute a negative correlation because the magnitudes of the 1 st and 2 nd derivatives, closely related to our measures of trajectory speed and instability, are negatively correlated in sinusoidal oscillations. Thus, the generally positive temporal correlations (relative to both types of phasescrambled controls) suggest that cross-site within-frequency and cross-frequency phase relations are both adjusted in such a way that EEG-state trajectories are less oscillatory in all regions. We additionally obtained some characteristic regional differences ( Figures 4E and 5E), but further research is necessary to elucidate the mechanisms underlying these effects.

Consideration of volume conduction effects
The results indicative of cross-site phase synchronization in EEG (oscillatory or non-oscillatory) could have been influenced by spurious contributions of volume conduction. The potentially affected results include increases in the temporal variability in trajectory instability as well as reductions in time-averaged trajectory speed and increases in its variability, all of which at least partially implicated cross-site within-frequency synchronization. For a standard 64-site configuration, volume conduction effects after applying the surface-Laplacian transform are minimal beyond the distance between adjacent sites (e.g., Cohen, 2014). We thus reexamined these effects for state trajectories of a global array comprising 16 non-adjacent sites evenly sampled from the whole head (see the map in Figure 10). If spurious volume conduction effects played substantial roles in the original analyses, the effects should minimize in this re-analysis.
All effects replicated. The temporal variability in trajectory instability was consistently elevated relative to the full phase-scrambled control across all timescales, but was largely at the same level as the crossfrequency phase-scrambled control ( Figure 10A), replicating the results obtained for the five regions ( Figures   4C and 5C). The time-averaged trajectory speed was consistently reduced relative to both the full and crossfrequency phase-scrambled controls across all timescales, and consistently more reduced relative to the full than cross-frequency phase scrambled control across all timescales ( Figure 10B), replicating the results obtained for the five regions ( Figures 4B and 5B). The temporal variability in trajectory speed was consistently elevated relative to both full and cross-frequency phase-scrambled controls across all timescales, and consistently more elevated relative to the full than cross-frequency phase scrambled control across all timescales ( Figure 10C), replicating the results obtained for the five regions ( Figures 4D and 5D). All these results were virtually identical for the odd and even numbered participants (the right panels accompanying all main panels in Figure 10). We also replicated all temporal association effects involving reductions in timeaveraged trajectory speed, increases the variability in trajectory speed, and increases in the variability in trajectory instability (insets in Figures 8 and 9). Thus, our conclusions on the effects of phase relations are unlikely to be explained by volume conduction .   Fig 10. An evaluation of potential contributions of volume conduction to the effects attributable to cross-site phase synchronization. Because volume conduction effects are negligible beyond adjacent scalp sites on a standard 64-site array after applying the surface-Laplacian transform (e.g., Cohen, 2014), we re-evaluated the effects that may be susceptible to volume conduction, (A) increases in the temporal variability in trajectory instability, (B) reductions in time-averaged trajectory speed and (C) increases in its temporal variability, for state trajectories of a global array comprising 16 non-adjacent sites (see upper left). The effects should minimize if volume conduction played substantial roles in the original analyses. The effects relative to the full phase-scrambled controls are shown with solid lines and those relative to the cross-frequency phasescrambled controls are shown with dotted lines. The error regions represent ±1 standard error of the mean. A. The temporal variability (measured as temporal standard deviation) in trajectory instability as a function of timescale Dt relative to the phasescrambled controls. It was consistently elevated relative to the full phase-scrambled control across all timescales, but was largely at the level of the cross-frequency phase-scrambled control, replicating the original analyses presented in Figures 4C  and 5C. B. Time-averaged trajectory speed as a function of timescale Dt relative to the phase-scrambled controls. It was consistently reduced relative to both the full and cross-frequency phase-scrambled controls across all timescales, and consistently more reduced relative to the full than cross-frequency phase-scrambled control across all timescales. These results replicate the original anlayses presented in Figures 4B and 5B. C. The temporal variability (measured as temporal standard deviation) in trajectory speed as a function of timescale Dt relative to the phase-scrambled controls. It was consistently elevated relative to both the full and cross-frequency phase-scrambled controls across all timescales, and consistently more elevated relative to the full than cross-frequency phase-scrambled control across all timescales. These results replicate the original analyses presented in Figures 4D and 5D. For each plot, the reliability of the data patterns is demonstrated by the accompanying plots on the right showing virtually identical patterns obtained from the odd and even numbered participants.

Concluding remarks
A major goal of the current study was to elucidate potential macroscopic rules governing global neural dynamics by examining EEG spatiotemporal state trajectories.
A first potential rule is the tight inter-region convergence of time-averaged state-trajectory instability at the level of ~100° on the timescale of ~25 ms ( Figure 3A-D). Both individual differences and temporal variability in trajectory instability minimized on the critical timescale of ~25 ms ( Figure 7B-D) providing indirect evidence suggesting that the convergence is actively maintained. Thus, it is possible that spectral-amplitude profiles (e.g., Figure 2) are globally adjusted to maintain (on average) EEG spatiotemporal dynamics at a nearly directionally unconstrained but slightly conservative level in all regions on a universal timescale of ~25 ms. We speculated that this convergence may establish a global "conduit" for flexible (nearly unconstrained) information exchanges in ~25 ms units.
A second potential rule is that the broadband coordination of cross-site phase synchronization globally (across all regions) generates largely timescale invariant increases in the temporal variability in trajectory instability, reductions in time-averaged trajectory speed, and increases in its temporal variability. The same cross-site phase synchronization mechanisms appear to generate all these effects based on theoretical considerations, simulation results, and the fact that the effects were closely temporally associated with one another ( Figures 8A and 9A-B).
A third potential rule is that the coordination of within-site cross-frequency phase relations, that likely generates sharper and more intermittent changes that are self-similar across timescales, globally (across all regions) generates largely timescale invariant reductions in time-averaged trajectory speed and increases in its temporal variability. The same within-site cross-frequency coordination mechanisms appear to generate both effects as they were closely temporally associated with each other ( Figure 8D).
In comparing the rules two and three, while both the broadband coordination of cross-site phase synchronization and the coordination of within-site cross-frequency phase relations generate largely timescale invariant reductions in time-averaged trajectory speed and increases in its temporal variability, the temporal association results suggest that these contributions are independent ( Figure 8B and 8C).
A fourth potential rule is that both the broadband coordination of cross-site phase synchronization and the coordination of within-site cross-frequency phase relations generate overall positive temporal correlations between trajectory instability and speed, globally (in all regions) nudging EEG state trajectories away from being oscillatory.
A deeper understanding of these potential global rules governing intrinsic neural dynamics will depend on future research investigating how the behaviors of the spectral-amplitude dependent and phase-relation dependent aspects of state-trajectory instability and speed change while people engage in different behavioral tasks. Some of the global rules of intrinsic dynamics may remain task invariant, implying homeostatic processes; disruptions of such rules may result in behavioral dysfunctions. Other rules may characteristically change according to cognitive requirements and/or sensory interactions, elucidating their potential functional relevance. Some preliminary results are suggestive of these distinctions. For example, the inter-region convergence of time-averaged trajectory instability at 100° at ~25 ms (obtained while participants rested with their eyes closed in the current study) appears to be maintained during viewing of a nature video as well as when people view and rate flickering displays for aesthetic preference. In contrast, the effects of phase relations on state-trajectory instability and speed appear to change in a task and stimulus specific manner (Menceloglu, Grabowecky, & Suzuki, unpublished results).
Most EEG research has focused on elucidating the roles of distinct neural sources or networks of sources in behavioral functions. Multi-electrode EEG signals are typically decomposed into estimated neural sources and then time-frequency decomposed so that functional networks can be identified based on spectral amplitude and phase associations. Even when spatiotemporal EEG trajectories are considered (as in the current study) they are typically dimensionally reduced based on linear associations (e.g., using a principalcomponent analysis) to find a relatively small number of functionally relevant components. Behaviors of the identified functional networks and components are examined in relation to performing various behavioral tasks.
These approaches based on identifying behaviorally relevant networks have generated much knowledge in terms of neural correlates of behavior. The prevalence of this approach is evident in the literature (see Introduction) and contents of advanced EEG analysis textbooks (e.g., Cohen, 2014). Here, we took an alternative approach of examining basic features (approximately the 1 st and 2 nd derivatives) of EEG spatiotemporal trajectories and identified simple rules that appear to govern intrinsic neural dynamics while people rested with their eyes closed. An effort to understand the mechanisms that maintain these simple macroscopic rules may help to advance our understanding of the mechanisms that maintain the integrity of global neural dynamics.