Typical and Aberrant Functional Brain Flexibility: Lifespan Development and Aberrant Organization in Traumatic Brain Injury and Dyslexia

Intrinsic functional connectivity networks derived from different neuroimaging methods and connectivity estimators have revealed robust developmental trends linked to behavioural and cognitive maturation. The present study employed a dynamic functional connectivity approach to determine dominant intrinsic coupling modes in resting-state neuromagnetic data from 178 healthy participants aged 8–60 years. Results revealed significant developmental trends in three types of dominant intra- and inter-hemispheric neuronal population interactions (amplitude envelope, phase coupling, and phase-amplitude synchronization) involving frontal, temporal, and parieto-occipital regions. Multi-class support vector machines achieved 89% correct classification of participants according to their chronological age using dynamic functional connectivity indices. Moreover, systematic temporal variability in functional connectivity profiles, which was used to empirically derive a composite flexibility index, displayed an inverse U-shaped curve among healthy participants. Lower flexibility values were found among age-matched children with reading disability and adults who had suffered mild traumatic brain injury. The importance of these results for normal and abnormal brain development are discussed in light of the recently proposed role of cross-frequency interactions in the fine-grained coordination of neuronal population activity.


Preprocessing
The MEG data underwent artifact reduction using Matlab (The MathWorks, Inc., Natick, MA, USA) and Fieldtrip routines (Oostenveld et al., 2011). Independent component analysis (ICA) was used to separate cerebral from non-cerebral activity using the extended Infomax algorithm as implemented in EEGLAB (Delorme and Makeig, 2004). The data were whitened and reduced in dimensionality using principal component analysis with a threshold set to 95% of the total variance (Delorme andMakeig, 2004, Escudero et al., 2011). Kurtosis, Ré nyi entropy, and skewness values of each independent component were used to identify and remove ocular and cardiac artifacts. A given component was considered artifactual if, after normalization to zero mean and unit variance, more than 20% of its values exceeded 2 SDs from the mean (Escudero et al., 2011, Dimitriadis et al., 2013a, Antonakakis et al., 2013. Axial gradiometer recordings were transformed to planar gradiometer field approximations using the sincos method implemented in Fieldtrip (Oostenveld et al., 2011).

Intra and Inter-Frequency Coupling Estimators and statistical filtering
The center of the sliding window moved forward at 0.1-sec steps and the set of, were re-estimated leading to a series of 900 functional connectivity graphs ( TV FCGs) each representing a snapshot of the spatial profile of each connectivity index during the recording session. Surrogate data were generated by cutting at a single random time point of each time series and swapping the two resulting time courses (Dimitriadis et al. 2013). This approach was chosen over the shuffling-based surrogate analysis which significantly alters the spectrum of the time series and also impairs the non-stationarity of brain activity. Here we assessed determinism/stochastisity and linearity/nonlinearity of the surrogate time series using the Delay Vector Variance (DVV) method which relies on signal predictability in phase space to characterize the time series (Dimitriadis et al., 2017b). Preliminary analyses using DVV indicated that the surrogate data produced in this method displayed similar frequency composition and stationarity as compared to the recorded time-series.

Amplitude Envelope Correlation (AEC)
Surrogate data analyses ensured that significant observed AEC values could not have emerged from time series with zero amplitude envelope coupling. The final AEC data set consisted of the highest, significant correlations across all frequencies and frequency pairs (if more than two frequencies or pairs of frequencies were significant for a given pair of sensors we selected the one with the maximum AEC value). Finally, the significant, dominant AEC values for each pair of sensors and across sliding windows were integrated into a 4D graph ( TV AEC 4D array of size [900 x 248/275 x 248/275]). In this array significant AEC interactions were indicated by a value of 1, with zeros indicating non-significant AEC interactions.

Phase-to-amplitude Cross-Frequency Coupling (CFC, cross-frequency iPLV)
Let x(isensor, t), be the MEG activity recorder at the isensor-th site, and t=1, 2,.... T the successive time points. Given two band-limited signals x(isensor,t) and x(jsensor,t), cross-frequency coupling was estimated by allowing the phase of the lower frequency (LF) oscillation to modulate the amplitude of the higher frequency (HF) oscillation. The complex analytic representations of each signal zLF(t) and zHF(t) were derived via the Hilbert transform (HT[.] Next, the envelope of the higher-frequency oscillation AHF(t) was bandpass-filtered within the range of the LF oscillation and the resulting signal was submitted to an additional Hilbert transform to derive its phase dynamics component φ'(t): which expresses the modulation of the amplitude of the HF oscillation by the phase of the LF oscillation. Phase consistency between the two timeseries, corresponding to PAC strength, was measured by means of the imaginary portion of the Phase Locking Value (iPLV; Lachaux et al., 1999) defined as follows: The imaginary portion of PLV is considered to be less susceptible to volume conduction in assessing CFC interactions and was used in all subsequent analyses. While iPLV is not affected by volume conduction, it may be sensitive to changes in the angle between two signals, which do not necessarily imply a PLV change. In general, iPLV is only sensitive to non-zero-phase lags and is thus resistant to instantaneous self-interactions associated with volume conduction (Nolte et al., 2004).
For each pair of sensors, in each sliding window and for each subject, the maximum PAC score was determined for each combination of modulating frequency bands (fφ; from δ to γ) and modulated frequencies (fa; ranging from 1 to 90Hz in 1Hz steps) using the following equation: In this manner we determined if the prominent CFC interaction (indexed by the maximum PAC value), and specifically the frequency of the low-frequency phase, was associated with the frequency of maximum power. If the two frequencies were identical this would imply that the observed prominent CFC interaction was driven by the power of the dominant frequency and not its phase. In order to ensure that we did not include in further analyses CFC interactions of this type, we required a minimum frequency difference of 1Hz distance between the two frequencies (frequency of the low-frequency phase, frequency of the higher power), the one identified by equation 1 (maximum PAC value) and the one associated with the maximum spectral power, within the range of high amplitude. Surrogate data analyses further ensured that a very low p value associated with a given observed PAC indicated that it could not have emerged from time series with zero crossfrequency phase-amplitude coupling. To be included in surrogate analyses a given PAC value would have to meet an additional criterion, namely that the frequency of the modulated HF should be distinct from the frequency associated with maximum power. Here, we adopted a criterion of at least 1 Hz difference between the two frequencies.
The dominant PAC values for each pair of sensors and across sliding windows were integrated across frequency bins yielding 28 possible pairwise PAC estimates among the eight frequency bands. For each participant the resulting TV PAC profiles constituted two 4D arrays of size [28 (pairs of frequencies) x 900 (time windows) x 248/275 (sensors) x 248/275 (sensors)]. The first array contained the PAC values, and the second array the corresponding p-value. In a third array of size [2 x 900 x 248/275 x 248/275], dICMs were integer-coded (e.g., 1 for δ, 2 for θ,…., 11 for δα, 18 for θβ1).

Intra-Frequency Phase-to-Phase Coupling (same-frequency iPLV)
Intra-frequency phase coupling was estimated using the Hilbert phase transform implemented in the following formula adapted from the following equation: Based on surrogate data analyses an p value was assigned to each observed iPLV value, with very low p values suggesting that the observed iPLV could not have emerged from time series with zero phase coupling. Finally, the dominant iPLV mode for each sensor pair was determined based on the highest, statistically significant iPLV value. The significant, dominant iPLV values for each pair of sensors and across sliding windows were integrated into two 4D arrays ( TV iPLV) of size [8 x 900 x 248/275 x 248/275]). The first array contained the iPLV values and the second array the corresponding p-value. In a third array of size [2 x 900 x 248/275 x 248/275], the dICMs were integercoded (e.g., 1 for δ, 2 for θ,…., 11 for δα, 18 for θβ1).

Delay Symbolic Transfer Entropy (dSTE)
In principle, asymmetric dependences between coupled systems can be detected with measures of mutual information (Dimitriadis et al., 2010) taking into account the dynamics of information transfer. Transfer entropy (Shannon & Weaver, 1949), which is related to Granger causality (Granger, 1969), was introduced in order to distinguish the driving and responding elements within a network. Through proper conditioning of transition probabilities this quantity has been shown to be superior to the standard time-delayed mutual information, which fails to distinguish information, which is actually exchanged, from shared information due to common history and input signals. Various techniques have been proposed to estimate transfer entropy from observed data. Most techniques, however, make great demands on the data, require fine-tuning of parameters, and are highly sensitive to noise contributions, which limits the use of transfer entropy to field applications (Granger, 1969, Verdes, 2005. The algorithmic steps employed to transcribe the temporal dynamics of each pair of sensors into two distinct symbolic timeseries that share a common codebook (set of symbols) are described below. The size and content of the codebook is data-dependent and estimated causal relationships are to be inferred from a different pair of recorded signals each time. The unsupervised nature of the algorithm reduces computational burden (Dimitriadis et al., 2016a,c,d). Given signals A xt and B xt recorded from sensors A and B, respectively, time-delay vectors are first reconstructed from each time series. These vectors take the form of where m is the embedding dimension,  denotes the time lag and t=1,2,..,T runs over the time points. Then, the two individual sequences of timedelay vectors are collectively gathered in data matrices in the form of: Next, the two trajectories are brought to a common reconstructed state space by forming the overall data matrix Partitioning all the tabulated m-dimensional vectors into groups of homogenous patterns is the most direct way to summarize the temporal variations in the two time series and describe them with a common vocabulary. In our approach, a codebook of k code vectors is designed by applying the NG algorithm to the data matrix described by Eq. 2. The NG algorithm is an artificial neural network model, which converges efficiently to a small number stochastic gradient descent procedure with a soft-max adaptation rule that minimizes the average distortion error (Martinetz et al., 1993). In the encoding stage, each of the 2T vectors was assigned to the nearest code vector. By replacing the original vectors with the assigned code vectors, the two vectorial time series can be reconstructed with a measurable error. If we denote the reconstructed (i.e. decoded) version of the vectorial time series as x , the fidelity of the overall encoding procedure can be estimated via the index described in the following equation, which reflects the total distortion error divided by the total dispersion of the original vectors: , the better the encoding. This index becomes smaller as k increases, and reaches a plateau at a relatively small value of k . In the present study, we considered encoding to be acceptable if it was produced by the smallest k and if Distortion n was less than 5%. The NG algorithm was repeatedly applied at progressively higher k values to determine the optimal 0 k , which in turn identified the codebook to use in the subsequent symbolization scheme. At the vector-quantization stage, each A X and B X vector was assigned (according to the nearest-prototype rule) to the most similar among the derived code-vectors , which in mathematical notation reads as follows: In the derived symbolic time series, the temporal dynamics of a pair of neural subsystems are encoded as transitions among adaptively-defined (i.e. data-dependent) symbols. We adopted the Ragwitz criterion for optimizing the embedding dimension d and the embedding delay τ. Optimality of the embedding dimension refers to minimal prediction error for future samples of the time series. The Ragwitz criterion predicts the future of a signal based on estimates of the probability densities of future values of its nearest neighbors after embedding. The adopted method was based on the minimization of mean squared prediction error (Ragwitz and Kantz, 2002;Lindner et al., 2011).
dSTE as an index of Effective Connectivity. Given a pair of symbolic sequences A st and B st, the relative frequency of occurrence of symbols can be used to estimate joint and conditional probabilities, and to define symbolic transfer entropy (STE) as follows where the sum runs over all symbols and δ denotes a time step.
Effective connectivity is defined as "the influence one system exerts over another" (Granger, 1969;Ito et al., 2011). In the context of brain networks, effective interactions are directed from one sensor location to another. To account for the time-delay between neuromagnetic signals recorded from different sensors, Eq. 3 was adapted as follows: where d is the time delay between the driving and the driven time series.
A log of base 2 is used so that STEBA is expressed in bits. STEAB is defined in complete analogy. The directionality index dSTEAB=dSTEAB-dSTEBA quantifies the preferred direction of information flow taking up positive values for a unidirectional coupling, with time series A as the driver, and negative values for time series B driving A. For symmetric bidirectional couplings dSTE approximates zero. The formulation of transfer entropy with a time delay has been validated in a recent study, which presented a robust method for neuronal interaction delays (Wibral et al., 2013).
The significance of a particular type of interaction between two sensors was assessed through surrogate data (1000 time series of sensor B created through permutation of the elements of the symbolic time series B st). A one-sided p-value was computed that corresponded to the percentage of surrogate dSTEBA values that were higher than the observed dSTEBA. The final dSTE data set contained the strength, direction and delay of the significant and dominant pair of frequencies for each sensor. If more than two dSTE frequencies or frequency pairs were significant, the one with the maximum dSTE value was selected. For each participant the resulting TV dSTE profiles constituted two 4D arrays of size [28 (pairs of frequencies) x 900 (time windows) x 248/275 (sensors) x 248/275 (sensors)]. The first array contained the dSTE values and the second array the corresponding p-value. In a third array of size [2 x 900 x 248/275 x 248/275], the dominant dSTE interactions were integercoded (e.g., 1 for δ, 2 for θ,…., 11 for δα, 18 for θβ1).

Directed Phase Lag Index (dPLI)
Τhe dPLI quantifies the phase difference between two time series A and B which can be classified as follows: a) Signal A is consistently leading signal B in the phase domain (dPLI > 0.5), if the majority of phase differences fall within  Inconclusive causality between A and B in either direction, if the phase difference of the two signals averages π radians c) Signal B is consistently leading signal A in the phase domain (dPLI < 0.5), if the majority of phase differences fall within Significant, causal phase relationships for every pair of sensors were determined via surrogate data analyses to ensure that a very low p value indicated that the observed dPLI could not have emerged from time series lacking directed phase interactions.

