BOLD signal variability and complexity in children and adolescents with and without autism spectrum disorder

Highlights • Resting-state BOLD signal variability and complexity were examined.• No significant group differences were observed in youth with and without autism.• A continuum of brain-behavior relationships was observed across diagnostic groups.• Positive correlations were found between brain measures, age and global efficiency.• Negative correlations were found between the brain measures and behavioral severity.


Introduction
Two key components of healthy brain functioning are variability of neural signaling and complexity of these signals. In the healthy brain, variability of neural signaling allows for the formation of functional networks (Fuchs et al., 2007) and the exploration of multiple stable functional states (Ghosh et al., 2008;McIntosh et al., 2008McIntosh et al., , 2010. Previous work has shown that higher variability of brain signals is associated with better behavioral performance (Garrett et al., 2011). Further, variability changes in response to task demands: BOLD variability has been shown to be lowest in the resting-state, increased during internally focused tasks, and highest during externally focused tasks (Grady and Garrett, 2018). It has been suggested that variability allows for greater environmental uncertainty during externally-directed compared to internally-directed tasks, and variable, flexible neural signaling allows the brain to adapt to such uncertainty (Grady and Garrett, 2018).
Several studies have examined age related changes in brain signal variability. Using fMRI, Garrett et al. (2010) found that during fixation periods of a task, the standard deviation of BOLD time series increased with age in 33% of voxels and decreased in 67% of voxels. The spatial pattern of age-related changes in variability was different than that of mean BOLD activity, suggesting that these metrics provide unique information about age-related changes in brain activity. Variability has also been examined using the mean square successive difference (MSSD;von Neumann et al., 1941), defined as the average of the sum of squared differences in amplitude between successive time points. Thus, MSSD measures variability from one time point to the next. MSSD reflects variability in a signal that is independent of drifts in the mean (Garrett et al., 2011). Using fMRI, Nomi et al. (2017) found that MSSD decreased from ages 6 to 85 in the majority of brain regions examined. In attention deficit hyperactivity disorder (ADHD), positive correlations have been reported between symptom severity and MSSD in default mode regions (Nomi et al., 2018).
Signal complexity is also important for optimal brain function. One way to measure the complexity of a signal is with sample entropy, which assigns high values to more complex signals, and low values to highly deterministic or random signals (Costa et al., 2005). It is important to measure the complexity of brain signals in addition to variability, because a variable brain signal may not necessarily be complex. A signal with higher entropy can be interpreted as having higher information processing capacity (Gatlin, 1972;Heisz and McIntosh, 2013;Shannon, 1948). Like variability, complexity is thought to reflect the ability of the brain to adapt to unpredictable environments (Goldberger et al., 2002). Using fMRI, it has been shown that signal complexity in various regions of the default mode network was positively correlated with multiple cognitive functions, including attention, language, and short-term memory . Higher multiscale entropy of EEG signals has been associated with greater knowledge representation (Heisz et al., 2012). It has further been shown that entropy of EEG signals increases across development (Lippé et al., 2009;McIntosh et al., 2008;Misic et al., 2010).
Abnormal levels of variability and complexity in the brain may be related to sub-optimal cognition. Too much variability can result in inefficient information processing and ineffective exploration of different network configurations in the brain (Ghosh et al., 2008;McIntosh et al., 2008McIntosh et al., , 2010. Autism spectrum disorder (ASD), a neurodevelopmental disorder that is characterized by atypical social communication and restricted, repetitive and stereotyped behaviors (American Psychiatric Association, 2013), is one condition that may be characterized by detrimental levels of noise. ASD is hypothesized to be characterized by an increased ratio of excitatory to inhibitory coupling in the brain, which could be associated with noisier and less stable signaling, reduced functional differentiation, and less efficient information processing (Rubenstein and Merzenich, 2003). These alterations may relate to abnormal myelination and synaptic development (Rubenstein and Merzenich, 2003). Inhibitory signaling is believed to be associated with improving the specificity of excitatory brain signals; thus, imprecise brain activity and increased noise can result from increased excitatory and decreased inhibitory signaling (Nelson and Valakh, 2015).
Studies of dynamic functional connectivity (FC) have shown that FC is also more variable in ASD (e.g. Chen et al., 2017;Falahpour et al., 2016;Wee et al., 2016;Zhang et al., 2016). Chen et al. (2017) reported greater variance of distributed functional connections in ASD, which was related to total scores on the Autism Diagnostic Observation Schedule (ADOS). Dynamic FC has been shown to improve classification of ASD compared to typically developing (TD) participants (Wee et al., 2016), and to be a better predictor of ASD behaviours compared to "static" FC measured over an entire time series (Jia et al., 2014). However, the nature of resting-state BOLD signal variability in ASD compared to TD individuals remains unclear. Entropy of EEG signals has also been shown to be reduced in ASD at rest (Bosl et al., 2011) and during both social and non-social tasks (Catarino et al., 2011), but the nature of the complexity of resting-state BOLD time series in ASD is currently unknown. Further, relationships between information processing capacity in functional compared to structural networks in ASD are not well-characterized. Anatomy plays a role in shaping functional networks in the brain (Deco et al., 2013;Honey et al., 2007Honey et al., , 2009Honey et al., , 2010. For instance, regions that exhibit high structural connectivity (SC) typically exhibit strong FC (Honey et al., 2007(Honey et al., , 2009. Previous work has revealed relationships between information processing capacity in structural networks and cognitive functioning (Berlot et al., 2016;Reijmer et al., 2013). Diffusion weighted imaging (DWI) can be used to estimate the strength of white matter connections between brain regions via tractography. As time series are not available for structural connectomes, a graph theory metric, global efficiency (GE), can be used to estimate information processing capacity (Bullmore and Sporns, 2012).
In the present study, we characterized BOLD signal variability and complexity in youth with and without ASD. Further, we characterized relationships between these measures and age, GE, IQ, and behavioral severity.

Participants
Twenty male participants with ASD (M = 13.25 years, SD = 2.87 years) and 17 male TD participants (M = 13.42 years, SD = 3.21 years) from the San Diego State University sample from the Autism Brain Imaging Data Exchange (ABIDE) II database  were included in this study. Informed consent or assent had been obtained for all participants and caregivers in accordance with the University of California, San Diego and San Diego State University Institutional Review Boards. Participants were matched for age, full-scale IQ and head motion (Table 1). Participants were excluded if they were less than 8 years of age, their full-scale IQ was below 75 or if their head motion during the scan exceeded a mean framewise displacement (FD) 0.2 mm. Further, participants were included only if they had an MPRAGE scan, resting-state fMRI scan, and DWI scan.

Resting-state fMRI preprocessing
Data were preprocessed using the Optimization of Preprocessing Pipelines for NeuroImaging (OPPNI) software (Churchill et al., 2012a(Churchill et al., ,b, 2015. First, motion correction was performed using AFNI's 3dvolreg function, followed by the generation of subject-specific non-neuronal tissue masks using the PHYCAA + algorithm (Churchill and Strother, 2013), then replacement of outlier volumes with interpolated values from neighbouring time points (Campbell et al. (2013); https://www. nitrc.org/projects/spikecor_fmri). This motion censoring step was implemented due to previous evidence suggesting that censoring is highly effective at removing residual motion-related artifacts in resting-state fMRI data. For instance, Ciric et al. (2017) noted that censoring mitigates motion artifacts as well as distance-dependent artifacts. They also found that methods that are less effective at denoising are not able to identify modular network structure in functional connectomes. They noted that global signal regression can mitigate motion artifacts, but has the drawback of introducing distance-dependent artifacts. Thus, in this study, censoring was used to remove residual motion artifacts. Importantly, the maximum number of time points that were censored was 7 out of 175, and the number of censored time points did not differ between ASD and TD groups. Nonetheless, we analyzed the effects of censoring on MSSD and entropy estimates. These analyses are described in the Supplementary Material. These steps were followed by slice-timing correction and temporal detrending using a second-order Legendre polynomial. Next, principal component analysis was performed on the motion parameters obtained from the motion correction step. Principal components that accounted for more than 85% of the variance of the motion parameters were regressed out of the fMRI data. The time series of the mean white matter and CSF signals were regressed out of the data, and finally, lowpass filtering was performed with a cutoff of 0.1 Hz.
The ROI atlas used in this study is described in Bezgin et al. (2012). Parcellation of each hemisphere into 48 regions was performed on a macaque brain surface using mapping rules established by Kötter and Wanke (2005). The parcellation was then transformed to the MNI human brain using landmarks from Caret (www.nitrc.org/projects/ caret; Van Essen et al., 2001;Van Essen and Dierker, 2007). As direct anatomical connections have been measured in the macaque brain, using this parcellation scheme helps to eliminate false positive connections found by probabilistic tractography used in human DWI analyses. This atlas is the default atlas used in The Virtual Brain (Sanz Leon et al., 2013) to simulate brain network dynamics from human neuroimaging data. Further, as noted in Bezgin et al. (2012), "although the detailed anatomy of the human brain is still poorly understood, there exist many homologies with the brains of non-human primates. Thus, a viable route to enhance the functional-anatomical framework in humans would be to map human brain regions to areas that have been well studied in invasive monkey experiments". For this study, we only used cortical regions from the atlas, resulting in 82 ROIs in total (41 per hemisphere; Table 2; Fig. 1). The atlas was transformed from MNI space to each subject's T1 space, then subsequently transformed to each subject's functional space using the inverse of each subject's functionalto-anatomical transform, which was obtained using linear registration with 6°of freedom (DOF). Next, the time series of each region was extracted using the mri_segstats function in Freesurfer.

DWI preprocessing
Preprocessing of dwMRI data included motion correction using eddy current correction in FSL, fitting diffusion tensor models at each voxel using the dtifit function in FSL, co-registering each subject's skullstripped T1 MRI to DTI space using linear registration with 12°of freedom (DOF), registering a standard MNI T1 image to each subject's T1 image, labeling the grey matter in T1 space with the 82 cortical regions using nearest neighbour interpolation, then transforming the ROI-labelled T1 image into diffusion space. BEDPOSTX (Bayesian Estimation of Diffusion Parameters Obtained using Sampling Techniques; the X refers to the modeling of crossing fibres) in FSL was used to fit the probabilistic diffusion model on the preprocessed dwMRI data. Markov Chain Monte Carlo sampling is run to build up distributions on diffusion parameters at each voxel. Probabilistic fiber tracking using the probtrackx2 function in FSL was then performed to define weights (fiber counts/number of streamlines) and anatomical distances between each pair of ROIs. Probtrackx2 takes repetitive samples from the distributions of voxel-wise principal diffusion directions. A streamline is computed through each of these local samples, thus generating a probabilistic streamline. By taking many samples, a histogram of the posterior distribution of the streamline location, or connectivity distribution, is created. All masks for tractography were interface masks, that is, masks of the boundary between the gray matter and white matter. This approach is referred to as anatomically-constrained tractography, and allows one to avoid the white matter seeding bias, that is, the tendency for major white matter structures to be over-defined (Smith et al., 2012). The seeds for tractography were masks of each of the 82 ROIs. The following settings were used: 5000 samples, 2000 steps per sample, step length of 0.5, curvature threshold of 0.2, loop check on, and correction of the path distribution for the length of the pathways. Targets were masks for each of the other 81 ROIs; streamlines were terminated once they reached the target mask. For a given ROI, the exclusion mask for connections with ipsilateral ROIs consisted of that ROI and ROIs in the contralateral hemisphere, and the exclusion mask for connections with contralateral ROIs was the mask for that ROI. Any pathways that entered the exclusion masks were discarded. For each subject, an SC weights matrix was defined for each pair of ROIs as the number of streamlines detected between the two ROIs divided by the total number of samples. A tract lengths matrix was defined for each pair of ROIs as the average distance between the two ROIs. Finally, both the weights and tract lengths matrices were thresholded based on tractography data for the same atlas from the Co-CoMac database (Stephan et al., 2001), such that non-existent connections in the CoCoMac data were set to 0 in the human weights and tract lengths matrices to control for false positives. Next, the GE of each participant's SC weights matrix was calculated. GE was used as a summary measure of structural networks because a single value is calculated for the entire network, instead of one measure for each ROI. GE is defined as the average inverse shortest path length in the brain network, and is a measure of the overall capacity of the brain network to perform "parallel information transfer and integrated processing" (Bullmore and Sporns, 2012). Therefore, it would be expected that brain networks that exhibit high GE in their structural networks would exhibit high overall entropy of functional networks. GE was calculated for each participant's SC matrix using the efficiency_wei.m function from the Brain Connectivity Toolbox (Rubinov and Sporns, 2010). We used a graph metric to estimate information processing capacity in each participant's structural network, since time series are not available for estimating variability or entropy.
Based on previous work examining correlations between GE and age that controlled for effects of head motion (Rudie et al., 2012), we examined the correlation between GE and mean FD that was calculated for the DWI scans, and found that this relationship was significant, r(35) = −0.39, p = 0.02. Therefore, the effects of head motion were regressed out of the GE values.

BOLD signal variability and complexity
Prior to calculating variability and complexity, the BOLD time series for each ROI was normalized to have a mean of 0 and standard deviation of 1. BOLD signal variability was defined using MSSD (von Neumann et al., 1941), which is defined as the sum of the squared differences of BOLD signal values between successive time points for a region of interest, divided by the number of time points minus one. A primary benefit of MSSD is that it prevents overestimates of dispersion in data that result from a shift in the mean (Garrett et al., 2011). The formula for MSSD is as follows: Complexity was defined as the sample entropy of the signal using the sample entropy function from PhysioNet (https://physionet.org/ physiotools/mse/). Sample entropy is used to determine the appearance of repetitive patterns in a time series, thus measuring the regularity of the time series. As such, sample entropy values are low for signals that are completely random or strongly deterministic, and are high for more complex signals. (Costa et al., 2005;Pincus, 1991;Richman and Moorman, 2000).
One issue with measuring sample entropy is that two parameters must be selected: the pattern length m, which is the number of data points used for pattern matching, and the tolerance factor r, which is the fraction of the standard deviation of the signal. Recently, Yang et al. (2018) presented a method for selecting optimal m and r values for an fMRI dataset based on the standard error of entropy estimates in CSF voxels. For different combinations of m and r, entropy was calculated for each voxel in each participant's CSF mask. Then, the standard error of this entropy estimate was defined as = SE 1.96* ( ) 2 σSampEn SampEn as in Yang et al. (2018). The median standard error was then calculated across participants, for each combination of m and r. Further, the number of sample entropy estimation failures due to an absence of pattern matches was summed across all CSF voxels and participants. "Acceptable" combinations of m and r were defined as those with a median standard error of less than 0.1, and with zero sample entropy estimation failures across all CSF voxels and participants. These values were used to calculate the sample entropy of the ROIs in the grey matter. Subsequently, for each participant, sample entropy values were averaged across the acceptable m and r combinations for each ROI.
Using this method, we tested m values of 1 to 4 with a step size of 1, and r values ranging from 0.05 to 0.80, with a step size of 0.05, which were the same r values used in Yang et al. (2018). These investigators also tested greater values of m, but found that for m > 4, the standard error was greater than 0.1 for all values of r; thus, we only test m values of 4 or less.  Lobaugh, 2004) is a multivariate statistical method that is used to determine optimal relationships between a set of brain variables and either a study design (mean-centering PLS) or a set of behaviour variables (behavioural PLS). In PLS, singular value decomposition is used to calculate orthogonal patterns that explain the maximal covariance between brain variables and design or behaviour variables. For each pattern, "brain saliences" are calculated for each brain region to indicate how strongly each region contributes to the relationship between brain and design/behaviour. In other words, brain saliences indicate which brain regions best characterize the relationship. In mean-centering PLS, design saliences indicate the group, condition, or group x condition profiles that best characterize the relationship between the brain and design variables. In the case of behavioural PLS, behaviour saliences indicate the profile of behaviour variables that best characterize the relationship between the set of brain variables and behaviour variables. "Brain scores" and "behaviour scores" are calculated for each pattern of relationships between brain and behaviour variables, which indicate the contribution of each participant's brain and behaviour variables, respectively, to the pattern. Brain scores are calculated by multiplying the matrix of brain variables by the brain saliences; behaviour scores are calculated by multiplying the matrix of behaviour variables by the behaviour saliences. Further, singular values show the proportion of covariance that each pattern (relationship between brain and design/behaviour) accounts for. In this study, mean-centering PLS was used to determine optimal contrasts in entropy between conditions (acceptable combinations of m and r), and optimal contrasts in MSSD and entropy between ASD and TD groups. Behavioural PLS was used to determine relationships each of these brain variables and a set of predictor variables, including GE of the structural networks, age, IQ, and scores on the Social Responsiveness Scale (SRS; Constantino and Gruber, 2005). The SRS is a parent or teacher report of ASD-related traits that was created for use in the general population in both clinical and educational settings. Behavioural PLS was also performed using scores on the Autism Diagnostic Observation Schedule 2 (ADOS2) instead of SRS scores in the ASD group. This analysis was only performed for ASD participants, as ADOS2 scores were not available for the TD group. Each behaviour variable was normalized to have a mean of 0 and standard deviation of 1.
The significance of each PLS pattern can be determined using permutation testing. The rows (participants) of the data matrix were reshuffled and the singular value was recalculated. This procedure was repeated 1000 times to obtain a distribution of singular values. Then, a p-value for the original singular value was obtained by calculating the proportion of singular values from the sampling distribution that are greater than the original singular value.
Further, the reliability of each brain salience can be determined using a bootstrapping procedure. Here, 500 bootstrap samples were generated by randomly sampling participants with replacement while maintaining group membership. Next, for each brain salience (brain region), a bootstrap ratio (BSR) was calculated as the ratio of the brain salience to the standard error of the salience from the bootstrap samples, which is a measure of the stability of each brain salience regardless of which participants are included in the analysis. Stable brain regions were defined as those that surpassed a BSR threshold of +2, which corresponds to approximately a 95% confidence interval.

Data visualization
Brain plots were created using the following Python packages: matplotlib (Hunter, 2007), nibabel (Brett et al., 2016), and nilearn (Abraham et al., 2014). All other figures were created in MATLAB (version R2015b). Fig. 2 shows the results of the entropy estimates in CSF voxels. This analysis revealed that for m = 1, all r values were acceptable. For m = 2, r values from 0.20 to 0.65 were acceptable. Hereafter, the term "condition" will be used to refer to an acceptable combination of m and r.

Effects of m and r on entropy estimates
To ensure that any group differences in MSSD and entropy in the gray matter ROIs were not confounded by residual effects of head motion, we first performed a behavioral PLS using mean FD as the "behaviour" variable. We found that there was no significant relationship with head motion for MSSD (p = 0.53) or entropy (p = 0.23). However, we found that the relationship between SD and head motion was significant (p = 0.002), even when additional ICA denoising was implemented using ICA-AROMA (p = 0.005; Pruim et al., 2015aPruim et al., , 2015b. These findings therefore provide more support for the use of MSSD over SD for measuring brain signal variability. Mean-centering PLS was then used to examine how different acceptable combinations of m and r affected entropy estimates in the 82 grey matter ROIs across all participants. There was one significant pattern from this analysis that showed a contrast between different combinations of m and r (p < 0.001, 99.97% covariance explained; Fig. 3). Therefore, entropy estimates differ significantly depending on which combinations of m and r are used for the analyses. For subsequent analyses, since there were differences in entropy that depended on the choice of m and r, we averaged entropy estimates across all acceptable m = 2 conditions except for r = 0.20, which showed a contrast with the other r values for m = 2.

Categorical analyses of MSSD and entropy using mean-centering PLS
Next, mean-centering PLS was performed to determine differences in MSSD and entropy between diagnostic groups. In other words, a categorical approach was used to analyze MSSD and entropy. For MSSD, no significant differences were observed between ASD and TD groups (p = 0.18). Similarly, groups were not significantly different in entropy (p = 0.15).

Dimensional analyses of MSSD and entropy using behavioral PLS
Next, behavioral PLS was performed to characterize relationships between MSSD and entropy in each brain region and the following measures: GE, age, IQ, and SRS scores, and between entropy and these measures. One main advantage of using PLS in this way is that relationships between different "behaviour" variables in relation to brain variables can be analyzed, for instance, to determine if two behaviour variables differ in terms of the strength or direction of their linear relationships with the brain variables. PLS was performed using SRS total scores as opposed to scores on the five SRS subscales, because scores on the subscales were highly correlated with each other as well as total SRS scores (r > 0.80) and therefore biased the PLS results.
The behavioral PLS analyses revealed one significant pattern for MSSD (p = 0.02, 58.11% covariance explained, Fig. 4) and for one significant pattern for entropy (p = 0.01, 48.88% covariance explained, Fig. 5). For both analyses, there was a distributed set of brain regions that exhibited positive correlations with GE and IQ, but negative correlations with SRS scores. For MSSD, 19 out of 82 (23.17% of) brain regions contributed reliably to the brain-behaviour pattern. For entropy, 17 out of 82 (20.73% of) brain regions contributed reliably to the pattern. These regions are listed in Supplementary Table 1. For both MSSD and entropy, a continuum of brain-behaviour relationships can be observed across ASD and TD participants, as shown in Figs. 4B and 5B. In other words, a broad range of brain and behaviour scores can be observed across all participants, but overall, brain and behaviour scores were higher in the TD group (MSSD brain scores: t(35) = −1.96, p = 0.06; MSSD behaviour scores: t(35) = −3.03, p = 0.005; entropy brain scores: t(35) = 2.40, p = 0.02; entropy behaviour scores: t (35) = 3.49, p = 0.001). Further, for both MSSD and entropy, behaviour scores were more variable for the ASD group compared to the TD group, which illustrates the highly heterogeneous nature of ASD, although the group difference did not reach significance for MSSD (MSSD, SD ASD = 1.18, and SD TD = 0.94, F(19, 16) = 2.29, p = 0.10; for entropy, SD ASD = 1.23, SD TD = 0.93, F(19, 16) = 3.64, p = 0.01).
Notably, as shown in Table 3, although the correlations between brain and behaviour scores were higher for the ASD group compared to the TD group, the correlations did not differ significantly between the diagnostic groups. The brain saliences for the MSSD and entropy analyses were significantly correlated, r = 0.46, p < 0.001 (1000 permutations), showing that there was a strong relationship between the brain regions that contributed to the MSSD pattern and the regions that contributed to the entropy pattern (Fig. 6).
Behavioral PLS analysis was also performed in the ASD group using ADOS2 scores in addition to the other predictor variables. As ADOS2 scores were only available for the ASD group, the TD group was not included in these analyses. When ADOS2 total scores were included, the PLS analysis was significant for MSSD (p = 0.045), but not for entropy (p = 0.18). However, in both cases, the behaviour saliences were not reliable for ADOS2 total scores, as the 95% CIs crossed the x-axis. A similar pattern was observed when the ADOS2 social affect (SA) and RRB scores were used instead of the total scores: PLS was significant for MSSD (p = 0.02) and not significant for entropy (p = 0.15), but the behaviour saliences for ADOS2 SA and RRB scores were not reliable.

Interactions between MSSD-behaviour and entropy-behaviour correlations
As brain-behaviour relationships were similar for both MSSD and entropy, we combined both of these brain measures in a single behavioural PLS analysis, with each brain measure being treated as a separate "condition", to determine if any additional patterns existed that showed a contrast between MSSD and entropy in terms of correlations with the predictor variables. However, the results of the analysis showed that there was one significant pattern (p < 0.001) showing the same pattern of brain-behaviour correlations as the individual analyses (i.e. positive correlations between brain variables and age and brain variables and GE, and negative correlations between brain variables and SRS scores). No additional significant patterns showing different Fig. 3. Contrast between entropy conditions and associated BSRs for each region, at a threshold of +2. Error bars show 95% confidence intervals determined through bootstrap resampling.
brain-behaviour correlations for MSSD and entropy were observed.

Correlations between behaviour measures
Finally, the relationships between the set of predictor variables were analyzed across all participants and within groups. As shown in Fig. 7, there was a moderate positive correlation between GE and age, a weak negative correlation between GE and IQ, and a moderate negative correlation between GE and SRS scores. SRS scores were weakly negatively correlated with IQ. Further, GE did not differ significantly between the ASD and TD groups, t(35) = −1.18, p = 0.25. A notable difference between groups was the correlation between age and GE: in ASD, this correlation was positive (r(18) = 0.46, p = 0.04), whereas the correlation was negative, but not significant, in controls (r(15) = Fig. 4. Behavioural PLS results for MSSD. A) Contrast in relationships for correlations between MSSD and predictor variables, B) associated brain and behaviour scores for each group, C) BSRs for each region, D) brain plot of BSRs. Regions with a BSR surpassing a threshold of +2 are shown in orange. Error bars show 95% confidence intervals determined through bootstrap resampling. Blue circles = ASD, red circles = TD. Fig. 5. Behavioural PLS results for entropy. A) Contrast in relationships for correlations between entropy and predictor variables, B) associated brain and behaviour scores for each group, C) BSRs for each region, D) brain plot of BSRs. Regions with a BSR surpassing a threshold of +2 are shown in orange. Error bars show 95% confidence intervals determined through bootstrap resampling. Blue circles = ASD, red circles = TD.
The correlation matrix for the predictor variables used in the ASDonly PLS analysis (i.e. including ADOS scores) is shown in Supplementary Fig. 4.

Overview
Variability and complexity of resting-state BOLD signals were examined in participants with and without ASD. For both MSSD and entropy, a continuum of brain-behaviour relationships was observed across diagnostic groups despite a lack of significant differences between groups. A distributed set of brain regions exhibited positive correlations between MSSD and GE and age in both groups, and negative correlations with SRS scores. A similar pattern was observed for entropy.

Variability and complexity in ASD
Categorical and dimensional approaches were used to characterize MSSD and entropy in participants with ASD compared to TD participants. The categorical approach involved analyzing differences between diagnostic groups (ASD and controls), whereas the dimensional approach involved the use of continuous measures: age, IQ, SRS scores, and GE. No group differences were observed when a categorical approach was used; however, using the dimensional approach, a set of brain regions exhibited relationships with the predictor variables.
The negative relationship between entropy and SRS scores supports the notion that greater severity of ASD behaviours is associated with decreased entropy, and therefore decreased information processing capacity, in the brain. Few studies have examined entropy in ASD; however, our results are in line with a previous EEG study that found reduced MSE in ASD during both social and non-social tasks (Catarino et al., 2011). Reduced entropy in participants with more severe ASD behaviors supports the "loss of brain complexity hypothesis" proposed by , which states that the complexity of neural signals reflects the capability of the brain to adapt to changing environments, and that in pathological states, this "adaptive capacity" of the brain is reduced.
Similarly, a continuum of brain-behaviour relationships was observed for MSSD, whereby SRS scores were negatively correlated with MSSD. A dimensional approach has been used to study MSSD in ADHD. Nomi et al. (2018) did not find significant differences in MSSD between children with and without ADHD, but found positive correlations between MSSD and ADHD behavioural severity across diagnostic categories. While Nomi et al. (2018) found positive correlations between MSSD and symptom severity, we found negative correlations between these measures. This finding may seem counterintuitive to theories suggesting that ASD is characterized by an increased ratio of excitatoryinhibitory coupling in the brain, which may lead to poor functional differentiation and noisier and unstable neural signaling, resulting in inefficient information processing (Rubenstein and Merzenich, 2003). One possible explanation for this discrepancy is that detrimentally increased levels of noise may manifest at smaller scales in ASD, but not at macroscopic scales as measured by fMRI. This hypothesis is in line with the theory that physical connectivity at the level of microcolumns is increased in ASD, but computational connectivity is reduced (Belmonte et al., 2004). In other words, within neural assemblies, connections between synapses and fiber tracts are increased, but long-range connections between functional brain regions are reduced. The authors hypothesized that increased physical connectivity could lead to undifferentiated neural regions, which would preclude the development of effective communication between long-range functional regions. Thus, at macroscopic scales, ASD-like symptoms may be associated with reduced variability, which could be associated with a reduced capacity to explore different functional network configurations. Further, we included a larger set of predictor variables in our study compared to Nomi et al.'s study. In addition to behavioral severity, we included age and a Table 3 Correlations between brain and behaviour scores.  measure of information processing in structural networks (GE). This analysis revealed important interactions between GE, age, and behavioral severity, whereby for both MSSD and entropy, a distributed set of brain regions exhibited positive correlations with GE and age, but negative correlations with behavioral severity. Thus, factors such as the efficiency of structural networks and age may be modulating the relationship between MSSD, entropy and measures of behavioral severity. The PLS analyses used in this study provided a data-driven, "systems" approach for understanding BOLD signal variability and complexity in ASD. For both MSSD and entropy, a set of distributed brain regions contributed reliably to the brain-behaviour patterns in the behavioural PLS analyses. Previous studies have also reported widespread alterations in measures of brain function in ASD. For instance, various studies have reported atypical FC in multiple resting-state networks (e.g. Abraham et al., 2017;Anderson et al., 2011;Rashid et al., 2018;Rudie et al., 2012;Supekar et al., 2013;Wee et al., 2016). Further, distributed abnormalities in white and grey matter development have been reported in ASD (Carper et al., 2002;Sparks et al., 2002). Müller (2007) suggested that ASD should be characterized as a "distributed disorder", or in other words, a network disorder, which is applicable to genetics, brain structure and function, and behaviour. It was noted that the many potential genetic risk factors for ASD can influence the development of multiple networks in the brain; therefore, "localizing models" may not be helpful for understanding ASD. As stated by Müller (2007), it is likely that there is a "distributed ontogenetic starting point affecting several emerging brain regions and functional systems… and each of these will in turn affect additional regions and functional systems throughout development". Therefore, data-driven, whole-brain analyses are important for understanding brain function in ASD.

Variability and entropy increase with age
Distributed brain regions showed increases in MSSD and entropy from childhood through adolescence, suggesting that across development, information processing capacity increases in the brain (Figs. 5 and 6). Previous work has suggested that increased information processing capacity over development results from a greater quantity of possible configurations of functional networks . Entropy of EEG signals has been shown to increase from childhood to adulthood, and is associated with higher accuracy and less variable reaction times during task performance (Lippé et al., 2009;McIntosh et al., 2008;Misic et al., 2010). Increases in entropy across development are thought to reflect increased integration and segregation in the brain (Garrett et al., 2013;McIntosh et al., 2008).
While developmental changes in variability and entropy have mostly been studied using EEG, Nomi et al. (2017) used fMRI to study changes in variability across the lifespan. In this study, central executive, default mode, visual, sensorimotor, and subcortical regions exhibited linear decreases in MSSD from ages 6 to 85, whereas nodes in the salience network and bilateral ventral temporal cortices exhibited increases in variability. The discrepancies between these findings and ours may relate, in part, to differences in age ranges and preprocessing strategies, which can affect estimates of BOLD signal variability (Garrett et al., , 2013. Additional research on developmental changes in BOLD signal variability across different age ranges using various preprocessing strategies is required to resolve these discrepancies.

Variability and entropy are associated with global efficiency in structural networks
Behavioural PLS also revealed positive correlations between MSSD of BOLD time series and GE in structural networks, and between entropy and GE (Figs. 5 and 6). This finding supports the notion that structure shapes function in the brain (Deco et al., 2013;Honey et al., 2007Honey et al., , 2009Honey et al., , 2010. Further, there was a contrast in the relationship between GE and behavioural severity, such that brain regions that exhibited positive correlations with GE exhibited negative correlations with behavioural severity. Previous work has shown relationships between GE in structural networks and cognitive abilities. Berlot et al. (2016) showed that mild cognitive impairment (MCI) patients exhibited reduced GE of brain networks, which was associated with reduced cognitive control. Further, in Alzheimer's disease, reduced GE has been associated with memory and executive function abilities (Reijmer et al., 2013). Collectively, these results suggest an important relationship information processing capacity in structural and functional networks and cognitive functioning.

Correlations between GE and age
While groups did not differ in GE, a significant positive correlation was found between GE and age in the ASD group, yet the correlation between GE and age in the TD group was negative and non-significant. A previous study of similarly-aged ASD and TD participants also did not find group differences in GE; however, they found that GE in structural networks increased with age (controlling for head motion) in TD, but not ASD, participants (Rudie et al., 2012). Another study reported a positive correlation between GE and age in TD individuals (Hagmann et al., 2010); however, the investigators studied a much broader age range (18 months to 18 years old). It is possible that developmental trajectories of GE are nonlinear, given that other properties of structural development do not follow linear trajectories. For instance, white matter integrity exhibits nonlinear increases from childhood to adulthood, with steeper increases in childhood (Lebel et al., 2008). Therefore, future studies should study larger samples and broader age ranges of individuals with and without ASD to examine potential nonlinearities and group differences in developmental trajectories of information processing capacity in structural networks. Another possibility is that efficiency trajectories may differ in different parts of the brain: Dennis et al. (2013) found that GE increased with age in the left hemisphere, but decreased with age in the right hemisphere, in a sample of 439 adolescents and adults aged 12 to 30 years. Thus, future studies should aim to further elucidate region-or network-dependent development of information processing capacity.

Limitations
There are several limitations of the current study. First, our sample size was small due to the limited number of high quality datasets with both DTI and resting-state fMRI data from the ABIDE II database. We chose to analyze data from a single site for this study. While previous work has demonstrated the feasibility of multisite DTI studies based on high concordance of fractional anisotropy and mean diffusivity measurements across five different scanning sites (Fox et al., 2012), the reliability of measurements of resting-state FC across scanning sites is less clear. FC measurements can be affected by differences in scanner manufacturers and acquisition protocols (Shinohara et al., 2017;Yu et al., 2018). Dansereau et al. (2017) analyzed multisite resting-state FC data across 8 scanning sites, and reported small to moderate betweensite effects. Jovicich et al. (2016) studied difference in DMN FC in 5 healthy elderly participants who were scanned at 13 different sites with harmonized protocols, and found significant differences in temporal signal-to-noise ratio between sites. These differences were hypothesized to relate to differences in hardware and pulse sequences between sites. Nielsen et al. (2013) found that classification of ASD compared to controls based on resting-state FC was lower when multiple sites were included as opposed to a single site. Further, Easson et al. (2018) defined subtypes of ASD and controls based on FC using 5 sites from the ABIDE database, and found that it was necessary to regress out the effects of scan site prior to k-means classification, otherwise, the clusters differed significantly in scan site. This finding shows that prominent FC clusters defined in an unsupervised manner are driven by differences in FC across scanning sites. Thus, it will be crucial to further elucidate the effects of scanning site on resting-state FC to allow for multisite studies, and thus, larger sample sizes.

Conclusions
This study reveals relationships between brain signal variability and complexity and GE, age, and behavioural severity across ASD and TD participants, and illustrates the importance of taking a dimensional approach to studying brain function in ASD. By analyzing brain variability and complexity in relation to a set of predictor variables, a continuum of relationships between brain variables and predictor variables was observed; however, when treating ASD and TD groups categorically, significant group differences were not observed. Further, increased GE in structural networks and increased age were associated with higher MSSD and entropy, revealing important information about structure-function relationships in the brain and the developmental trajectories of variability and complexity.

Conflict of interest
The authors declare no competing financial interests.