Altered Temporal Structure of Neural Phase Synchrony in Patients With Autism Spectrum Disorder

Functional connectivity, quantified by phase synchrony, between brain regions is known to be aberrant in patients with autism spectrum disorder (ASD). Here, we evaluated the long-range temporal correlations of time-varying phase synchrony (TV-PS) of electrocortical oscillations in patients with ASD as well as typically developing people using detrended fluctuation analysis (DFA) after validating the scale-invariance of the TV-PS time series. By comparing the DFA exponents between the two groups, we found that those of the TV-PS time series of high-gamma oscillations were significantly attenuated in patients with ASD. Furthermore, the regions involved in aberrant TV-PS time series were mainly within the social ability and cognition-related cortical networks. These results support the notion that abnormal social functions observed in patients with ASD may be caused by the highly volatile phase synchrony states of electrocortical oscillations.


INTRODUCTION
Autism spectrum disorder (ASD) is a group of neurodevelopmental disorders with sociocommunicative impairments, restricted and repetitive behaviors and interests as core diagnostic features (1). Using modern neuroimaging techniques, scientists found that the symptomatology of this disease was not caused by specific cortical regions, but was associated with the hyper-or hypo-connectivity between cortical regions (2,3). In most published ASD studies, researchers investigated the undirected or directed functional connectivity (FC) between brain regions using measures such as correlation, phase synchrony, and Granger causality, and assumed that these measures were stationary over time (4). However, recent studies have confirmed that FC metrics exhibit variation over time, even during the task-free resting state (5). Moreover, temporal variation in FC metrics can be modulated by various factors, such as brain disorders (6). Surprisingly, the majority of previous studies concentrated on the altered magnitude of functional connectivity between brain regions; few studies have investigated the temporal structure of the time-varying connectivity in the autistic brain.
Long-range temporal correlation (LRTC) is a well-established property of temporal structure in many levels of nervous signals. This property suggests that cortical oscillations have scale-free structures over multiple time scales and work near the critical state (7). This phenomenon offers certain functional benefits for the brain system, such as balanced propagation of external/internal perturbations and maximization of information storage and transfer (8). Detrended fluctuation analysis (DFA) is commonly conducted to quantify LRTCs by amplitude modulation of cortical oscillations (9). For human brain activities, the calculated exponent of DFA is usually larger than 0.5 but <1.0, implying that the neural processes under investigation have LRTCs (10). Empirical studies in both humans and animals have revealed that the scale-free dynamics of brain activity (i.e., its LRTCs or DFA exponent) could be modulated by neurodegenerative diseases, such as ASD, major depressive disorder, and Alzheimer's disease (11)(12)(13). It is assumed that LRTCs are crucial to the efficiency of sustained cognitive operations, including language and social communication, which are impaired in patients with ASD. This theory has been proven by previous studies, which showed that the LRTCs during amplitude modulation of cortical activities were attenuated in patients with ASD (13)(14)(15). However, to the best of our knowledge, no study has tested whether the LRTCs of time-varying phase synchrony (TV-PS) of intrinsic brain oscillations could be modulated by ASD.
Here, we investigated whether the LRTCs of TV-PS between electrocortical signals were altered in patients with ASD using DFA on the resting electroencephalogram (EEG) datasets from the healthy brain network (HBN) (16). We hypothesized that the DFA exponents of the TV-PS time series were altered in these patients, especially the phase synchrony between brain regions known to be associated with autism.

Participants
In the HBN database, 25 ASD patients and 27 typical developing (TD) individuals were selected as participants; note that some of the EEG datasets were used in our previous studies (15,17). All participants were male, aged between 5 and 18 years (ASD: mean age = 11.5 years, SD = 4 years; TD: mean age = 9.2 years, SD = 2.1 years), right-handed, with an IQ higher than 66. Significant differences between two groups on age and IQ scores were not found (ps > 0.05). The diagnosis of ASD was based on the Schedule for Affective Disorders and Schizophrenia-Children's version (18) and the Autism Diagnostic Observation Schedule (ADOS) (1). Note that, the ADOS scores of participants were not available in the current HBN database, thus their ADOS scores were not reported here.