Topological Filtering based on Orthogonal Minimal Spanning Tress (OMSTs)
The general concept of the OMST method can be described as follows: The initial (1 st ) MST connects all the V sensors through V-1 edges. Then, the V-1 connections of the 1 st MST were substituted with zeros and a 2 nd MST was estimated that connected all of the V sensors with minimal total distance, satisfying the constraint that it was orthogonal-i.e. it shared no common edges-with the 1 st MST. Next, the V-1 connections of the 2 nd MST were substituted with zeros and a 3 rd MST was estimated that connected the sensors with the minimal total weight, subject to the constraint that it was orthogonal to the previous two constructed MST's (1 st and 2 nd ). In general, an m-MST is orthogonal to all the previous (m-1) MST's, having exactly m(N-1) edges. The OMST method was applied as follows: For each added N-1 edges correspond to a single OMST across multiple round OMTS, the objective function of Global Cost Efficiency (GCE) = GE-Cost was estimated, where Cost denotes the ratio of the total weight of the selected edges, over multiple iterations of OMST, divided by the total strength of the original full-weighted graph. The values of this formula range within the limits of an economical small-world network for healthy control participants. The network which is considered as functionally optimal is the one associated with the maximum value of the following quality formula:

= −
Topological filtering was applied next in order to reduce the spatial density of the sensor network featuring dICM weights. The curve maximum in the following Figure illustrates the optimization of the OMST algorithm for a healthy participant from the first temporal segment.

