Abnormal autonomic and associated brain activities during rest in autism spectrum disorder

Skin conductance response acquisition

SCR was acquired according to the procedure described by Fan et al. (2012). Briefly, SCR was recorded using the GSR100C amplifier (BIOPAC Systems), together with the base module MP150 and the AcqKnowledge software (version 3.9.1.6). The GSR100C measures skin conductance by applying a constant voltage of 0.5 V between two electrodes that are attached to the skin. This allows for the measurements of both skin conductance level and SCR, which vary with sweat gland activity due to stress, arousal or emotional excitement. Skin conductance (measured in μS) was recorded using a 2000-Hz sampling rate (gain = 2 μS/V, both high-pass filters = DC, low-pass filter = 10 Hz). After cleaning the skin with alcohol swabs, two EL507 disposable EDA (isotonic gel) electrodes were placed on the palmar surface of the distal phalanges of the big and second toes of left foot. The electrode leads were shielded and the signal was low-pass filtered (using the MRI-Compatible MRI CBL/FILTER System MECMRI-TRANS) to reduce radiofrequency interference from the scanner. BIOPAC recording was synchronized to the E-Prime (Psychology Software Tools) program through the parallel port of the computers to enable precise time alignment of skin conductance recording with scan onsets.

Image acquisition

All brain images were obtained using a 3 T Siemens Allegra MRI system at the Icahn School of Medicine at Mount Sinai, and were acquired parallel to the anterior-posterior commissures axis (AC-PC). Foam padding was used to reduce head motion. Whole-brain anatomical T₂-weighted images were acquired in high-resolution using a turbo spin-echo plus sequence: 40 axial slices of 4 mm thickness; skip = 0 mm; repetition time = 4050 ms; echo time = 99 ms; flip angle = 170°; field of view = 240 mm; matrix size = 448 × 512; voxel size = 0.47 × 0.47 × 4 mm. After the anatomical image acquisition, one run of T₂*-weighted images was obtained during rest, corresponding to the T₂-weighted images localization, using a 6 min gradient echoplanar imaging sequence for resting-state functional MRI: 40 axial slices, 4-mm thick; skip = 0 mm; repetition time = 2500 ms; echo time = 27 ms; flip angle = 82°; field of view = 240 mm; matrix size = 64 × 64; in-plane resolution = 3.75 × 3.75 mm. The resting-state functional MRI run started with two dummy volumes before the onset of the fixation to allow for equilibration of T₁ saturation effects, followed by 144 image volumes. Refer to the Supplementary material for a detailed procedure of the eyes-open resting-state scan.

Skin conductance analysis

Skin conductance level was calculated as the average of all data points on the skin conductance waveform for each of the participants. A t-test was conducted between the skin conductance level values of the groups to examine possible differences in basal skin conductance levels. AcqKnowledge software (version 4.2) was used in order to identify and count SCRs (Supplementary material). The number of valid SCR events was determined for each participant, and a t-test was then conducted to examine possible group differences in the number of SCRs.

Regression analysis