EEG Recording
During the 5 min resting-state EEG collection, the participants alternately kept their eyes open and closed. Previous study have shown that EEG signal length of ∼5 min was enough to produce stable measurements of the DFA exponent (19). A 128-channel Hydro-Cel Geodesic system (EGI Inc., Eugene, Oregon, USA), with an electrode over the vertex of the head as an online reference channel, was used. The ground electrode for the EEG recordings was placed on the forehead. The electrode impedances were kept below 40 kΩ throughout the data recording. The EEG data were recorded using a sampling rate of 500 Hz and a band-pass filtering 0.1-100 Hz. This study was conducted in accordance with the Declaration of Helsinki; ethical approval was obtained from the Chesapeake Institutional Review Board. Written consent was obtained from the participants and their legal guardians (16).

EEG Data Preprocessing
The EEG data were preprocessed using EEGLAB v13.0.0 (20). The preprocessing consisted of the following steps. Firstly, the electrodes over the neck/face and data portions with large drift were excluded. Secondly, the electrodes with low signalto-noise ratio (SNR) and large drift longer than 1 min were identified through visual inspection, and were interpolated using a spherical spline method in EEGLAB (function pop_interp.m). The large drifts <1 min were deleted after interpolating the above "bad electrodes." Thirdly, EEG signals were resampled to 250 Hz, band-pass filtered between 0.5 and 80 Hz using Hamming windowed finite impulse response (FIR) filters (function pop_eegfiltnew.m). The order of the FIR was chosen to be 1,500, i.e., three cycles of the lower edge of the band considered. A notch filter was used to eliminate 60 Hz line noise. Fourthly, the Infomax independent component analysis (ICA) algorithm (function pop_runica.m) was applied to correct physiological artifacts (i.e., eye movements & blinks, electromyography, electrocardiography) and non-physiological artifacts. Lastly, the EEG data were re-referenced to a common average reference.
In the previous data filtering, the alpha band limits were determined to be 8-13 Hz for all participants in both groups. This was based on the fact that the peak alpha frequencies of the occipital electrodes were between 9 and 11 Hz for all participants.
For the gamma band, we studied the low-gamma band (30-55 Hz) and high-gamma band (65-80 Hz). It's well-known that the power of gamma band, especially the high-gamma band, was much lower than other frequency bands and can be easily contaminated by some kinds of artifacts (such as the electromyography), which may lead to inaccurate phase estimation. However, many previous studies have shown that with proper EEG preprocessing and data analysis, it's possible to accurately estimate the power and phase of high-gamma band (22,23).