Supplementary metrics of brain activity & connectivity
For comparison purposes two well-established measures of brain activity and sensor interdependence were also computed, namely relative power spectrum and coherence. In contrast to dICM indices they provide static representations of the strength of oscillatory activity and withinfrequency phase coupling, respectively.

Imaginary Coherence between sensors (ImCOH)
Coherence is a widely used measure of linear dependence and phase consistency between two time series. Let si(f) and sj(f) represent the complex Fourier transforms of two time series xi(t) and xj(t). For each of N epochs of the segmented time series, the cross-spectral density function Sij(f) of i and j is defined as: = ( ) * ( ) and the single epoch power spectrum (auto-spectral density) as: where * indicates the complex conjugate. The coherency across n epochs is defined as the crossspectrum normalized by the auto spectral density functions xi(t) and xj(t): Coherence is defined as the absolute value of coherency, i.e. the normalized amplitude of the complex cross-spectrum value: In the present study we adopted the imaginary portion of coherence (Nolte et al., 2004) corresponding to the perpendicular presented to the real axis values of a complex number, positively increasing in magnitude to the left. In this manner, when two time series are completely simultaneous, i.e. they have zero-lag oscillatory activity, the value of ImCoh is zero: Thus ImCoh reflects only non-zero time-lag correlations which are presumably caused by interactions between distinct neuronal populations rather than common interference sources.