To correlate SCR with brain fluctuations, the skin conductance waveform was down-sampled by averaging the data points in each 2.5-s bin matching the repetition time (2.5 s) of the echoplanar imaging scan of functional image acquisition. Because relaxation causes skin conductance level to decrease slowly in a linear fashion, the SCR waveform was detrended and band-pass filtered with the same frequency range (0.01–0.08 Hz) that is used in a typical resting-state functional MRI analysis. The similarity between the SCR and haemodynamic response function curves (Bach et al., 2010a, b) allows for using SCR as a regressor in the model directly without transformations, using the same filtering band for both SCR and blood oxygenation level-dependent signals (0.01–0.08 Hz) (Fan et al., 2012). SCR was then entered as a regressor in a general linear model (Friston et al., 1995) using statistical parametric mapping package (SPM8; Wellcome Trust Centre for Neuroimaging, London, UK). Numbers of SCR were not equated between the groups to avoid specification error (Supplementary material). The echoplanar imaging scans were realigned to the first volume, timing corrected, coregistered to the T₂ image, normalized to a standard template (MNI, Montreal Neurological Institute), resampled to a 2 × 2 × 2 mm voxel size, and spatially smoothed with an 8 mm full-width at half-maximum Gaussian kernel. The general linear model was then conducted with the SCR time series as a predictor of the observed blood oxygen level-dependent signals. Low-frequency drifts in signal were removed using a standard high-pass filter with a 128 s cut-off. Serial correlation was estimated using an autoregressive AR(1) model. To remove non-neural noise from the data, ventricle and white-matter signals were extracted using corresponding masks and were entered as covariates. In addition, the six parameters generated during motion correction were also entered as covariates. Head motions (for all participants) did not exceed 2.5 mm of displacement or 2.5° of rotation in any direction. Finally, mean voxel value was used for global calculation and grand mean scaling, applied with global normalization to further remove non-specific noise (Van Dijk et al., 2010). Detailed justifications for using the global mean correction can be found in the Supplementary material.

The contrast images from the participants in each group were entered into a second-level random effect group analysis. The resultant voxel-wise statistical maps were thresholded for significance using a cluster-size algorithm that protects against an inflation of the false-positive rate of multiple comparisons; an uncorrected P-value of 0.05 for the height (intensity) threshold of each activated voxel and extent threshold of k = 120 was used based on Monte Carlo simulation (Slotnick and Schacter, 2004). Assuming an individual voxel type I error of P < 0.05, a cluster extent of 120 contiguous resampled voxels (2 × 2 × 2 mm) was indicated as necessary to correct for multiple voxel comparisons at P < 0.05. Statistical results were mapped onto the standardized surface of the cerebral cortex.

Functional connectivity analysis of the default mode network

To investigate the relationship between non-specific SCR and brain connectivity during rest, a functional connectivity analysis was conducted using a seed region within the posterior cingulate cortex as in previous studies (Koshino et al., 2005). Specifically, time-series volumes of functional MRI scan images were preprocessed for each participant using the Data Processing Assistant for Resting-State fMRI (DPARSF) toolbox (Yan and Zang, 2010). This included slice timing correction, realignment, coregistration, normalization and spatial smoothing (using an 8 mm full-width at half-maximum Gaussian kernel) as in the regression analysis. In addition, de-trending (to remove the systematic drift) and temporal filtering (band-pass, 0.01–0.08 Hz, to reduce the effect of low-frequency drift and high-frequency physiological signal or noise) were applied. Time course of the posterior cingulate cortex (left and right combined) was then extracted using the automated anatomical labelling template (Tzourio-Mazoyer et al., 2002), and a voxel-wise linear correlation between the mean time course of the posterior cingulate cortex and the time course of each voxel in the whole brain was calculated using the resting-state functional MRI data analysis toolkit (Song et al., 2011). The six head motion parameters, global mean signal, white matter signal, and CSF signal were included as covariates. A two-sample t-test was then conducted to examine between-group connectivity differences.

To test the effects of SCR on brain connectivity, changes in functional connectivity of the posterior cingulate cortex before and after regressing out SCR effects were examined and compared between the groups in an ANOVA model. SCR signals were regressed out as covariates in a second voxel-wise connectivity analysis, followed by transformations of the correlation coefficients using Fisher’s r-to-z’ transformations. Posterior cingulate cortex time course was also extracted after regressing out the SCR signal. In addition, to examine SCR contributions to the posterior cingulate cortex functional connectivity within each group, a paired t-test was conducted for each group before versus after regressing out SCR signals for the posterior cingulate cortex functional connectivity maps. Posterior cingulate cortex connectivity analysis is reported here because it allows comparisons with previous studies investigating DMN connectivity in ASD that used a posterior cingulate cortex seed (Monk et al., 2009; Weng et al., 2010), and because the posterior cingulate cortex has a more focal anatomical definition than alternative seed regions. However, the ventromedial prefrontal cortex, also commonly used as a DMN seed, is potentially more relevant to ASD abnormalities, emotional processing and autonomic pathways. We, therefore, conducted an additional functional connectivity analysis for the time course of the ventromedial prefrontal cortex, using the coordinates (−1, 47, −4) reported by Fox et al. (2005), and group differences were examined. To further examine possible differences in the connectivity of areas that are related to emotional and autonomic signal processing, we also tested the functional connectivity of the anterior insular cortex (Supplementary material). It is important to note that there were no significant differences in head motion between the two groups, as was indicated by the root mean squares of both overall head motion displacement and rotation and the temporal derivatives (Supplementary material).

