{'fraction_non_alphanumeric': 0.03752707225831857, 'fraction_numerical': 0.01744437881472731, 'mean_word_length': 5.197047340165935, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Early childhood developmental functional connectivity of autistic brains with non-negative matrix factorization

Highlights • NMF-based method to identify abnormal connectivity patterns of ASD brains.• ASD children were characterized by interhemispheric under-connectivity.• Long-range connections modulated by slower rhythm; local-range by faster frequency.• Alpha-band interhemispheric pattern play pivotal role in ASD neurodevelopment.


Introduction
Autism spectrum disorder (ASD) is a complex neurodevelopmental disorder that is characterized by repetitive behaviors, restricted interests, and impairments in social interaction and communication (Keown et al., 2017). Evidence that ASD could be a result of altered patterns of neural connectivity is growing (Geschwind and Levitt, 2007;O'Reilly et al., 2017;Minshew and Keller, 2010). The developmental-disconnection model of ASD indicates that changes in functional connectivity in individuals with ASD follow a pattern of shortrange over-connectivity (Belmonte et al., 2004) and long-range underconnectivity (Geschwind and Levitt, 2007) or even global under-connectivity (Hughes, 2007;Just et al., 2004). Recent findings, however, emphasize a mixture of hypo-and hyper-connectivity (Kana et al., 2014). This abnormal connection pattern may then reflect inadequacy in the fine-tuning of networks during development in ASD, which is a process that is typically characterized by increased integration within local networks and increased segregation across separate networks (Bos et al., 2014). However, several findings remain unreplicated, and conclusions regarding the nature of altered patterns of connectivity in ASD have been divergent.
A substantial body of research has partially addressed questions regarding brain connectivity in ASD, including approaches that have investigated structural connectivity by diffusion imaging (Anagnostou and Taylor, 2011;Müller, 2014) and correlated activity by functional magnetic resonance imaging (fMRI) (Thai et al., 2009;Uddin et al., 2013). Specifically, fMRI studies have identified atypicalities in functional-connectivity patterns in individuals with ASD (Hahamy et al., 2015), which include reductions in various functional connections (Belmonte et al., 2004;Just et al., 2004), particularly in interhemispheric connections (Hahamy et al., 2015;Dinstein et al., 2011;Di Martino et al., 2014). Moreover, fMRI studies have reflected a shift in the understanding of the biological basis and cognitive deficits of ASD toward the patterns of functional connectivity between areas of https://doi.org/10.1016/j.nicl.2020.102251 Received 18 September 2019; Received in revised form 17 March 2020; Accepted 18 March 2020 the so-called default mode network (DMN), which are neural networks that are activated at rest (Raichle, 2009).
Brain connectivity can also be evaluated by multiple electroencephalography (EEG) signals. According to its high temporal resolution and direct relationship with neuronal potential activity, EEGconnectivity analyses can provide abundant information about the dynamic activation and deactivation of functional networks (O'Reilly et al., 2017). Recently, aberrant gamma activity has been reported in autistic children, which has been interpreted as the supporting hypothesis of aberrant brain connectivity (Brown et al., 2005). EEG coherence has been directly used as a measure of connectivity in several EEG studies, and the reports support the evidence of both underand over-connectivity in different frequency bands in ASD populations (Murias et al., 2007). However, a more complete understanding of the patterns of connectivity in ASD remains lacking. Thus, our goal in this study is to identify connectivity patterns and effective brain connectivity in ASD, with a focus on factors impacting over-versus underconnectivity during the early stage of neurodevelopment (ages: 3-6 years old).
Given the complexities and conflict in the ASD literature, datadriven techniques provide exploratory approaches that are appropriate for identifying connectivity patterns, even in the absence of strong directional assumptions (Keown et al., 2017). To address these conflicting reports, a data-driven approach that can extract the parts of functional brain networks that encode resting states associated with altered patterns of connectivity-and track their expression over time-is required. Such a capability would improve our understanding of those patterns of brain networks that are significant for different groups of ASDs and how these patterns change over time.
In this work, we identified four frequency-band-specific (δ, θ, α, and β bands) patterns of functional brain networks associated with autistic resting states using an unsupervised machine-learning technique known as non-negative matrix factorization (NMF) (Lee and Seung, 1999). In recent years, this new tool has been applied to identify time-varying functional connectivity patterns in the brain (Marimpis et al., 2016). Intuitively, NMF decomposes functional brain networks into the following: (1) a set of subnetworks (patterns) overlapping in space and time and (2) corresponding coefficient time series that quantify the contribution of each subnetwork (pattern) at each time point (Chai et al., 2017;Khambhati et al., 2017Khambhati et al., , 2018a. One possible definition for an EEG subnetwork-and the one that we use in this study-is that it reflects the pattern of connectivity. This computational tool allows us to track how groups of networks interacting across brain areas are dynamically expressed during rest states. As compared to hard-partitioning schemes, the advantage of this method is that it provides information about brain-network dynamics in a continuous, overlapping manner in space and time rather than in discrete partitions. Furthermore, owing to the parts-based nature of the technique, we obtained subnetworks that resembled the localized features of largescale brain networks rather than generalized patterns of the overall network (Lee and Seung, 1999;Chai et al., 2017). Based on prior work suggesting distance-dependent organization of brain networks into local, function-specific interactions (characteristic of clusters and modules (Rubinov and Sporns, 2010;Telesford et al., 2011)), we hypothesize that subgraphs are selectively sensitive to functional interactions over different distances. Specifically, we applied this technique to a control group of early childhood neurodevelopmental autistic and age-matched children to address the following hypotheses: the first hypothesis considered autism to be characterized by long-range underconnectivity (particularly in interhemispheric connections) potentially combined with short-range over-connectivity, reflecting atypical development of functional networks. The second hypothesis is that we would find differences between groups (i.e., between autistic and control children) in these canonical patterns of functional connectivity. The third hypothesis is that there would be an age-related trend in the measured patterns of connectivity, potentially owing to slight differences in early childhood developmental trajectories between autistic and control children.

Participants
In this study, we collected EEG data from 256 children, including 115 children with ASD (age range 3-6 years old, mean age: 4.7 years, SD: 1.53 years) and 141 age-matched typically developing (TD) controls (age range: 3-6 years old, mean age: 4.85 years, SD: 1.47 years). All children in the ASD group were recruited because they had a prior diagnosis of ASD, according to the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition (DSM-V) (Association, 2013). The research protocol was approved by the ethics committee of Beijing Normal University. Informed consent was obtained from all parents before the start of the experiment. All protocols of this study conform to the Declaration of Helsinki guidelines.

EEG collection and preprocessing
The EEG was continuously recorded from a 128-channel EEG system (Electrical Geodesics Inc., Eugene, OR). Scalp impedance was checked online via the Use Net Station (EGI, Inc.) and was maintained below 50kΩ. EEG data were referenced online to Cz and digitized at a sampling rate of 1,000 Hz. Open-eyes resting-state EEGs were recorded for at least 5 min. Children were instructed to sit comfortably and relax during the recording and were usually accompanied by their caregiver in a quiet room.
The resting-state EEG data were preprocessed offline using EEGLAB, v14.0.0b (Delorme and Makeig, 2004) and MATLAB 2017b (Mathworks, Inc., USA). According to the 10-10 standard international system, the same 62 electrodes were selected from the 128 channels of the Geodesic Sensor Net (GSN) for all participants to ensure the broadest spatial coverage of the whole brain (frontal, central, temporal, and occipital). The data were re-referenced to a common average reference and down-sampled to 512 Hz. A notch filter centered at 50 Hz was applied to minimize line-noise artifacts, and the data were then band-pass filtered (0.5-45 Hz). Independent component analysis (ICA) was performed on cleaned data using FastICA (Rogasch et al., 2014) to decompose 62 spatiotemporal features of the EEG-independent components (EEG ICs), and then the corresponding blink/oculomotor, muscle, or transient electrode artifacts were subtracted from the data via visual inspection in terms of channel-based scalp topography measures and power spectral density (PSD) measures (median, variance, and kurtosis). To ensure the quality of the resting-state data, the middle portion of the EEG data (of 2 min duration) was used for subsequent analysis, as it was less noisy than the entire time series. The EEG data were then cut into 30 non-overlapping segments of 4 s each.

Time-varying network construction
After preprocessing the rest-state EEG data for each participant, we obtained a 120-s time series from each electrode and divided it into 30 non-overlapping 4-s windows. Thereafter, we examined the timevarying brain functional connectivity, for which the connectivity between all pairs of brain regions (i.e., channels or electrodes) was required to be estimated as a function of time and frequency using the sliding-window technique (Cai et al., 2018). For each time window, we measured functional connectivity between each pair of brain regions based on the imaginary part of wavelet coherency within four frequency bands: δ (2-4 Hz), θ (4-8 Hz), α (8-13 Hz), and β (13-30 Hz).
In this study, we only analyzed connectivity at the sensor level owing to the lack of MRI data. Certain known limitations exist to scalplevel analysis of functional networks (Colclough et al., 2016;Brookes et al., 2014): (i) the locations of EEG channels are not related T. Zhou, et al. NeuroImage: Clinical 26 (2020) 102251 directly to the location of the underlying sources and (ii) spurious estimates of functional connectivity can possibly occur between channels owing to the effects of field spread and volume conduction. To overcome the problems of sensor-level connectivity, Stam et al. (2007) and Nolte et al. (2004) came up with two methods-the phase lag index (PLI) and the imaginary part of coherency, respectively. These methods discard zero-lag connectivity and are, therefore, insensitive to volume conduction. In this study, we adopted the imaginary part of coherency (Nolte et al., 2004;Helfrich et al., 2016) based on wavelet transforms. The main steps are detailed below. Let x t and y t denote two stationary multivariate time series. The continuous wavelet transforms, denoted as X ω and Y ω , are defined as where b is a time index and a denotes the wavelet scale. The function ψ 0 is chosen to be the complex Morlet wavelet: where the parameter ω 0 governs the relative time and frequency resolution.
We assume throughout that X ω and Y ω each have zero mean. Let denote complex-valued covariance matrices, where the superscript'*'denotes vector/matrix transposition and complex conjugation. These matrices correspond to the wavelet cross-spectral density matrices.
Subsequently, we use the imaginary part of the covariance matrices to calculate general lagged coherence, proposed by Nolte et al. (2004), with the zero-lag effect removed: where s yyw and s xxw denote pure real variances and s yxw denotes the complex-valued covariance, with the real and imaginary parts denoted as Re(s yxw ) and Im(s yxw ), respectively. For each participant, coherence values were stored in an N × N × T connectivity tensor A, where N = 62 channels and T = 30 time windows. Thereafter, these tensors were converted into matrices (vectorizing the upper diagonal of each N × N connectivity matrix of a time window). For real data with multiple participants, we built the matrix V by concatenating the columns from all participants (from all conditions) together to form a single time-varying network configuration matrix V.
where S = 115 + 141 = 256 participants. These pooled neural data, which facilitated network sharing by all participants could be identified ( Fig. 1A−E illustrates these matrices). In our weighted-network analysis, we retained and analyzed all possible connection weights between nodes without applying any threshold. This matrix was then used as input to the NMF algorithm and for subsequent processing (Fig. 1E, F).

Decomposing dynamic network into network patterns
To identify network patterns, we applied an unsupervised machine-learning algorithm, NMF (Lee and Seung, 1999), to the constructed time-varying network (Fig. 1D). This technique allowed us to pursue a near-decomposition of the whole network into topological subnetworks and the corresponding dynamically temporal coefficients (Fig. 1E, F). Each pattern is an additive component of the original network, weighted by its associated time-varying expression coefficient, and represents a pattern of functional interactions between brain regions. The NMF-based connectivity pattern learning paradigm is a basis decomposition of a collection of dynamic networks that separates co-varying network edges into subnetworks, or basis functions, with associated temporal coefficients, or basis weights. Unlike other graph clustering approaches that seek a hard partition of nodes and edges into clusters (Mucha et al., 2010;Bassett et al., 2013), the temporal variability coefficients provide a soft partition of the network edges, such that the original functional network of any time window can be reconstructed through a linear combination of all of the connectivity patterns weighted by their associated temporal variability coefficients in that time window (Leonardi et al., 2013;Leonardi et al., 2014). This computation implies that at a specific time window, connectivity patterns with a high temporal variability coefficient contribute their pattern of functional interactions more than subgraphs with a low temporal variability coefficient (Khambhati et al., 2018a).
Mathematically, NMF is a multivariate-decomposition technique that approximates a data matrix V (here, of size E × M) as the product of two non-negative matrices W and H. Next, we solved the matrixfactorization problem V ≈ WH, W ≥ 0, H ≥ 0 by optimizing the following cost function: 1] is the number of patterns to decompose, β is a penalty weight that enforces sparsity on the temporal variability coefficients, and α is a regularization parameter that provides an upper bound on the connection strengths within the patterns (Kim and Park, 2011). To solve the NMF equation, we performed 100 iterations using alternating non-negative least squares with the block-pivoting method to quickly and efficiently decompose the original matrices (Kim et al., 2014). We initialized W and H with non-negative weights drawn from a uniform-random distribution on the interval [0, 1].
To select the parameters k, β, and α, we used a random-sampling scheme, demonstrating effective optimization of high-dimensional parameter spaces (Khambhati et al., 2018a;Bergstra and Bengio, 2012), where we re-ran the NMF algorithm for 1000 parameter sets, where k was randomly selected from 2 to 20, β was randomly selected from 0.01 to 1, and α was randomly selected from 0.01 to 1 (Fig. 2). We evaluated subnetwork-learning performance based on a fivefold cross-validation scheme. Iteratively, four folds were used to extract subnetworks and the one fold left out was used to calculate the cross-validation error ( − V WH F 2 ). The optimal parameter set yielded subnetworks that minimized the cross-validation error and reliably spanned the observed network topology space (Khambhati et al., 2018a). Based on these criteria, we chose an optimal parameter set (k, β, ᾱ ) that yielded a low residual error according to our scheme ( Fig. 1G−I).
Owing to the non-deterministic nature of this approach, we integrated patterns estimates over multiple runs of the algorithm using consensus clustering-a general method of testing robustness and stability of clusters over several runs of one or more non-deterministic clustering algorithms (Monti et al., 2003). Our adapted consensus clustering procedure involved the following steps: (i) run the NMF approach R times on each network configuration matrix, (ii) concatenate the subnetwork matrix W across R runs into an aggregate matrix with dimensions × E R k ( *¯), (iii) apply NMF to the aggregate matrix to determine a final set of subnetworks W consensus and the dynamic coefficient H consensus . In this study, we set R = 100.
As a sanity check, we compared our NMF-based results of decomposing the α-band network with the classic and sparse NMF schemes (Lee and Seung, 1999;Berry et al., 2007;Hoyer, 2004), which in theory should provide very similar results. The results of this comparison are displayed in Figure S2. Overall, the three NMF derivative approaches yield very similar patterns. We evaluated their performance based on a fivefold cross-validation scheme. The classification module classified the feature vector, i.e., extracted six network components energy, using the support vector machine (SVM) classifier to identify the children with ASD (ASD vs TD). The classification accuracies of our NMF-based, original NMF, and sparse NMF approaches are 72.5%, 70.1%, and 70.25%, respectively. The results reflect that the proposed NMF-based approach performs better than the original and sparse NMF approaches.

Test-retest reliability of connectivity patterns
It is essential to consider the reproducibility of patterns measured using NMF. To quantify the reproducibility of the connectivity patterns, we measured the extent to which the pattern of the subgraph edge weights measured in one dataset predicts the pattern of the subgraph edge weights measured in a second dataset. Specifically, we first applied NMF using the optimal parameter set to two different datasets (V 1 and V 2 ), resulting in two subnetwork matrices (W 1 and W 2 ). It is noted that the connectivity patterns along the columns of W 1 may not necessarily be ordered similarly to the subgraphs along the columns of W 2 because of the stochastic nature of the NMF algorithms. To reorder the connectivity patterns from the second dataset such that they correspond to the same order as the patterns from the first dataset, we sought a Next, we unfolded matrices for each time window into column vectors and concatenated time windows from all participants. We then applied nonnegative matrix factorization (NMF) algorithm, which decomposes concatenated matrix into matrix W of subnetworks (E) and matrix H of time-dependent coefficients (F) that quantify activation of each subnetwork in each time window. (G-I) NMF-based subnetwork detection requires optimizing three parameters-number of subnetworks, k; temporal sparsity of subnetwork expression, β; and regularization of subnetwork-edge weights, α. Kernel-density estimate of each bivariate distribution is represented by contour plot, where darker shades of blue indicate higher probability mass of random-sampling distribution. The best parameter is average parameter value that produces a cross-validation error at bottom of sample distribution and is indicated by green dashed line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) T. Zhou, et al. NeuroImage: Clinical 26 (2020) 102251 Fig. 2. Age-related differences in patterns of homotopic−interhemispheric connections. Among the 23 patterns of networks, four patterns were highlighted as examples of age-related variation in network connectivity. For visualization, the thresholds of the subnetwork maps were set at the 90th percentile, and the connections are depicted in a circular format. Significant differences for each age group in terms of the energy (middle panel) and temporal variability (denoted by entropy; right panel) of the four patterns as well as statistical analysis of the two groups without considering age factors (inset) were observed. Only significant group effects are marked. * indicates the network pairs with significant connectivity differences between groups (*, **, and *** indicate p < 0.05, p < 0.01, and p < 0.001, respectively). mapping X ij of the subnetwork W i 1 to the subnetwork W j 2 , where X is a Boolean matrix that prescribes whether the ith subnetwork from the first dataset is uniquely assigned to the j th subnetwork from the second dataset. The cost C i, j associated with assigning W i 1 to W j 2 is equal to − W W i j 1 2 . To determine a unique X, we minimized the cost function ∑ ∑ X i j i j , using the Hungarian algorithm (Kuhn, 1955). After calculating an optimal assignment between the connectivity patterns of the two datasets, we measured the similarity in the pattern of the edge weights between the assigned subnetwork pairs (i, j) by calculating the Pearson correlation coefficient (Khambhati et al., 2018a, b). This approach would enable us to assess the reproducibility of each individual subnetwork based on the magnitude of the Pearson correlation similarity measure relative to that expected by chance (Khambhati et al., 2018a).

Statistics
To track the time-varying changes in the expression or activation of the subnetworks, we first calculated the signal energy of the time-dependent coefficients for each subnetwork and child. Signal energy is defined as ∑ = H n L n 1 2 , where H n is the activation coefficient in time window and L is the length of the signal for a participant (Chai et al., 2017). Second, to quantify the dynamic switching behavior of the network component expression, we calculate the signal entropy of the activation coefficient time series to characterize temporal variability. The signal entropy is defined as ∑ where P(x i ) is a probability mass function on the activation coefficient computed using the histogram-based entropy estimator (Nelson et al., 2010).
To facilitate between-subject comparisons, we standardized the energy and entropy values for each subnetwork-we computed the means of the time-dependent coefficients for each participant and each network component and then divided all time-dependent coefficients for each subject and each subnetwork by their respective means. A twoway analysis of variance (ANOVA) was performed to evaluate the energy and entropy values in different bands, including the factor groups (ASD, TD) and age (3-6 years), using Bonferroni's post-hoc test for multiple comparisons between groups. The α-level of the significance was set at 0.05. Furthermore, the Pearson correlation was calculated between the energy (or entropy) values and age (3-6 years) to explore the trajectory of brain development in early childhood. Only significant group effects are reported. All statistical analyses and graphical representations were produced using GraphPad Prism 7.0 Software. Data are expressed as mean ± standard error unless otherwise specified.

Extract patterns of connectivity from the resting-state EEG via NMF
To identify the patterns of resting-state connectivity network, we first measured the resting-state EEG from a group of 115 children with ASD and 141 age-matched typically developing controls (TD; Fig. 1A). In particular, we selected 62 sensors of interest to ensure the maximum spatial coverage of the frontal, central, temporal, and occipital regions (Fig. 1B). We then constructed time-varying brain networks for each participant, in which network nodes represented brain regions and network edges between the nodes represented the wavelet-based coherence coefficient among the regional EEG time series (Fig. 1C). Specifically, we computed a 62 × 62 adjacency matrix for each participant within each EEG time window (corresponding to a 120-s time series that was divided into 30 non-overlapping 4-s time windows) for four specific-frequency bands (δ, θ, α, and β). This process resulted in 30 window-level adjacency matrices per participant. Finally, we aggregated all functional brain networks into a network-configuration matrix (Fig. 1D) with a size of 1891 × 7680. The first dimension of size 1891 corresponded to all unique, pairwise edges between the 62 channels; the second dimension of size 7680 corresponded to all combinations of the 30 repeated time windows and 256 participants.
To uncover the patterns of connectivity and their dynamic expressions from group-level brain networks, we applied the NMF technique to the network-constructed configuration matrix. This technique enabled us to pursue a parts-based decomposition of network edges into additive functional subnetworks (i.e., patterns; Fig. 1E) with accompanying expression coefficients over time and over subjects (Fig. 1F) (Chai et al., 2017;Khambhati et al., 2018a). Each subnetwork was composed of a 62 × 62 adjacency matrix, and the expression coefficients of each subnetwork were composed of a vector of length 7680. Thus, the subnetworks comprised the topological components of the functional brain network, and the temporal coefficients quantified their connection strengths and dynamic expressions over time.
A critical step in using NMF is the optimization of model parameters that identify a robust set of connectivity patterns. We apply a random sampling scheme to characterize the rich parameter space of the number of subnetworks k, temporal sparsity of subnetwork expression β, and regularization of subnetwork-edge weights α across the timevarying networks of all participants to ensure the generalizability of patterns without overfitting the model to the constructed configuration matrix. In this paper, we only present the optimization results of the αband connectivity (Fig. 1G−I).
By designing a fivefold, leave-one-fold-out cross-validation scheme, we measured the relationship between the fivefold cross-validation error and each parameter k, α, and β ( Fig. 1G−I). We observed a weak relationship between the cross-validation error and temporal sparsity (Pearson r = −0.15, p = 0.15). We also observe a significant increase in cross-validation error as the subnetwork regularization parameter (Pearson r = 0.36, p < 0.001). Similarly, we observe a significant positive relationship between the cross-validation error and number of subnetworks (Pearson r = 0.66, p < 1× − 10 16 ). Collectively, these results suggest a potential strategy for choosing parameters that achieve a balance between the spatial and temporal generalizability and specificity of the subgraphs. Therefore, we averaged the randomly sampled parameters associated with the lowest 25% cross-validation error and found the optimal number of subnetworks to be six, the temporal sparsity to be 0.65, and the regularization of the subnetwork edge weights to be 0.45 (Fig. 1G−I). Thus, we performed parameter optimization on four frequency-band-specific connectivity−configuration matrices (Table 1). The results presented a near-linear increase in the number of patterns with increasing bandwidth from the δ to β bands. These results can be explained by the fact that a wider bandwidth yields a greater capacity to hold information and, thus, results in a larger number of patterns. To summarize, 23 patterns of connectivity were identified in the autistic group-level EEG data over four frequency-band-specific connectivity patterns (patterns A-D for δ connectivity, patterns E-I for θ connectivity, patterns J-O for α connectivity, patterns P-W for β connectivity). We then visualized the resulting adjacency matrices as circular-ring graphs ( Figure S1). The pairwise correlation matrix between the spatial topology of network components identified with NMF in four specific frequency bands is displayed in Figure S3. Overall, the four frequency-band-specific network components are similar. The α-band-specific network components were topologically more similar than the other frequency-based counterparts. The autocorrelation of spatial topology is on the diagonal.
Next, we considered the time-varying expression of patterns across Table 1 Result of parameter optimization on four frequency-band-specific connectivity−configuration matrices. individuals. To measure how strongly a pattern was expressed over time windows, we computed the energy of pattern expression for the temporal coefficients of each pattern for each participant. Similarly, to measure how transiently the pattern was expressed, we also computed the temporal variability of pattern expression using a histogram estimator. We then statistically compared the energy and temporal variability of pattern expression between the group of 114 autism children and the group of 141 typically developing children (Table S1). The development trajectory of all connectivity patterns was evaluated in Table S2.

Homotopic interhemispheric patterns
We first analyzed the homotopic interhemispheric-connection patterns, as they offered a relatively simple yet robust subset of brain connectivity. The homotopic interhemispheric pattern is defined as interhemispheric connectivity between homologous areas. Interhemispheric connections were apparent between the bilateral temporal lobes across four frequency bands (Fig. 2). We observed that in the α-band connectivity pattern (pattern K), the energy at different ages displayed a significant decreasing trend in the ASD group compared to the TD group (see Table S1), indicating an under-connectivity in the homotopic interhemispheric-connection patterns in ASD. In the θ-band connectivity pattern (pattern I), the energy was significantly lower in the ASD group compared to the TD group at the ages of 3 (p < 0.05) and 4 (p < 0.05). However, in the δand β-band connectivity patterns (patterns B and Q), we found no significant differences in energy between the groups (ASD vs TD) across all ages. We then considered the time-varying expression of subnetworks-or subnetwork dynamics-across individuals. In the four frequency-band connectivity patterns, a slightly decreasing trend of temporal variability was observed in the ASD group, but this did not reach statistical significance. The results suggest that the α rhythms are preferentially associated with the interaction processes involving homotopic interhemispheric connections. Fig. 3 illustrates group differences between ASD and TD children in subnetwork energy in different age groups, demonstrating differences in early childhood developmental trajectories of interhemispheric connectivity. In Fig. 3, the homotopic interhemispheric pattern of the ASD group displays a significant trend of increasing network energy with age in the δ band (pattern B; r = 0.24, p = 0.0075) and α band (pattern K; r = 0.18, p = 0.03), whereas the trends were not found to be significant in the TD group (pattern B; r = 0.15, p = 0.06; pattern K; r = 0.12, p = 0.13). In the θ-band-based pattern I and β-band-based pattern Q, the TD group illustrated a weak correlation between subnetwork energy and age (pattern I; r = 0.05, p = 0.5; pattern Q; r = 0.11, p = 0.17); a similar trend was observed in the ASD group.
Taken together, the control group exhibited a faster rate of incremental increases in interhemispheric connectivity during early childhood development compared to that of the ASD group. Particularly, the α-band connectivity exhibited the largest developmental significant difference in terms of long-range interhemispheric connections; ASDrelated abnormalities in the α band have been reported fairly consistently and, thus, may represent a biomarker of ASD (Han et al., 2017). Interesting relationships might be hypothesized linking this potential initial under-connectivity with the early abnormal development of the brain in ASD individuals.

Heterotopic interhemispheric and intrahemispheric patterns
The heterotopic interhemispheric connection was considered as a characteristic distortion of idiosyncratic connectivity. To address this question, we compared heterotopic interhemispheric-connectivity patterns and intrahemispheric-connectivity patterns of autistic and control children across four frequency bands (Fig. 4). The heterotopic interhemispheric pattern is defined as interhemispheric connectivity between heterogeneous areas. Extending our findings of idiosyncratic connectivity, the ASD group demonstrated a consistent reduction in long-range heterotopic interhemispheric-connectivity patterns in comparison to that of the control groups (e.g., patterns A, C, E, N, and U). In δ-band connectivity patterns (pattern A), the energy was significantly lower in ASD group compared to that in the TD group at ages of 3, 4, and 5, whereas the δ-band connectivity pattern C displayed no significant difference between the groups. Similarly, in the θ-band connectivity pattern E, the energy illustrated a significant decreasing trend in the ASD group compared to the TD group. The same phenomenon was evident in terms of the connectivity within the bilateral short-range hemispheric patterns (e.g., patterns L and O). These qualitative assessments of patterns of connectivity reveal a generally heterotopic bilateral and lateralized organization across hemispheres. These results indicate that ASD is characterized by heterotopic interhemispheric patterns.
In addition to the significant difference in connection strength, we also observed a significant difference between the groups in terms of the temporal variability (denoted by the entropy value) of the patterns A, C, E, and U, which represent lateralized long-range interhemispheric connections. Specifically, temporal variability was higher in the control group, suggesting that there was a higher tendency for these patterns to change over time in the controls compared to that in ASD children. Conversely, the patterns of α-band connectivity-including the shortrange, within-frontal lobe (pattern L), within-occipital lobe (pattern O), and prefrontal−occipital junction in both hemispheres (pattern N)-were not significantly different between the two groups and were less prone to change. Fig. 5 depicts group differences between ASD and TD children in subnetwork energy at different ages. Pattern L comprised relatively dense connectivity, including both short-range, within-frontal-lobe and longer, between-frontal-lobe connections. The strength of pattern L followed a weak increasing trend with age in the TD group (r = 0.09, p = 0.28). However, the ASD group reflected an opposite trend with respect to age in the α band (r = −0.07, p = 0.4). Similar trends in connectivity were seen for the pattern O, which comprised relatively dense connectivity that included both short-range, within-occipital-lobe and longer, between-occipital-lobe connections. The strength of pattern O followed a weak increasing trend with age in the TD group (r = 0.03, p = 0.65), whereas the opposite trend was observed in the ASD group (r = −0.04, p = 0.64). Interestingly, the variance in pattern U in the ASD group appears to be nearly zero (approximately 0.1); meanwhile, there was a significance difference (p < 0.001) in variance between groups (ASD vs TD). A significant decreasing trend was observed for the TD group (r = −0.04, p = 0.64). Pattern U appears to be an early brain development network pattern sensitive to typically developing children, which may be useful for tracking development trajectory.

Generalization of short-range connectivity patterns
Thus far, we have established that the interhemispheric-connectivity patterns were more under-connected in the ASD group compared to those of the control group. However, it is also possible that there were common patterns of short-range connectivity in ASD groups across the four frequency bands of connectivity. Indeed, the patterns of over-connectivity within the ASD group have been identified in α-band and β-band connectivities ( Fig. 6; patterns M and S). These maturational patterns appeared to differ in the ASD group from lower frequencies (δ and θ bands) to higher frequencies (α and β bands), as the short-range over-connectivity patterns were only observed in higher frequency bands (α and β bands). An interesting relationship may be hypothesized, linking this potential initial over-connectivity with the early overdevelopment of the brain in ASD (Courchesne et al., 2011(Courchesne et al., , 2007 and the early maturation of white-matter tracts previously reported in young children with ASD (Bashat et al., 2007;Weinstein et al., 2011;Billeci et al., 2012). Some recent evidence suggests that seemingly divergent under-and over-connectivity findings in ASD may be reconciled in a model of reductions, both in network integration (under-connectivity within neurotypical networks) and in network differentiation or segregation (over-connectivity with atypical regions outside neurotypical networks) (O'Reilly et al., 2017). Fig. 7 depicts a linear correlation of network energy with age in patterns M and S. We did not observe a clear developmental trajectory in patterns M and S. Specifically, the TD group illustrated a slight increasing trend of energy with age in pattern M (r = 0.05, p = 0.54), whereas the opposite trend was observed in the ASD group (r = −0.04, p = 0.65).

Discussion
Using novel dynamic machine-learning techniques that decompose matrices into time-varying subnetworks, we studied how dynamic patterns altered and evolved with early childhood development in a group of 115 children with ASD and 141 age-matched typically developing controls. Considerable changes occurred in the patterns of functional brain networks during an early stage of neurodevelopment. We demonstrated that the homotopic interhemispheric patterns of resting-state connectivity in the ASD participants were significantly decreased relative to those in the control children. Moreover, we compared heterotopic interhemispheric-connectivity patterns as well as within-hemisphere connectivity patterns of the autistic and control participants. The ASD group demonstrated a consistent reduction in long-range connectivity in comparison to that of the control group in terms of heterotopic interhemispheric-connectivity patterns. Finally, we observed a pattern of over-connectivity in the ASD group in terms of α-band and β-band connectivities.

Functional connectivity via resting-state EEG
Resting-state brain networks have revealed architectures closely related to underlying anatomical connections (Deco and Corbetta, 2011). The notion of globally coordinated and dynamically competing resting-state networks directly motivates studies that can identify overlapping cognitive systems and explain and predict their Fig. 3. Development of homotopic interhemispheric subnetwork properties of ASD group (red) and TD group (blue) at each frequency band. Scatter plot depicting correlation between age and network energy in each frequency band. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 4. Age-related differences in patterns of heterotopic−interhemispheric connections. The following seven patterns were highlighted as examples of age-related variation in the heterotopic−interhemispheric network across four frequency-band-specific connectivities-patterns A and C for δ-band connectivity; pattern E for θband connectivity; patterns L, N, and O for α-band connectivity; and pattern U for β-band connectivity. Significant differences for each age group in terms of the energy (middle panel) and entropy (right panel) of the four patterns as well as statistical analysis of the two groups without considering age factors (inset) were observed. * indicates network pairs with significantly different connections between groups (*, **, and *** indicate p < 0.05, p < 0.01, and p < 0.001, respectively).
Previous work in spatially and temporally overlapping community detection has uncovered cohesive structures in resting-state brain networks (Fortunato, 2010). In the present study, we employed an NMF approach, which identifies the subnetworks of a brain network that dynamically vary across participants and across time, to obtain the connectivity for each subject in each time window as a nonnegative combination of basis subnetworks (Chai et al., 2017;Eavani et al., 2013). NMF offers prominent advantages over other dynamic graphmodeling approaches. First, the non-negativity constraint means that subnetworks can be interpreted as additive components, which together form the original network. Second, NMF does not make any explicit assumptions (including orthogonality or independence) with respect to the resulting subnetworks, which provides extra flexibility in extracting network components that may correspond to brain processes that overlap spatially (Khambhati et al., 2018c). Third, compared with the hard-partitioning of nodes into discrete modules in dynamic-community detection, NMF pursues a soft-partitioning of the network and allows one to track how functionally interacting brain areas can be dynamically expressed in a continuous and overlapping manner (Khambhati et al., 2018a).
Thus, from the perspective of network-based encoding of the ASDrelated networks, we observed salient connectivity patterns in restingstate EEG data acquired from a group of 115 children with ASD and 141 age-matched typically developing controls.

Abnormal connectivity patterns in ASD
Our study revealed that the ASD group was characterized by shortrange over-connectivity and long-range under-connectivity, particularly in interhemispheric connections. Both aberrant intrahemispheric and interhemispheric connectivity was observed, with a general trend toward under-connectivity, but with probable local short-range overconnectivity as well. Long-range under-connectivity was observed in lower frequency bands, which we hypothesized to be preferentially involved in long-range integrative networks (Figs. 2 and 4), whereas a stronger tendency for over-connectivity was observed (see Fig. 6) in higher frequency bands (hypothesized to be generally associated with more localized processes).
Our present finding regarding the frequency-band-based modulation of over-versus under-connectivity is consistent with that of slower oscillators involving more neurons in larger volumes (Von Stein and Fig. 5. Development of heterotopic interhemispheric and intrahemispheric patterns of ASD group (red) and TD group (blue). Scatter plot illustrating correlation between age and network energy. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) T. Zhou, et al. NeuroImage: Clinical 26 (2020) 102251 Sarnthein, 2000). Indeed, coordination of faster frequency required lager axons for any particular distance between two regions. Based on this fact, sharing a short range of high-frequency-band activity is more efficient for bias connectivity, whereas long-range interactions are required for lower frequency bands. Several studies have recently reported that high EEG frequencies are preferentially associated with more localized neural activity, whereas slower rhythms (δ, θ) are preferentially associated with more widespread integrative cross-region Fig. 6. Age-related differences in generalized patterns of localized connections. Following two patterns are highlighted as examples of age-related variation in localized network across four frequency-specific bands of connectivity: pattern M for α-band connectivity and pattern S for β-band connectivity. Significant differences for each age group in terms of the energy (middle panel) and entropy (right panel) of the four patterns as well as statistical analysis of the two groups without considering age factors (inset) were observed. * indicates network pairs with significantly different connections between groups (*, **, and *** indicate p < 0.05, p < 0.01, and p < 0.001, respectively). T. Zhou, et al. NeuroImage: Clinical 26 (2020) 102251 interaction processes that involve long-range connections. The synchronization of short-range cortical networks tend to be associated with higher frequencies (α, β) (Uhlhaas and Singer, 2010;. These recent findings are consistent with our present results.

Neurodevelopmental changes in salient connectivity patterns
Considerable changes in the patterns of functional brain networks occur throughout neurodevelopment. The modularity of brain structures evolve with age, as does inter-regional connectivity (Meunier et al., 2009). Several studies (Teipel et al., 2009) have reported a correlation between EEG α-band coherence and the structural integrity of white matter and the fact that the connections among resting-state networks strengthen with age (Betzel et al., 2014). Consistent with these findings, we observed age-related differences in the strength and flexibility of the α-band patterns between autistic and control children. Specifically, the α-band patterns of interhemispheric homotopic long-range connections (pattern K) displayed the largest developmental differences. Interesting relationships may be hypothesized between the α-band interhemispheric long-range connections and neurodevelopment to link this pattern with ASD-related idiosyncratic early childhood development markers.
Group differences were observed in the developmental trajectories of specific brain connectivity patterns. With respect to the brain connectivity pattern that may display the largest developmental differences, our findings demonstrated that interhemispheric connectivity patterns are particularly vulnerable to ASD. In general, our findings demonstrated an age-dependent ASD connectivity pattern with interhemispheric long-range under-connectivity and local short-range overconnectivity. This observation is supported by findings from several past studies-a delay in the onset of axonal remodeling could drive brain overgrowth; thus, short-distance connections are favored over long-distance connections via typical developmental mechanisms (Lewis et al., 2013;Lewis and Elman, 2008).

Limitations
The present work had some limitations. First, we analyzed only sensor-level signals, which suffer from the effects of spatial leakage. In addition, the estimation of connectivity between sensors may be difficult to interpret, which prompted us to adopt source-level connectivity with the NMF algorithm. There are several advantages of source space projection in connectivity assessment (Schoffelen and Gross, 2009). First, the results can be overlaid directly onto structural brain images, enabling direct interpretation with respect to underlying anatomy. Second, source localization (via adaptive techniques such as beamforming) reduces artifacts in M/EEG data, indicating that the signal-tonoise ratio (SNR) of projected data is higher than the SNR of raw data in channel space (Lai et al., 2018). This suggests to us that the combination of source-level connectivity and NMF could be complementary and provides a tool for analyses of cortical communication patterns from the EEG.
Moreover, the ICA-based artifact-rejection technique that we used may have distorted the network topology we obtained. ICA is widely regarded as a robust and powerful method for artifact removal from EEG data and has recently received increasing attention (Castellanos and Makarov, 2006;Zafar et al., 2019). However, during subtraction (ICA is a linear operation applied to all sensors) of an artifactual waveform from segments of artifactual data, several challenges are encountered: (a) it is unclear what part of the artifactual waveform is artifact and what part is signal, (b) ICA does not completely remove the entire signal from any channel, even if the entire signal was artifactual, (c) ICA contaminates other channels with the artifactual component, and (d) ICA applies the same transformation throughout the recording and, therefore, mixes the non-artifactual time segments using the same spatial coefficients computed for that component. Therefore, ICA should only be applied to artifacts such as eyeblink, muscle, and cardiac rhythms that persist throughout the recording. Future studies should elucidate the effects of ICA artifact rejection on subsequent connectivity and network analyses.

Conclusions
We applied a novel extracting-dynamic-network framework for uncovering a set of interpretable patterns that highlighted significant group-level differences between ASD and control children. Our results support the framework that ASD is characterized by a general trend toward a long-range under-connectivity of lower band frequencies, particularly for interhemispheric connections, combined with shortrange over-connectivity. These results provide a context for more clearly understanding how early childhood developmental factors induce ASD-related abnormalities in brain connectivity throughout early childhood development.
Few studies have investigated individual variability within ASD and neurodevelopmental status. Hence, future studies must pay more attention to individually distinct (idiosyncratic) distortions in participants with ASD. Further, the EEG analyses of ASD individuals typically involve the computation of spectral power over predefined frequency bands (i.e., 2-4 Hz for θ, 4-8 Hz for δ, 8-13 Hz for α, and 13-30 Hz for β bands). In fact, these frequency bands change over time and among participants, particularly in children. In future work, individual frequency bands must be considered. Moreover, our NMF approach may provide invaluable insights on a variety of other neurodevelopmental disorders such as schizophrenia and attention-deficit hyperactivity disorder.

Declaration of Competing Interest
The authors state that they have no conflict of interest.