Multiscale Entropy (MSE)
Multiscale entropy (Costa et al., 2002(Costa et al., , 2005) is a novel method to estimate the complexity of a finite length time series. For the purposes of the present study an MSE algorithm was developed for neuromagnetic time series using sample entropy (Richman and Moorman, 2000). Sample Entropy reflects the probability that sequences that match each other on the first two data points will also match on the next point (Costa et al., 2002 In contrast to traditional entropy indices, MSE provides a more direct and accurate measure of complexity by overcoming key limitations of the former, such as their susceptibility to bias introduced by surrogate data comparisons. The MSE method involves two steps: 1. A coarse-grained procedure is first applied to the time series. For a given time series, multiple coarse-grained time series are constructed by averaging the data points within non-overlapping windows of increasing length τ (see Figure S1). Each element of the coarse-grained time series, ) ( j y  , is calculated according to the equation: where τ represents the scale factor and 1≤j≤N/τ. The length of each coarse-grained time series is N/τ. For scale value = 1, the coarse-grained time series simply corresponds to the original time series. 2. Sample Entropy is then plotted as a function of scale values. As a "regularity statistic" SampEn searches for patterns within a time series and quantifies its degree of predictability or regularity (see Figure S2). Here, we estimated MSE separately for each participant, frequency band, and lobar section of the sensor array.  , to determine the ratio between the total number of 2-component template matches and the total number of 3-component template matches. Sample Entropy is the natural logarithm of this ratio and reflects the probability that sequences that match each other for the first two data points will also match for the next point (Costa et al., 2002). Drucker et al. (1997) extended the Support Vector Machine method to SVR in order to make continuous real-valued predictions. SVR retains some of the main features of SVM classification, although in the former method misclassified cases are penalized, whereas in SVR a penalty is served only for cases associated with extreme distances from the regression line in high-dimensional space (Dosenbach et al., 2010). Epsilon-insensitive SVR defines a tube of width ε, which is user defined, around the regression line in high-dimensional space. Any points within this tube carry no loss. In essence, SVR performs linear regression in high-dimensional space using epsilon-insensitive loss. The C parameter in SVR controls the trade-off between how strongly points beyond the epsiloninsensitive tube are penalized and the flatness of the regression line (larger values of C allow the regression line to be less flat; Dosenbach et al., 2010). SVR predictions described in the present analyses used epsilon-insensitive SVRs carried out in the Spider Machine Learning environment (Weston et al., 2005), as well as custom scripts run in MATLAB (R2010a; MathWorks, Natick, MA, USA). The parameters C and ε were tuned using a holdout subset of the respective dataset. An important feature of SVR is that it allows determination of the most important model features in predicting participant age, according to their corresponding weights.

Prediction of participant age group using Multi-class Support Vector Machines (MC-SVM)
Multi-class Support Vector Machines represent an extension of Support Vector Machines which are appropriate for binary classification (k = 2). Multi-Class pattern recognition (k>2 classes) are usually solved using a voting scheme based on combining many binary classification functions simultaneously. Here, we used a linear MC-SVM which involved training on all six age groups simultaneously, in order to build a unique decision function for each group (Crammer and Singer, 2001 and implemented in (http://www.cs.cornell.edu/people/tj/svm_light/svm_multiclass.html). Briefly, this method constructs

Cross Validation of SVRs and MC-SVMs
Cross validation of the age prediction and classification, respectively, was performed using three complementary schemes with each tool. In the first two, data from one MEG system were used to train the SVR and data from the second MEG system served as the testing set. The third scheme was applied to data aggregated across both MEG datasets. The prediction/classification capacity of each scheme was assessed using two alternative procedures: 5-fold cross-validation and leave-one-outcross-validation (LOOCV; a method associated with the most unbiased test error estimates; Hastie et al., 2001). On each LOOCV round, one case is left out and the remaining cases are used as the training set. During each SVR round (LOOCV) or fold, the correlation between each feature and the dependent variable (age) was computed, and the features that had the highest absolute correlation values were selected for inclusion in the final prediction model.

Identifying age-related neuromagnetic features
S. Table 1 presents results of various cross-validation schemes in the context of Support Vector Regression in predicting participant age based on the combination of 22 optimal features characterizing the short-term persistence of dICMs in 12 subnetworks (as defined by Mean Subgraph Strength) and 10 subnetworks (as defined by Fractional Occupancy). S. Table 1 illustrates the 22 optimal features that classify correct each participant to each age group based on Mean Subgraph Strength (MSS) and Fractional occupancy (FO). This model was able to predict participant age with a reasonable degree of accuracy according to the following regression equation: Predicted Age = 1.01 x Actual Age -0.78 (R 2 =0.893, p < 2× 10 −9 ; Fig. S3). Corresponding results using Multi-Class Support Vector Machines to predict the correct age group of each participant are shown in S. Table 2. The most robust age prediction was achieved using the 5-fold cross validation method with data aggregated across both systems (as indicated by the capacity of the Support Regression model to account for 90% of individual variability in participant age and by 89% correct classification of participants to the appropriate age group). Moreover, the comparability among the two data sets (Magnes-248, CTF-275) is evident by the high rates of successful age and age-group prediction when the Support Vectors were trained with the same set of dICM features from one system and tested with corresponding features from the other system.  Figure S3. (Upper panel) Least squares regression of predicted over actual participant age in years based on Mean Subgraph Strength and Fractional Occupancy in 12 and 10 lobar subnetworks, respectively (listed in S. Table 2). (Lower panel) The distribution of model regression residuals as a function of participant age.

Age-Related Differences based on Dominant Intrinsic Coupling Modes (dICMs)
The spatial layout of dICMs which contributed significantly to the prediction of participant age is shown in the upper panels of Figures 1 and S4-14. The relative contribution of each of the 22 dICM features in the significant classification of participants to their correct age group through SVR is shown in the upper panel of S16. The prediction accuracy of participant age for each of the 22 features is displayed in the lower panel of S16. Figure S4. Dominant inter-and intra-hemispheric Phase-Amplitude Coupling (PAC) between frontal and temporal sensors in the θ and γ2 bands, respectively. A) Topographical layout of statistically significant sensor pairs for the six age groups. B) Mean Subgraph Strength and C) Fractional Occupancy (FO) derived from PAC across the six age groups. Significant differences between successive age groups are marked by brackets (p < .0001). Figure S5. Dominant inter-and intra-hemispheric coupling between frontal and parietal sensors indexed by delay Symbolic Transfer Entropy (dSTE) in the θ and α2 bands, respectively. A) Topographical layout of statistically significant sensor pairs for the six age groups. B) Mean Subgraph Strength and C) Fractional Occupancy (FO) derived from dSTE across the six age groups. Significant differences between successive age groups are marked by brackets (p < .0001). Figure S6. Dominant inter-hemispheric coupling between parietal and occipital sensors indexed by imaginary Phase Locking (iPLV) in the α1 band. A) Topographical layout of statistically significant sensor pairs for the six age groups. B) Mean Subgraph Strength and C) Fractional Occupancy (FO) derived from iPLV across the six age groups. Significant differences between successive age groups are marked by brackets (p < .0001).          Figure S16 displays average RP values at each frequency band across the six age groups. The prediction of age though a linear combination of 12 RP features was significant, albeit considerably poorer than the prediction using dICM features (R 2 =0.430, p < 3.1 x 10-4, see Figure S17). The most consistent age-related features corresponded to data from 12 sensors located over Frontal and Parieto-Occipital regions, the individual contribution of which to the prediction of participant age is shown in Figure S18.