Functional connectivity analysis of the whole brain

To examine whether SCR contributes to the whole brain connectivity, a voxel-wise whole brain functional connectivity strength analysis was performed both before and after regressing out the SCR signal. Images were preprocessed using DPARSF, similar to the posterior cingulate cortex functional connectivity analysis, with two exceptions. First, the resolution of resultant echoplanar imaging was 3 × 3 × 3 mm after normalization to reduce computational load. Second, to reduce artificial local correlations between voxels introduced by smoothing, a 4 mm full-width at half-maximum Gaussian kernel was used (instead of 8 mm as in the posterior cingulate cortex connectivity analysis). For each voxel, functional connectivity strength was measured as the summed weights of all connections linking this voxel and every other voxel. Pearson correlation coefficients for the time series of every possible pair of voxels were calculated to obtain the whole brain correlation matrix for each participant. The calculation was constrained within a customized mask (n voxels = 48 159) including all voxels with grey matter tissue probability >20% on the averaged grey matter map of all participants. The functional connectivity strength (FCS) was computed for each voxel as follows:

where r_ij is the correlation coefficient of voxel i and voxel j, z_ij is the normalized r_ij value using Fisher’s r-to-z transformation, and 0.3 is the threshold set to eliminate the potential contributions of weak connections arising from noise. To evaluate the reproducibility of our results, additional functional connectivity strength maps were calculated using both 0.2 and 0.4 correlation thresholds; these computations led to no major changes in our primary results. The connectivity map was then standardized by converting to Z scores so that maps across participants could be averaged and compared. The Z score transformation is:

where are mean and SD of the functional connectivity strength across all the voxels in the whole-brain map.

Notably, in graph theory, the functional connectivity strength is referred to as ‘degree centrality’ or ‘degree’ of weighted networks, and voxels with high functional connectivity strength usually play important roles in information transformation (Buckner et al., 2009; Cole et al., 2010; Zuo et al., 2012; Liang et al., 2013; Wang et al., 2013). To examine the SCR contribution to whole brain connectivity, functional connectivity strength maps before and after regressing out the SCR signals were compared within each group using paired t-tests. Between-group differences of SCR signal contribution to whole brain functional connectivity (i.e. functional connectivity strength maps before versus after regressing out SCR signals) were examined in an ANOVA model. All statistical maps were corrected for multiple comparisons to a significance level of P < 0.05 by combining the individual voxel P-value < 0.05 with cluster size >67 voxels, based on Monte Carlo simulation (Ledberg et al., 1998).

Relationship between autism symptom severity, skin conductance and imaging data

Correlations between the ADI-R subscale scores of each of the ASD participants and number of SCRs and overall skin conductance level were examined. In addition, correlations between the individual ADI-R subscale scores and SCR-brain association, posterior cingulate cortex connectivity, ventromedial prefrontal cortex connectivity, and whole brain connectivity maps were examined using regression analyses. The ADI-R four subscales include ratings of social interaction (subscale A), communication (subscale B), repetitive and restricted behaviour (subscale C), and early development (subscale D). Higher scores on each scale indicate higher symptom severity in each of the specific domains. Note that these exploratory correlation analyses were not corrected for multiple comparisons.