DFA on the Region-to-Region, Time-Varying, Phase Synchrony Time Series
In most of the previous literatures, the term, "phase synchrony, ", refers to a fixed phase difference between two signals for a certain duration (24). Full synchronization (i.e., the phase difference between neural oscillators is completely consistent over recording time) or full desynchronisation (i.e., the phase difference between neural oscillators varies randomly over recording time) may indicate a pathological state in humans, whereas the non-fixed yet non-random phase relationships between neurophysiological signals could easily emerge within the normal cortical networks (25). Thus, in the current study, the term, "phase synchrony, " is used to describe any pattern of phase relationship between two neuronal oscillations, which could either be fixed or non-fixed values.
Here, the LRTCs of time-varying fluctuations of region-toregion phase synchrony were estimated using the DFA. The procedures for each frequency band and each participant are as follows: 1. The phase time series of neuronal oscillations of each cortical ROI was extracted using the Hilbert transform. Assuming X H (t) is the Hilbert transform of the original band-pass filtered source activity X(t), the time-varying phase φ(t) can be calculated as: . The calculated phase time series are limited to a range [-π, π]; discontinuity (e.g., the phase difference between two consecutive time points is larger than or equal to π) can be found in the phase time series. To convert the phase time series into a continuous time series, the phases were unwrapped by adding multiples of ± 2π when discontinuity occurred. 2. Determine the phase differences between the neuronal oscillations in the two cortical ROIs. Note that, since 84 ROIs were defined in the present study, a total of 3,486, that is, C 2 84 , phase difference time series were calculated for each frequency band and participant. 3. Compute the rate of change of the phase-difference time series. Because the unwrapped phase time series will continue to evolve as time increases, their phase difference time series are unbounded. However, a DFA is typically applicable to bounded processes. Consistent with previous studies, phase synchrony was quantified as the rate of change of the phase difference time series, that is, its first-time derivative (25). 4. The signal profile of the rate of change of the phase difference time series (i.e., its cumulative sum) was computed for each pair of neuronal oscillations.
5. Each signal profile was divided into dozens of windows with length τ and 50% overlap. The set of window length τ was between 1 and 15 s (for alpha, beta, low-gamma, and highgamma bands) or between 2 and 15 s (for delta and theta bands) equidistantly on a logarithmic scale. The linear trend was detrended by subtracting a trendline determined by a least-squares fit. The standard deviation of each detrended window was then calculated. Finally, the fluctuation function for window length τ was computed, that is, the mean value of standard deviations across all detrended segmented signal profiles with length τ .
The scatter plot between logarithmic-transformed window lengths and logarithmic-transformed fluctuation functions was defined as a DFA fluctuation plot. The slope of the least-squares line for this scatter plot, which is referred to as the DFA exponent in the literature, provides an estimation of the LRTCs or scalefree property of a given region-to-region TV-PS time series. The linear scaling in log space was confirmed by the high value of the coefficient of determination (i.e., R 2 > 0.9 for most of the TV-PS time series) and stringent validation tests, as shown below. The above DFA exponent calculation was conducted using the Neurophysiological Biomarker Toolbox (https://github.com/ NBT-Analytics/NBTpublic).

Validating the Presence of Scale-Invariance
The

Statistical Tests
To examine whether the LRTCs of the TV-PS time series were significantly altered in patients with ASD, the networkbased statistic (NBS) was performed for the DFA exponent at each frequency band. Using an approach similar to the cluster-based permutation test used in the activation tests, the NBS is a non-parametric, statistical method to avoid the multiple comparison problem encountered when conducting mass univariate statistical testing in functional networks (27). This method consists of the following steps (27): First, Fstatistics for all the 3,486 (C 2 84 ) edges on the graph based on the differences in DFA exponent values between groups, with the age of participants included as a covariate, were computed. In NBS, the statistical model is specified in terms of the general linear model (GLM). Here, the covariate was included as nuisance regressor in the design matrix of GLM. Second, a primary threshold (p < 0.005, uncorrected) was used to identify all edges displaying potential differences in the DFA exponent value.
From these identified edges, the so-called, "connected graph components, " which were defined as a set of supra-threshold edges for which a path could be found between any two cortical ROIs, were identified. Third, the size of each connected graph component was defined as the total number of edges it comprises (i.e., "component extent"). In NBS, the component size could also be measured by the sum of test statistic values across all connections comprising the component. This is referred to as the component intensity. In the current study, these two approaches revealed exactly the same results. Fourth, the null distribution of the size of the connected graph component was empirically derived using a permutation approach with 5,000 permutations. For each permutation, all participants were randomly allocated into two groups; the above three steps were then conducted. The component with a maximum cluster size was noted for each permutation, which yielded an empirical null distribution for the largest component size. Lastly, the one-sided, family wise error rate (FWER)-corrected p-value for an originally identified component was estimated as the proportion of permutations for which the largest component was of equal size or greater. The final results controlled the FWER (weak sense) at the cluster level (p < 0.05).

The Presence of Scale-Invariance
The results of the ML-DFA conducted on the region-to-region TV-PS time series are shown in Figure 1. We found that for most ROI pairs, the percentage of participants for which the presence of scale-invariance was accepted was higher than 90%. Moreover, for most ROI pairs in the beta, low-gamma, and high-gamma bands, this percentage reached 100% in the TD group. Additionally, this percentage was much higher in higher frequency bands (i.e., alpha, beta, low-gamma, and high-gamma bands) than in the lower frequency bands (i.e., delta and theta bands), and was much higher in the TD group than in the ASD group.

The Group Effect on DFA Exponents
For the high-gamma band, the statistical tests on the DFA exponents showed that those in 28 ROI pairs in the TD group were significantly larger than those in the ASD group (Figures 2-4). No significant results were found for the other frequency bands. The NBS also did not detect any significant results when the alternative hypothesis was ASD > TD.
Because the magnitude (i.e., power spectra) of the neuronal oscillations of each frequency band could influence the successful estimation of the phase time series and of the DFA exponents, we tested whether significant group differences in the power spectra of each frequency band could be detected. The results showed that significant group differences could not be detected in all six bands and ROIs.

DISCUSSION
Here, we investigated whether the LRTCs of the region-to-region TV-PS time series were significantly altered in patients with ASD. First, the results of the ML-DFA confirmed the presence of scale invariance in the TV-PS time series. Second, the LRTCs of TV-PS within the high-gamma band were found to be reduced in these patients.

The Scale-Free Dynamics and LRTCs in Brain Oscillations
Exploring the datasets in behavioral science and neuroscience, previous researchers found that scale invariance/LRTCs are omnipresent in nature, which supports the theory that cortical activities may sit in a dynamic state close to the criticality (28). The advantages of criticality demonstrated through computational simulations and experimental studies include the following aspects: (1) the maximum dynamic range of information processing/communication and its efficiency within systems, and (2) a readiness to respond to environmental changes (10,25). The presence of scale-invariance, which was validated by ML-DFA, was observed in most of the TV-PS time series across all the region-to-region pairs, the six frequency bands investigated, and the participants in both groups. Moreover, in Figure 1, we found that the percentage of participants in the TD group for which the presence of scale-invariance was accepted was larger than that of the ASD group for all six frequency bands and nearly all the ROI pairs. Thus, the results detected may suggest that although scale-invariant dynamics were inherent in cortical-cortical TV-PS time series of autistic brains, this property may be vulnerable in these patients.
Scale-free dynamics and LRTCs in neuronal oscillations are associated with the efficiency of certain operations (e.g., learning, memory processes, and information transfer) within the brain system (25). Several studies have demonstrated the presence of LRTCs and scale invariance for the TV-PS between pairs of neural oscillations and found that this property of phase synchrony could be influenced by some experimental operations, such as finger movement (25,29). Previous studies have indicated that the magnitude of LRTCs is positively correlated with the ability of the human brain to maintain transiently stable oscillations in support of active neuronal representations during sustained cognitive operations (e.g., language, social cognition, and communication), which are impaired in patients with ASD (30). Moreover, it has been widely revealed that aberrant functional connectivity is inherent in an autistic brain (31). Thus, we investigated whether the LRTCs of phase synchrony were disrupted in autistic brains using source-level electrocortical oscillations, which could provide rich temporal, spectral, and spatial information about brain activities.

The Attenuated LRTCs of Time-Varying Phase Synchrony in Autistic Brain
The statistical tests on the DFA exponent showed that those in 28 ROI pairs in the ASD group were significantly smaller than those in the TD group for the high-gamma band. Cortical regions involved in the aberrant phase synchrony time series were mainly located in the right hemisphere.
Here, the ROIs that explained 26 out of 28 aberrant connections and with more than four aberrant phase synchrony time series were defined as "crucial hubs." These brain areas were the precuneus, retrosplenial cortex (RSC), primary somatosensory cortex (S1), temporal pole (TP), and ventrolateral prefrontal cortex (VLPFC). The precuneus, located in the posteromedial portion of the parietal lobe, is involved in visuospatial imagery, episodic memory retrieval, and self-processing operations (32). The RSC is closely associated with a wide range of brain functions, such as episodic memory, navigation, imagination, and planning for the future (33). The TP plays a crucial role in both social and emotional processes, such as face recognition and theory of mind (34)(35)(36). Specifically, the right TP is believed to be the storehouse or site of recollection of personal, episodic memories, and is involved in high-level sensory representations with emotional responses and social memory (34). These three regions are the main hubs of the default mode network (DMN) and are suggested to be involved in the pathophysiology of autism (2,13). For S1, previous studies have shown that this region often exhibits atypical responses to touch in patients with ASD (37). Moreover, since S1 is a crucial region in the mirror neuron system (MNS), it may also closely modulate social cognition and interaction in patients with ASD. The VLPFC, especially the right VLPFC, plays an important role in emotional regulation in the social context (38). It is also an important cortical region of the MNS that underlies the dysfunction of social cognition in autism (39). Almost all the aberrant DFA exponents (i.e., 26 of 28) could be explained by the connections between these five hubs and other brain regions. It should be noted that these results may also be caused by the fact that these five brain regions manifested as certain specific patterns (e.g., periodic activity) in the autistic brain, which requires further investigation.
The presence of scale-free dynamics has been found in the moment-to-moment fluctuations of phase synchrony; the healthy resting brain state is accompanied by weak and variable neural synchrony which shows slow fluctuations over time (29). The LRTCs of the TV-PS time series may reflect a state of readiness that facilitates rapid task-dependent shifts toward and away from full synchronization or full desynchronisation (25). Thus, the attenuated LRTCs of phase synchrony between certain regions in patients with ASD may indicate that the autistic brain, especially those regions involved in autism, deviated from an optimum state of readiness to external or internal events.
In a study conducted by our team (15), the DFA exponents of the instantaneous amplitude of the beta and low-gamma oscillations of DMN, MNS, and salience network (SN) were significantly attenuated in patients with ASD, as compared with the TD participants. First, the regions with altered DFA exponents (i.e., LRTCs) were highly overlapped between DFA conducted on the instantaneous amplitude time series and DFA conducted on the phase synchrony time series. Second, compared with DFA conducted on the instantaneous amplitude time series, the EEG band with significant results turned into a high-gamma band when TV-PS was investigated. (40) found that temporal changes in phase synchrony exhibit less temporally structured and more complex correlations, indicating fast and flexible coding (40). This could be the reason why the frequencies showing aberrant LRTCs occur at much faster rates (i.e., high-gamma) when the phase synchrony time series were investigated.
Here, nearly all the aberrant phase synchrony time series were caused by cortical regions over the right hemisphere. Atypical laterality using neuroimaging techniques has been widely studied in the ASD group (41). For the EEG gamma band, previous studies showed that the resting gamma power was reduced in patients with ASD and was specific to the right lateral hemisphere (42). Note that none of the previous studies have investigated the LRTCs of phase synchrony time series in patients with psychiatric disorders. Thus, this is the first study to show that abnormal laterality of LRTCs could be seen in the phase synchrony between cortical regions in patients with ASD.
It should be mentioned that the largest window size used in the current study was 15 s. If the largest window size was increased (e.g., 20 s), the percentage of fluctuation plots being accepted as linear and the value of DFA exponents may be decreased. However, in an exploratory analysis, we found that if the largest window size was set to 20 s, the results of the group comparison were similar to those reported here.

Comparison With Other Studies Focused on the Time-Varying Features in Autistic Brain
In these years, more and more studies investigated various aspects of time-varying neural features in ASD. For example, Sase and Kitajo investigated the metastability in the brain oscillations through the phase-phase coupling (PPC) and phaseamplitude coupling (PAC) of resting-state EEG signals in typically developing individuals with different levels of autisticlike traits (43). They found that the metastable dynamics of synchronization and amplitude modulation is correlated with levels of autistic traits. Moreover, fewer transitions between states occurred in individuals with relatively longer attention span. In the other study, using resting-state fMRI signals, Watanabe and Rees characterized the cortical dynamics in autism through an energy-landscape analysis (44). They found that patients with ASD show fewer neural transitions between two major brain states, and this abnormal transition can predict the severity of autism. Here, we investigated the scale-invariant characters of resting-state EEG signals in autism, and found that the LRTCs of region-to-region, time-varying, phase synchrony time series within the high-gamma band were reduced in the autistic brain. Due to the complexity of time-varying neural features in human brain, more studies are needed in order to obtain a comprehensive understanding of neural dynamics in autism.

The Limitations of the Current Study
Several limitations of the current study should be mentioned. Firstly, the ADOS scores of participants were not available in the current HBN database, thus their ADOS scores were not reported and analyzed here. Without these scores, we could not assess the associations between the severity of symptoms and DFA exponents, which limited the reliability and interpretation of the findings as well as the guidance for clinical treatment and intervention. Secondly, during the validation of the presence of scale-invariance of TV-PS time series, we only showed some descriptive measures (e.g., the percentage of participants for which the presence of scale-invariance) in Figure 1. We did not conduct any statistical tests between groups, which may limit the interpretation of the findings between groups.

CONCLUSION
In the present study, we found that the LRTCs of region-to-region, time-varying, phase synchrony time series within the high-gamma band were reduced in the autistic brain. Furthermore, the cortical regions involved in aberrant LRTCs are mainly located in well-recognized brain networks associated with autism. These results indicate that reduced LRTCs of time-varying phase synchrony within high-gamma oscillations may play an important role in the dysfunction of social, communication, and emotional abilities commonly observed in patients with ASD.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: http://fcon_1000.projects.nitrc.org/indi/cmi_ healthy_brain_network.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Chesapeake Institutional Review Board. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
The work presented here was carried out in collaboration between all authors. HJ and FG initiated and analyzed the data and wrote the manuscript. DY wrote and revised the manuscript. All authors have read and approved the final published manuscript.