Imaginary Coherence
Figure S19 displays Mean Subgraph Strength values within and between lobar regions in each frequency band. Figure S20 illustrates the significant prediction of age via linear SVR (R 2 =0.53, p < 1.9 x 10 -6 ) based on 10 ImCOH features. The relative contribution of individual features significantly contributing to the SVR classification model is shown in Figure S21.    Figure S22 demonstrates the lobe-averaged MSE profile for each age group and frequency band. Figure S23 shows the 4 features that significantly correlated with age and, combined, resulted in the prediction of age via linear SVR (R 2 =0.71, p < 2.9 x 10 -6 ) demonstrated in Figure S24. The relative contribution of individual features significantly contributing to the SVR classification model is shown in Figure S25.

Prediction of participant age using RP, ImCoh and MSE features combined
By comparison the capacity of supplementary metrics (RP, ImCoh, and MSE combined) to predict age via regression analysis demonstrated similar results (R 2 =0.812, p < 2.2× 10 −7 ; Fig. S26 & S. Table 1).

MEG system equivalence
Results from the cross-validation schemes using dICM-based features derived from one MEG system to train the Support Vector Regressor and Multi-Class Support Vector Machine, revealed performance that was comparable, albeit somewhat lower, than the performance obtained when using mixed data sets for training and testing. The prediction of age was similar across systems using the same set of features each time: R 2 =0.876, p < 1.8× 10 −7 when using data from the Magnes-248 system as training-set and R 2 =0.869, p < 3.5× 10 −7 when using data from the CTF-275 as training set (S. Table  1).
Moreover, classification of participants to the correct age group averaged % 77 . 5 65 . 86  when using data from the Magnes-248 system as the training set, and % 03 . 6 09 . 86  when using data from the CTF-275 system as the training set (S. Table 2). The similarity of subnetwork topography and age trends across systems is illustrated in Figure S27 using the dCIM between Frontal lobes as example. Estimates of FI were also comparable across systems as demonstrated by the close similarity of the maturation curves in Figure S28. Figure S27. dICM reflecting interhemispheric interactions between frontal sensors as indexed by amplitude envelop correlation in the δ band (AEC). A) Topographical layout of the statistically significant sensor pairs for 4 age groups (overlapping age range between the two MEG systems). B) Mean subgraph strength and C) Fractional Occupancy (FO) across the four age groups. The upper half of the figure presents data recorded on the Magnes-248 system and the lower half displays data recorded on the CTF-275 system. Figure S28. Functional brain maturation curves based on the Flexibility Index (FI) computed for data obtained on the Magnes-248 (A; n=81 aged 6-59 years) and CTF-275 MEG systems (B; n=97 aged 18-60 years). Chronological age is shown on the x axis. The best-fitting curves for the data for each system are shown by the blue lines.

Reliability of the Flexibility Index
One-week test-retest assessment of FI values indicated excellent reliability (r = 0.94) with average test-retest differences = 0.027±0.015, range 0.012-0.029 (see Figure S29). Figure S29. Test-retest Flexibility Index (FI) values as a function of participant age (n=10).