Flexible nonlinear modeling reveals age-related differences in resting-state functional brain connectivity in autistic males from childhood to mid-adulthood

Abstract

Background

Divergent age-related functional brain connectivity in autism spectrum disorder (ASD) has been observed using resting-state fMRI, although the specific findings are inconsistent across studies. Common statistical regression approaches that fit identical models across functional brain networks may contribute to these inconsistencies. Relationships among functional networks have been reported to follow unique nonlinear developmental trajectories, suggesting the need for flexible modeling. Here we apply generalized additive models (GAMs) to flexibly adapt to distinct network trajectories and simultaneously describe divergent age-related changes from childhood into mid-adulthood in ASD.

Methods

1107 males, aged 5–40, from the ABIDE I & II cross-sectional datasets were analyzed. Functional connectivity was extracted using a network-based template. Connectivity values were harmonized using COMBAT-GAM. Connectivity-age relationships were assessed with thin-plate spline GAMs. Post-hoc analyses defined the age-ranges of divergent aging in ASD.

Results

Typically developing (TD) and ASD groups shared 15 brain connections that significantly changed with age (FDR-corrected p < 0.05). Network connectivity exhibited diverse nonlinear age-related trajectories across the functional connectome. Comparing ASD and TD groups, default mode to central executive between-network connectivity followed similar nonlinear paths with no group differences. Contrarily, the ASD group had chronic hypoconnectivity throughout default mode-ventral attentional (salience) and default mode-somatomotor aging trajectories. Within-network somatomotor connectivity was similar between groups in childhood but diverged in adolescence with the ASD group showing decreased within-network connectivity. Network connectivity between the somatomotor network and various other functional networks had fully disrupted age-related pathways in ASD compared to TD, displaying significantly different model curvatures and fits.

Limitations

The present analysis includes only male participants and has a restricted age range, limiting analysis of early development and later life aging, years 40 and beyond. Additionally, our analysis is limited to large-scale network cortical functional parcellation. To parse more specificity of brain region connectivity, a fine-grained functional parcellation including subcortical areas may be warranted.

Conclusion

Flexible non-linear modeling minimizes statistical assumptions and allows diagnosis-related brain connections to follow independent data-driven age-related pathways. Using GAMs, we describe complex age-related pathways throughout the human connectome and observe distinct periods of divergence in autism.

Background

Divergent brain connectivity in autism spectrum disorder (ASD) has been extensively studied in humans using resting-state functional magnetic resonance imaging (rs-fMRI) [1,2,3,4]. Resting-state fMRI can be used to describe functional brain connectivity, or correlations of between-region brain activity, across brain networks that are associated with various cognitive functions. Previous work has reported widespread patterns of hypo- and hyper-connectivity both within and between these functional networks in ASD [1,2,3,4]. Moreover, these functional differences have been shown to relate to symptom severity and symptom profiles in autism. For instance, the degree of default mode network (DMN) hypoconnectivity correlates to social symptom severity in ASD [3]. Likewise, increased local connectivity of the dorsal posterior cingulate cortex (PCC) has been correlated with cognitive inflexibility and behavioral rigidity in ASD [1]. Despite these findings, functional connectivity results and their relationships to cognition have varied, even among modern large-scale mega-analyses, calling into question the reliability and reproducibility of functional connectivity differences in ASD [5,6,7].

Even more elusive has been understanding the effects of age on brain connectivity trajectories in ASD. Some studies have reported specific default mode network (DMN) within-network hyperconnectivity at a young age which resolved in adolescence to adulthood [8–9]. Another study found childhood DMN, frontoparietal, and salience between-network hyperconnectivity which either persisted into adulthood or decreased to hypoconnectivity in adulthood [10]. These connectivity changes may have been influenced by anatomical distance both within and between networks, with long-range connectivity specifically subject to decreased connectivity with age in ASD [10–11]. However, the specific findings are inconsistent across studies, possibly due to small sample sizes. Indeed, more recent mega-analyses have reported no significant impact of age on divergent connectivity in ASD [5–6].

These inconsistent findings may be due, in part, to the use of varying statistical approaches to relate age and brain connectivity in ASD. Many previous reports have binned participants into age groupings, such as child (< 11), adolescent (11 < 18), and adult (> 18), when analyzing brain connectivity changes with age in ASD [8, 10, 11]. However, prior work indicates that developmental connectivity changes occur within these age ranges [12,13,14]. As such, it may be critical to maintain age as a continuous variable when examining age impacts on brain connectivity.