Modelling central generation and representation of skin conductance response

The central generation and representation (feedback) of SCR were modelled by shifting the SCR vector in relation to the haemodynamic response. Two models were created: a generation model in which brain activation preceded the SCR signal by one repetition time (2.5 s), and a representation model in which SCR preceded brain activation by one repetition time (Critchley et al., 2000; Patterson et al., 2002) (see Supplementary material for details of the modelling).

Posterior cingulate cortex seed-based analysis results

Using the posterior cingulate cortex as a seed region, the neurotypical control group had stronger connectivity between the seed and areas of the DMN, compared to the ASD group, mainly with the ventromedial prefrontal cortex and the left inferior parietal lobule (Table 5 and Fig. 3A). Stronger connectivity was found in ASD between the posterior cingulate cortex and the superior and inferior occipital gyri, temporal pole, lingual gyrus, fusiform gyrus, amygdala, hippocampus, posterior insular cortex, supplementary motor area, and precentral and postcentral gyri (Table 5 and Fig. 3B). The general linear model revealed a significant interaction between group (ASD versus neurotypical control subjects) and SCR (before versus after regressing out SCR), wherein SCR had greater positive effects on the connectivity between posterior cingulate cortex and medial prefrontal and orbitofrontal cortices, precentral and postcentral gyri, supramarginal gyrus, and superior and inferior parietal gyri in neurotypical control subjects than in ASD (Table 6 and Fig. 3C). In contrast, SCR had a greater positive impact on the connectivity between the posterior cingulate cortex and the posterior insular cortex, lingual gyrus, amygdala, hippocampus, parahippocampal and fusiform gyri, superior and middle occipital gyri, and superior and inferior temporal gyri in ASD than in neurotypical control subjects (Table 6 and Fig. 3D). For paired t-test results within each group (before > after regressing out SCR signal), see Supplementary Table 1 and Supplementary Fig. 1. For ventromedial prefrontal cortex and anterior insular cortex functional connectivity results, see Supplementary Tables 2 and 3, and Supplementary Figs 2 and 3.

Voxel-wise whole brain connectivity

A between-group comparison revealed that the neurotypical control group had higher functional connectivity strength in the posterior cingulate cortex, precuneus, medial prefrontal cortex, anterior cingulate cortex, gyrus rectus, posterior insular cortex, superior temporal gyrus, and inferior temporal, parietal, and occipital gyri (Table 7 and Fig. 4A). The ASD group had higher functional connectivity strength in the hippocampus, inferior temporal gyrus, temporal pole, lingual gyrus, calcarine cortex and amygdala, all on the left, as well as the cerebellum bilaterally (Table 7 and Fig. 4B). A significant interaction (Group × SCR) was also found; SCR signal contributed significantly to the connectivity of the posterior cingulate cortex, medial prefrontal cortex, anterior cingulate cortex, inferior parietal lobule, and the right anterior and posterior insular cortex in neurotypical control subjects, compared with ASD (Table 8 and Fig. 4C). In ASD there was significant contribution of SCR signal to the connectivity of the hippocampus, amygdala, nucleus accumbens, calcarine cortex, and the fusiform, lingual, and inferior occipital gyri, compared with neurotypical control subjects (Table 8 and Fig. 4D). For within-group paired t-test results (before > after SCR regression), see Supplementary Table 4 and Supplementary Fig. 4.

Central generation and representation of skin conductance responses

One-sample t-tests of both generation and representation models revealed comparable brain activation associated with SCR to the original data in both groups (i.e. anterior brain regions in neurotypical controls, and posterior brain activation in ASD). When the two models were compared within each group, no significant differences were found between representation and generation (representation > generation) in neurotypical control subjects, or between the two models in ASD. There was, however, significant activation in the right thalamus and the left cerebellum in the generation model, as compared with the representation model (generation > representation), in neurotypical control subjects (Supplementary Table 9 and Supplementary Fig. 9).