Beyond age binning, many previous reports have employed linear models to describe brain age-related changes in ASD [5, 6, 15]. However, substantial evidence suggests that some age-related network functional connectivity changes may be nonlinear during typical development (TD) [13, 16,17,18,19,20,21]. Some prior studies that examined age-related functional connectivity differences in ASD with non-linear analyses limited their modeling to quadratic regression approaches [9, 25, 26]. Polynomial regression studies have reported that TD connectivity-age relationships can vary in polynomial degree based on functional network [19,20,21]. This diversity in polynomial fit across networks suggests the need for flexible modeling to best describe connectivity-age relationships.

The present study employs generalized additive models (GAMs) to flexibly fit linear and nonlinear typically developing connectivity-age relationships [13, 16, 17]. GAMs are data-driven predictor functions whose curvatures are derived during data fitting [22,23,24]. As such, GAMs minimize assumptions about data curvature prior to modeling, allowing non-linearities to vary across multiple comparisons. Further, GAMs can resolve data inflection points, allowing identification of developmental periods where participant groups diverge. Here, we revisit the ABIDE datasets considering non-linear, flexible modeling to better describe the relationship between functional connectivity and age in ASD.

Methods

Participants

All 2,226 participants with rs-fMRI imaging included in the ABIDE I & II datasets were initially considered for inclusion in the present analysis. Cross-sectional data from our prior University of Utah longitudinal study (R01MH080826) are included in this dataset. For site-specific details on both inclusion/exclusion criteria and more detailed characterization of clinical and behavioral phenotypes, see the ABIDE I & II releases [27, 28].

After MRI preprocessing, scrubbing, visual data inspection, and exclusion of sites with fewer than 10 participants (see MRI preprocessing and analysis), 1111 male participants from 25 unique sites remained. Four additional male TD individuals were excluded as statistical age outliers (z > 4) to maintain an age-matched dataset. 291 female participants also remained eligible for analysis following data quality control. However, females were not included in the analysis due to disproportionate sex distribution and known sex differences in ASD [29]. Further details on data quality control and participant inclusion are described elsewhere [7].

MRI preprocessing and analysis

Structural MPRAGE data from the ABIDE dataset were processed using FreeSurfer (v6.0.0), which is documented and freely available for download online (http://surfer.nmr.mgh.harvard.edu/). Briefly, this processing includes motion correction and averaging [30] of volumetric T1 weighted images, removal of non-brain tissue using a hybrid watershed/surface deformation procedure [31], automated Talairach transformation, segmentation of the subcortical white matter and deep gray matter volumetric structures [32–33], intensity normalization [34], tessellation of the gray matter white matter boundary, automated topology correction [35–36], and surface deformation following intensity gradients to optimally place the gray/white and gray/cerebrospinal fluid borders at the location where the greatest shift in intensity defines the transition to the other tissue class [37,38,39].

Preprocessing of the ABIDE fMRI blood-oxygen-level-dependent (BOLD) data was performed in MATLAB (MathWorks, Natick, MA, USA) using SPM12 (Wellcome Department of Imaging Neuroscience, London, UK). All images were corrected for motion using a realign and unwarp procedure. Each participant’s BOLD images were coregistered to their individual MPRAGE anatomic image sequence. Phase-shifted soft tissue correction (PSTCor) [40] was used to regress participant motion parameters, eroded white matter, eroded cerebral spinal fluid, and soft tissues of the face and calvarium. PSTCor, an alternative to global signal regression (GSR), does not explicitly regress global fMRI signal but rather related signals, as GSR has been shown to artifactually inflate anti-correlations between functional brain networks [40, 41]. Eroded masks were obtained by removing all voxels from white matter and CSF masks that were adjacent to a voxel not in the mask. Volume censoring (scrubbing) was performed with removal of volumes before and after mean framewise displacement head motion greater than 0.3 mm [42]. Only participants with ≥ 50% volumes remaining after scrubbing were considered for further analysis. Additionally, quality control was completed by J.K. that consisted of visual inspection of each participant’s neuroimaging data for successful completion of the preprocessing pipeline and image quality. Individual research sites with less than 10 participants remaining after quality control were removed from further analysis.

Following preprocessing, rs-fMRI images were parcellated using the Yeo et al. 17 × 17 functional network schema [43]. Functional network connectivity was calculated using the Pearson correlation coefficient between pairs of time series from the parcellated region. This yielded a 17 × 17 connectivity matrix comprising 136 unique inter-network correlations.

Data harmonization

Connectivity data was then harmonized across sites using ComBat-GAM [44]. ComBat is an empirical Bayesian harmonization technique aimed to mitigate additive and multiplicative site-effects of multicenter neuroimaging data while maintaining biological variability [45, 46]. ComBat traditionally relies on a multiplicative linear regression framework and further assumes a normal distribution of site effects. However, the relationship between brain connectivity and age is presently hypothesized to be non-linear and non-parametric. As such, we implement ComBat-GAM, a generalized additive model ComBat framework that allows preservation of nonlinear covariate effects [44].

Generalized additive models

GAMs, unlike parametric regressions, do not contain a priori specification of predictor relationships [22,23,24]. Rather, GAM prediction functions are derived during model estimation [22,23,24]. This process eliminates regression model refinement and predictor selection, which is burdensome and error-prone when analyzing hundreds to thousands of unique models during large connectivity analysis [22].

Thin-plate splines (TPS) are the suggested optimal basis function to enable independent smooth functions of multiple predictors in GAMs [24, 45]. TPS use a penalized least squares method to mitigate over-smoothing while maintaining model flexibility [24, 47]. The final GAM curve is the weighted summation of a series of TPS fits in the data [22,23,24]. In doing so, GAM relationships can span a spectrum of linearity and non-linearity. Further, GAMs can resolve sharp data inflection points, which are often obscured in polynomial regression.

The GAMs presently used follows the mathematical framework:

$$y{\rm{ }} = {\rm{ }}{\beta _o}{\rm{ }} + {\rm{ }}\sum \left( {{\beta _i}{\rm{ }}*{\rm{ }}{x_i}} \right){\rm{ }} + {\rm{ }}f\left( {{x_1},{\rm{ }}{x_2}} \right)$$

Where y is the predicted variable, β₀ is the intercept, βi is the linear coefficient for the linear predictor variable xi, and f(x) is the smooth function of predictor variables. The smooth function was defined by the thin-plate spline:

$$\phi \left( r \right){\rm{ }} = {\rm{ }}{r^2}{\rm{ }}*{\rm{ }}log\left( r \right)$$

Where r is the Euclidean distance between a data point and a given knot, following the penalized least squares optimization:

$$\sum {\left( {{y_i}{\rm{ }} - {\rm{ }}f\left( {{x_i}} \right)} \right)^2}{\rm{ }} + {\rm{ }}\lambda {\rm{ }}*{\rm{ }}\smallint \smallint {\rm{ }}{\left( {{\nabla ^2}f\left( {x,y} \right)} \right)^2}{{ dx\,dy}}$$

Where λ is the smoothing parameter and ▽² is the Laplacian operator. Taken together, our general additive model follows the equation:

$$\eqalign{ Brain\;Connection & & \cr & \sim {\beta _o}\; + {f_1}(Age) + {f_2}(Age\;Group) \cr & + {\beta _{Group}}{}^*Group \cr & + {\beta _{Head\;Motion}}{}^*Head\;Motion\; + \in \cr} $$

Where ​f₁(Age) is a TPS of TD age, f₂​(Age, Group) is a contrasted TPS of ASD age, β_Group*Group is a predictor of Group, and β_{Head Motion}* Head Motion is treated as a nuisance regressor.

• A. Statistical analysis and visualization.

All statistical analyses and visualizations were conducted in R [48] (Version 4.4.1). Data were visualized using the ggplot2 R package [49]. Following harmonization, the relationship between each brain connection and participant age was analyzed using a TPS GAM with a maximum of 10 knots using the R mcgv package [50] (version 1.9.1).

Model significance was defined by age-related variables and was corrected for multiple comparisons using a Benjamini-Hochberg correction. Significance was defined as a connectivity-age relationship with corrected p < 0.05. For models in which the TD and ASD functional connectivity development trajectories both significantly changed with age, post-hoc pairwise analysis of the predicted models was used to define age ranges of group differences in aging using the tidy-gam package [51] (version 1.0.7). Specifically, pairwise significant differences were defined as where the 95% confidence interval of the difference between predicted curves was bound above or below 0.

Model evaluation

All significant models were evaluated to ensure the statistical assumptions of our GAM, namely residual homoscedasticity and residual normal distributions, were not violated. Residual homoscedasticity was assessed by plotting model residuals against predicted values. Residual distribution was visualized with a Q-Q plot and quantified with skewness, kurtosis, and Shapiro-Wilks normality tests.

Model reliability

Model reliability was evaluated with a stratified 5-fold split-sample validation. The mean correlation and standard deviation between the split-sample curve fits and the entire model fits at each year (5–40) were calculated to quantify model similarity. TD and ASD curves were independently analyzed and compared with a two-sample t-test across curve correlations to determine if reliability was group dependent.

Comparison to prior methods

Previous studies have commonly implemented both linear models and age-binning approaches to describe the relationships between functional connectivity, age, and ASD [5, 6, 8, 10, 11, 15]. To compare our model to these previous approaches, we implement both a linear model and age-binned ANCOVA on our harmonized dataset. Ages were binned as < 11, 11–18, and > 18, mirroring previous studies [8, 10, 11]. The linear model followed the equation:

$$\:Brain\:Connection\hspace{0.17em}\sim\hspace{0.17em}Age*Group\hspace{0.17em}+\hspace{0.17em}Head\:Motion$$

while the ANCOVA followed the equation:

$$\:Brain\:Connection\hspace{0.17em}\sim\hspace{0.17em}Binned\_Age*Group\hspace{0.17em}+\hspace{0.17em}Head\:Motion$$

where mean motion was treated as a nuisance regressor. P-values for the main effects of group, Age, and their interaction were extracted and Benjamini-Hochberg corrected for 136 comparisons. All statistical analyses were conducted in R [46] (Version 4.4.1).

Replication with Yeo 7 parcellation

Brain parcellation schemas have recently been shown to impact relationships between functional network connectivity and age [52, 53]. To this end, we presently attempted to replicate our findings using the Yeo et al. 7 × 7 functional network parcellation [43]. Importantly, the Yeo 7 × 7 schema does not allow for within-network analyses, limiting reproducibility comparisons to between-network analyses. Replication data processing, harmonization, modeling, and analyses follow previously described methods.

Results

Participant demographics

In the ABIDE I data release, 1112 participants were analyzed and 419 were removed (220 due to motion, 199 due to other quality issues). In the ABIDE II data release, 1114 participants were analyzed and 405 were removed (217 due to motion, 188 due to other quality issues). After excluding the remaining females and four age outliers (z > 4), a total of 1107 age-matched males ranging from 5 to 40 years old (N[Mean, STD], TD = 610[15.18, 6.38], ASD = 497[14.90, 6.01]) were included in the final analysis (Fig. 1). Of the final 1107 participants included, 15 ± 12% of volumes were censored for TD participants and 18 ± 13% of volumes were censored for ASD participants across sites (t (1105) = 4.16, p < 0.001, 95% CI [0.02–0.05]). There was no significant difference in mean age between the groups. Age distribution in both groups followed a positive-skewed Gaussian distribution.

Fig. 1

Participant age range. Bars are a histogram of participant age binned into integer age ranges. Curves are density plots of the age distribution by group

Diverse trajectories of functional brain development

Of the 136 analyzed brain functional connections, 55 changed significantly with age in the TD group (p-FDR < 0.05) (Fig. 2). In the ASD group, 19 functional connections changed significantly with age (p-FDR < 0.05), of which, 15 overlapped with TD. Nine of these 15 connections displayed significantly different connectivity-aging trajectories in ASD compared to TD. Curvature of these connectivity-age trajectories were diverse both between and within networks, highlighting model flexibility enabling analysis of differential age-related pathways throughout brain networks.

Fig. 2

Trajectories of functional brain development. Each individual curve displays connectivity strength (y-axis) with increasing age (x-axis, range 5–40). Connections that varied significantly with age after multiple comparison corrections are shaded in gray. *Common significant connections shared by both ASD and TD. Significant between group differences highlighted in gold

Fully disrupted linear age-related pathways in ASD

Post-hoc analyses of overlapping significant trajectories revealed six disrupted age-related pathways in ASD (Fig. 3). These six overlapping FDR-corrected between-network age-related trajectories displayed significantly different model fits and shapes between the ASD and TD groups. In all six models, non-linear age-related changes are seen in the TD group with increases in childhood and adolescence followed by stabilization or decreases in connectivity in adulthood (p-FDR < 0.05). Contrarily, the ASD group connectivity increases linearly with age throughout the analyzed age range (p-FDR < 0.05). These disrupted linear pathways are seen between the somatomotor network (SMN) and the dorsal attention (DAN), limbic, and medial superior parietal central executive (MSP-CEN) networks. Likewise, both the anterior ventral attention (aVAN) and MSP-central executive networks display these disrupted models in their connectivity to posterior dorsal attention (pDAN) network, different from the congruent MSP-central executive to VAN age-related pathways as seen in Fig. 4. Lastly, limbic-aVAN connectivity also displayed a disrupted linear path.

Fig. 3

Fully disrupted linear age-related pathways in ASD. Between-network connections displayed fully disrupted age-related pathways in ASD, demonstrated by entirely different model fits compared to TD. TD between-network SMN development increased in childhood and stabilized or decreased in young adulthood. Contrarily, ASD SMN connectivity was best modeled by an increasing linear model, displaying group-by-age interactions in ASD. Similarly, both the anterior VAN and MSP-CEN had disrupted connectivity-age-related pathways with the posterior DAN. A similar trend was seen in limbic-VAN connectivity. All results are FDR-corrected. Shaded areas denote regions of pairwise differences between groups. Acronyms: somatomotor network (SMN), dorsal attention network (DAN), medial temporal limbic (MT-limbic), medial superior parietal (MSP), ventral attention network (VAN)

Fig. 4

Default mode between-network chronic underconnectivity in ASD. The lateral temporal DMN subcomponent showed decreased between-network connectivity to both the anterior VAN and somatomotor networks throughout the analyzed period. The general trajectories of the groups followed a similar model curvature. Shaded areas denote regions of pairwise differences between groups. Acronyms: somatomotor network (SMN), ventral attention network (VAN), default mode network (DMN)

DMN between network chronic underconnectivity in ASD

Between network connectivity of the lateral-temporal DMN and the anterior VAN varied significantly with age in both the ASD and TD group (TD: p-FDR < 0.001; ASD: p-FDR = 0.005) (Fig. 5). Both groups followed an inverse parabolic connectivity-age path. The ASD group had significantly decreased DMN-VAN connectivity that persisted from age 8–26 (pairwise 95% C.I. > 0) and continued non-significantly throughout the analyzed age range. A similar chronic underconnectivity was seen between the lateral-temporal DMN and the ventral SMN (TD: p-FDR < 0.001; ASD: p-FDR = 0.004), with significant hypoconnectivity ranging from ages 5–31 (pairwise 95% C.I. > 0).

Fig. 5

Within somatomotor network divergent age-related trajectories in ASD. Both ASD and TD within-somatomotor connectivity followed an increasing sinusoidal aging pathway through childhood to mid adolescence. Beginning in late adolescence, the ASD group diverged to sustained hypoconnectivity as the curves horizontally stabilized. Shaded area denotes region of pairwise differences between groups. Acronyms: somatomotor network (SMN)

Within somatomotor network divergent age-related differences in ASD

Within-network somatomotor connectivity varied significantly with age in both ASD and TD groups (Fig. 6). Both ASD and TD SMN connectivity followed a sinusoidal age-related pathway (TD: p-FDR < 0.001; ASD: p-FDR < 0.001). Connectivity paths between the two groups were similar through childhood but significantly diverged in late adolescence with ASD showing decreased within-network connectivity from ages 18–27. This decreased connectivity trend continued from between ages 27–40, though not statistically significant.

Fig. 6

Congruent age-related group trajectories between ASD and TD. CEN between-network connectivity both increased and decreased depending on separate sub-components of the DMN. MSP-CEN to VAN connectivity increased at a young age and horizontally stabilized beginning in early adulthood. A similar trajectory was seen in CEN-ATTN and SMNa-ATTN connectivity. Acronyms: medial superior parietal central executive network (MSP-CEN), default mode network (DMN), ventral attention network (VAN), somatomotor association network (SMNa)

Congruent aging between ASD and TD

Six between-network connections followed congruent age-related paths with no significant group differences between ASD and TD (Fig. 4). In both groups, both connectivity between the ventral and dorsal DMN and MSP-CEN changed significantly with age. Connectivity between the MSP-CEN and ventral DMN increased linearly with age in both TD and ASD groups (TD: p-FDR = 0.001; ASD: p-FDR = 0.013). In contrast, connectivity between MSP-CEN and dorsal DMN decreased in a curvilinear manner throughout the analyzed period (TD: p-FDR = 0.005, ASD: p-FDR = 0.004).

Both subcomponents of the VAN varied significantly with age with the MSP-CEN in both the ASD and TD groups. VAN-CEN connectivity changes with age followed similar pathways for both posterior and anterior VAN networks, with an increase in childhood that leveled horizontally in early adulthood (CEN-pVAN: TD: p-FDR < 0.001; ASD: p-FDR = 0.004) (CEN-aVAN: TD: p-FDR < 0.001; ASD: p-FDR = 0.003). Lastly, both the MSP-CEN and pVAN networks curvilinearly increased significantly with the somatomotor association (SMNa) network (CEN-SMNa: TD: p-FDR < 0.001; ASD: p-FDR = 0.002) (pVAN-SMNa: TD: p-FDR < 0.001; ASD: p-FDR = 0.003).

Model evaluation

Given the positive-gaussian age skew of our datasets, all models were evaluated for statistical validity. Model evaluations showed that all 15 reported models did not violate the assumptions of residual homoscedasticity. Further, residual distributions all had a Shapiro-Wilks value W > 0.992, a kurtosis between 2.922 and 3.305, and a skewness between 0.022 and 0.351. Further, residual distribution did not differ between ASD and TD groups for any significant model (two-sample t-test, p > 0.05). Connection-level model evaluations and visualizations can be seen in Supplementary Table 1 and Supplementary Figures S1-S15 [Additional File 1].

Model reliability

A stratified 5-fold split-sample evaluation revealed high levels of curve reproducibility across all significant models. On average, the curve correlation for TD was 0.966 ± 0.024 while the curve correlation was 0.962 ± 0.028 in the ASD group. There were no differences in reliability between ASD and TD groups (t(28) = 0.779, p = 0.442). Reliability results on the connection-specific level can be seen in Supplementary Table 2 [Additional File 1].

Comparison to prior methods

Linear models analyzing the main effects of age, group, and their interaction yielded 57 brain connections that changed significantly with age, regardless of group (p-corrected < 0.05). No linear model group differences nor model interactions were significant following correction for multiple comparisons.

Binned age ANCOVAs, with ranges < 11, 11–18, and > 18, analyzing age, group, and their interaction yielded 30 brain connections that had a significant effect of binned age (p-corrected < 0.05). Additionally, 47 brain connections had a significant main effect of group (p-corrected < 0.05). However, there were no significant age-by-group interactions following correction for multiple comparisons.

Replication with Yeo 7 parcellation

Brain parcellation schemas may impact relationships between functional network connectivity and age [52, 53]. As such, we report a repeated GAM analysis with the larger-scale Yeo 7-network parcellation schema to examine the reproducibility of between-network findings. Of the 21 analyzed brain functional connections represented in the 7-network parcellation, seven changed significantly with age in the TD group (p-FDR < 0.05). Six of these seven represented connections between the limbic network and all other networks in the brain (visual, SMN, VAN, DAN, CEN, and DMN). Connectivity between these regions and the limbic network non-linearly increased in childhood and adolescence followed by stabilization or decreases in connectivity in adulthood, mirroring the more granular TD between-network limbic connectivity findings reported with the 17-network parcellation. Additionally, VAN-DAN between-network connectivity curvilinearly increased with age in the TD group (p-FDR < 0.05).

No functional brain connections significantly changed with age in the ASD group following statistical correction for multiple comparisons (p-FDR > 0.05). Interestingly, most of these non-significant ASD curve fits were linear, also reflecting increased linearity in ASD compared to TD reported with the 17-network parcellation.

Discussion

Using a flexible GAM framework, we found varying nonlinear connectivity-age relationships in typically developing (TD) individuals both within and between functional brain networks. Moreover, we found network-dependent coherence and deviations from these non-linearities in autism spectrum disorder (ASD). Specifically, the ASD group often followed typically developing connectivity-age-related pathways in between-network connectivity to the central executive network (CEN). Contrarily, we found chronic hypo-connectivity of multiple functional networks to the default mode network (DMN). We identified within-network somatomotor (SMN) hypoconnectivity in ASD, which specifically emerged in late adolescence. Lastly, we described fully disrupted age-related pathways in ASD localized to somatomotor, ventral attention (VAN), and limbic networks. ASD and TD coherence and deviations depend on network sub-parcellations. Taken together, our findings highlight the importance of flexible, non-linear modeling to describe divergent age-related trajectories in ASD and describing group differences that may emerge during certain periods of development.

Our findings replicate some from previous literature. Typically developing within-network connectivity has been reported to largely follow inverted U shape aging curves that decline or remain stable beginning in late adolescence [16, 18, 19, 21]. We report similar trends with ventral attention, dorsal attention, limbic, and central executive within-network connectivity. We break from previous findings in TD between-network connectivity changes with age. While previous studies have reported initial decreases followed by increases in DMN-CEN and DMN-dorsal attention (DAN) connectivity from childhood to early adolescence in typical development, we report steady increases during this period [16, 18]. These differences may be due to lack of harmonization despite the use of multisite data in these previous studies. Indeed, harmonization has been shown to improve power, sensitivity, and reliability in detecting connectivity-age relationships [54, 55]. Additionally, differences in our findings may be due to different parcellation schemas used in these studies which did not subdivide networks into multiple components as reported in the present study. Despite these differences in parcellation, our typically developing SMN between-network connectivity results mirror those reported by others, following cubic and sinusoidal curves [16, 17, 20, 21].

While two previous studies examined quadratic regression models of functional connectivity aging in ASD, their second-order polynomial approach may have resulted in limited findings [9, 25]. Indeed, previous reports have found complex non-linear divergent age-related changes in structural brain volume and cortical thickness in ASD [56,57,58]. Here, we report nonlinear functional connectivity across a multitude of functional networks in ASD. In doing so, we describe congruent and divergent age-related trajectories of functional connectivity in ASD.

DMN between network hypoconnectivity has been reported in numerous previous studies [5, 6, 8, 11, 59]. In line with these reports, we show between-network DMN-VAN and DMN-SMN hypo-connectivity in ASD that persists from childhood to adulthood. Importantly, the ASD group follows a similar inverse U-shaped aging trajectory in these connections with a constant decreased connectivity. Contrary to previous findings, we also report that age-related trajectories of between-network CEN connectivity with two other subcomponents of the DMN are congruent with TD individuals. These findings suggest that DMN between-network differences in ASD are not universal but may be localized to specific processes within the DMN. Indeed, previous studies of children with ASD have shown hyperconnectivity between the posterior cingulate and temporal cortex compared to TD, with simultaneous hypoconnectivity between the precuneus and visual cortex and basal ganglia [60]. Our results, combined with previous work, suggest between-network DMN differences in ASD may be sub-network specific.

The default mode network is implicated, among other processes, in self-referential thinking, passive sensory processing, and ‘theory of mind’ -- the ability to infer emotions and mental states of others [61,62,63]. Common ASD symptomatology includes sensory sensitivity, disruption of sensory integration and processing, and motor deficits [64, 65]. DMN hypo-connectivity to the somatomotor network implies a disconnect between active and passive sensory processing which is potentially relevant to sensory processing differences in ASD. Likewise, chronic DMN hypoconnectivity to the salience network may relate to theory of mind deficits in ASD, as detection of salient external events in reference to internal self-thought would be disrupted [66].

We also report divergent connectivity both within and between somatomotor network subcomponents in ASD. Within SMN connectivity follows a similar sinusoidal path in both ASD and typical development through childhood, but the ASD group departs from this path in adolescence and displays underconnectivity. Two recent mega-analyses reported within-network sensory-motor hypoconnectivity [5, 6]. Further, these hypoconnectivity directly correlated with self-reported sensory deficits and restricted social behavior [5]. While neither study reported age effects on these findings, they also implemented linear models in their analysis. Here, we similarly find within-SMN hypoconnectivity in ASD but report a non-linear aging pathway with hypoconnectivity specifically emerging in adolescence. Further, we find divergent aging of SMN between-network connectivity to the DAN, limbic, and CEN networks in ASD. Cubic-resembling SMN connectivity changes with age in typical development have been previously reported [16, 17, 20, 21]. Our results replicate this finding, while simultaneously revealing a completely divergent positive linear age-related pathway in ASD. As previously noted, common ASD symptomatology includes sensory sensitivity, disruption of sensory integration and processing, and motor deficits [65, 66]. Disrupted connectivity from the SMN to the DAN and CEN may contribute to decreased multi-sensory integration and sensory filtering, potentially contributing to sensory sensitivities and overstimulation often reported in ASD [67].

Divergent linear age-related pathways in ASD are also evident in between-network connectivity to the DAN. Interestingly, these divergent linear relationships are only present in connection with the posterior DAN containing the middle frontal gyri, supramarginal gyri, extrastriate cortex, and parietal cortex. In contrast, between-network connectivity to the dorsal-attention related somatomotor association network containing the fusiform and angular gyri shows congruent aging to TD. These results, along with DMN parcel dynamics, indicate that network parcellation size and schematics may influence the ability to detect specific brain connections exhibiting divergent age-related trajectories with GAMs. Indeed, prior studies have shown that parcellation schemas may impact relationships between functional network connectivity and age [52, 53].

To this end, we repeated our GAM analyses using a 7-network parcellation, which collapses the network subcomponent in our main analysis into seven large-scale network parcellations: visual, VAN, DAN, SMN, CEN, and DMN [43]. Using this large-scale parcellation, we replicate some TD aging pathways, largely pertaining to between-network limbic connectivity. Likewise, in the 7-network parcellation, the ASD aging paths and ASD-TD group relationships resembled those subsets of fully disrupted pathways described with 17-network framework. However, ASD replications and group differences with the 7-network parcellations failed to reach statistical significance after correction for multiple comparisons.

Our results highlight the need for flexible, non-linear modeling when analyzing the relationship between age and functional connectivity. To further demonstrate this, we repeated common linear model and age-binned ANCOVA analytical approaches to our harmonized 17-network dataset [5, 6, 8, 10, 11, 15]. Linear models revealed 42% of brain connections that changed significantly with age but yielded no significant group differences or group interactions following statistical corrections. Differently, the ANCOVA approach, where ages were binned at < 11, 11–18, and > 18, yielded fewer significant age-related changes in connectivity but resulted in significant group differences in at least one age bin in 35% of the connections.

These results highlight strengths, pitfalls, and the need for improvement in these common approaches. Linear models examine age as a continuous variable, allowing detection of connectivity changes across the age range, but are rigid and inadaptable to varying group differences that appear in some age ranges, but not others. Meanwhile, binned ANCOVAs are less powered to detect age-related changes in connectivity, as they compare changes across few binned groups. However, ANCOVAs allow group differences to be parsed during different age ranges, albeit pre-defined bins in which connectivity is known to dynamically change with age [12,13,14].

Generalized additive models address and improve upon both previous methods by maintaining age as a continuous variable while simultaneously allowing non-linear group differences to be parsed at any age range along the fitted curve. In doing so, GAMs can resolve the timing in which aging curvatures depart. Combined with increased sensitivity to inflection points compared to regression models, GAMs serve to define critical periods of connectivity divergence among groups.

Recent studies of cortical brain aging have employed general additive models for location scale and shape (GAMLSS) to explore brain changes throughout the human lifespan [68,69,70,71]. GAMLSS, an extension of GAMs, increases model flexibility by incorporating variance (scale) and distribution (shape) terms [72]. GAMLSS are particularly useful in cases of complex distribution with inhomogeneous variance and have most often been applied to define centiles in normative aging [68,69,70,71].

Given that the models currently do not violate residual homoscedasticity and normality assumptions, we currently employ the GAMs, not GAMLSS, for model simplicity and ease of interpretation. Further, the additional parameters in GAMLSS pose challenges with data sparsity, which is present at older ages in the ABIDE sample. In GAMLSS, errors are calculated on both mean and variance functions, which in sparse data can artificially inflate or deflate confidence intervals. For instance, if variance is underestimated due to sparsity, GAMLSS prediction intervals can be misleadingly tight which would lead to false detections of group differences in sparse regions. GAMs, whose errors depend on the mean function alone, are more transparently interpretable in this sparse scenario. It is important to note that GAMs’ penalized least squared splines limits model curvature inflection in sparse data, combatting overfitting at the potential loss of detecting small effects. Detection of smaller effects at older ages will be enabled by future data on aging in autism, which we are currently collecting (R01MH132218).

We presently apply GAMs for a conservative, interpretable approach to analyzing aging differences in autism. Nevertheless, GAMLSS and other extended GAM frameworks should be considered when analyzing complex data distributions. Generalized additive model statistical frameworks can be extended to many biological use cases where multiple comparisons are common. GAMs offer a data-driven modeling approach which minimizes assumptions about unique biological features relationship to a variable of interest. In doing so, GAMs ease the burden of regression model optimization, which is often unrealistic and time consuming when performing hundreds to thousands of models in large biological analyses.

Limitations

The present analysis contains several limitations. First, it is important to study autistic females as well as autistic males as divergent connectivity in ASD has been reported to be sex-dependent [73, 74]; the ABIDE samples did not have a sufficient number of females for our analyses. Second, we analyzed a restricted age range, ranging from childhood to mid-adulthood, thereby not allowing analysis across the full lifespan, including very early childhood and later in life brain aging. Future functional brain aging studies of older autistic adults, one of which is currently underway (R01MH132218), will help in the understanding of the biological bases of aging comorbidities associated with autism, such complex mental and physical health comorbidities which are currently understudied and ill-defined [75, 76]. Third, our analysis is limited to large scale network cortical functional parcellation. Subcortical functional connectivity may be fundamental to understanding aberrant aging in ASD [1, 9, 15]. To parse more specificity of brain region connectivity, a fine-grained functional parcellation including subcortical areas may be warranted. Fourth, post-hoc analyses of fitted curves were not corrected for multiple comparisons within each model. GAMs do not enable multiplicative group interactions as multiple trajectories are fit to different degrees of curvature. As such, between group differences rely on post-hoc comparisons of predicted values, which limits statistical power [72]. Finally, our age-related analyses used the cross-sectional data from ABIDE, rather than longitudinal data. As such, we focused on age-related cross-sectional ASD and TD group mean trajectories and were not able to quantify individual-level deviation. In the future, the quantification of an individual’s deviation from normative trajectories will be informative [77]. Using longitudinal data which we are currently collecting (R01MH132218), the GAM framework can extend to longitudinal within-individual analyses to parse deviations during aging across adulthood in autism.

Conclusion

Resting-state functional connectivity changes non-linearly with age in diverse curvatures across typically developing functional networks. Connectivity in ASD also follows nonlinear cross-sectional age-related trajectories that can be both divergent or congruent with typical development. Generalized additive models are an advanced statistical modeling tool to enable flexible non-linear aging models while simultaneously parsing group differences that appear during brain developmental and maturational periods in autistic individuals.