Autistic traits in neurotypical adults: correlates of graph theoretical functional network topology and white matter anisotropy patterns.

Attempts to explicate the neural abnormalities behind autism spectrum disorders frequently revealed impaired brain connectivity, yet our knowledge is limited about the alterations linked with autistic traits in the non-clinical population. In our study, we aimed at exploring the neural correlates of dimensional autistic traits using a dual approach of diffusion tensor imaging (DTI) and graph theoretical analysis of resting state functional MRI data. Subjects were sampled from a public neuroimaging dataset of healthy volunteers. Inclusion criteria were adult age (age: 18-65), availability of DTI and resting state functional acquisitions and psychological evaluation including the Social Responsiveness Scale (SRS) and Autistic Spectrum Screening Questionnaire (ASSQ). The final subject cohort consisted of 127 neurotypicals. Global brain network structure was described by graph theoretical parameters: global and average local efficiency. Regional topology was characterized by degree and efficiency. We provided measurements for diffusion anisotropy. The association between autistic traits and the neuroimaging findings was studied using a general linear model analysis, controlling for the effects of age, gender and IQ profile. Significant negative correlation was found between the degree and efficiency of the right posterior cingulate cortex and autistic traits, measured by the combination of ASSQ and SRS scores. Autistic phenotype was associated with the decrease of whole-brain local efficiency. Reduction of diffusion anisotropy was found bilaterally in the temporal fusiform and parahippocampal gyri. Numerous models describe the autistic brain connectome to be dominated by reduced long-range connections and excessive short-range fibers. Our finding of decreased efficiency supports this hypothesis although the only prominent effect was seen in the posterior limbic lobe, which is known to act as a connector hub. The neural correlates of the autistic trait in neurotypicals showed only limited similarities to the reported findings in clinical populations with low functioning autism.

Introduction

The American Psychiatric Association defines autism as a pervasive developmental disability with impairments in social interaction, verbal, nonverbal communication and limited behavioral flexibility [1]. Conventionally, if an individual with autism has an intelligence quotient in the normal range (or above), the term high-functioning autism (HFA) is used, while meeting all of the criteria for HFA except communicative abnormality or history of language delay defines Asperger's syndrome [2]. Neuroanatomic observations provided evidence about micro- and macrostructural abnormalities in the autistic brain [3] and in the quest for characterizing such malformations, neuroimaging gained prominent attention. MRI volumetry revealed abnormal developmental pattern in autism, which is characterized by an early brain overgrowth that diminished in adults [4], [5]. Further evidence on localized structural abnormalities was provided by voxel-based morphometry studies [6]–[8] and cortical thickness measurements [9], [10]. Distinct anatomic correlates were linked with Asperger's Syndrome [11] and high-functioning autism [12]. Technological developments in the last two decades enabled the use of more advanced structural and functional magnetic resonance approaches. Here we shortly summarize the contributions of diffusion magnetic resonance and functional MR (fMRI) to autism research.

Diffusion-weighted imaging (DWI) is based on the sensitization of the magnetic resonance signal to molecular displacements while diffusion tensor imaging (DTI) acquires information on the directional preference of the anisotropic diffusion in the brain [13]. The driving force behind DWI's 25-year success story is that opened a new window on the characterization of tissue microarchitecture by probing the microscopic, diffusion-driven molecular motions. This approach gave rise to diverse studies that elucidated the correlations between regional changes of diffusion anisotropy and various psychopathologies [14] or linked the diffusion characteristics with cognitive phenotypes [15]. Lowered quantitative measures, such as the DTI-derived diffusion anisotropy, were reported in the white matter of autistic children's brains [16], [17] or in adults with Asperger's syndrome [18]. Inconsistently, some studies report regionally increased fractional anisotropy (FA) in diverse brain regions [19], [20].

Functional magnetic resonance imaging (fMRI) is an optimal modality for the in vivo observation of the working brain and to infer patterns of brain activity from the changes of the blood oxygen-level dependent (BOLD) MR signal. In autism research putative functional impairments can be characterized by using experimental tasks such as executive function tests [21], cognitive control [22] or verbal fluency tasks [23]. The endeavor to study the normally functioning or pathologically altered macroscopic functional networks of the human brain was recently revitalized by the resting state fMRI concept [24]. In such studies the fundamental hypothesis is that the temporal coherence of BOLD signal time-courses marks the degree to which two regions are functionally connected [25], although this simple model does not imply causality and further inconsistency of such reports is caused by the large number of signal and image processing approaches that are in use [26]. More complex models can be used to examine the brain as a network of interconnected processing units rather than interpreting individually measured functional connections. One such mathematical approach is graph theory where brain regions define graph nodes and edges represent the strengths of structural or functional connectivity. This emerging network-based technique allows the calculation of graph topological measures that hallmark the properties of information exchange. The most commonly used indicators are small-worldness, efficiency, degree, modularity. Such quantitative studies are increasingly recognized – we refer to review articles on the background and possible clinical applications of the graph theoretical analysis of the human brain connectome [27], [28].

Resting state fMRI measurements, diffusion tensor imaging and tractography are potentially applicable to indirectly explore the disturbed neuronal communication in autism. Our study aims to stand in the line of studies that assume abnormal formation of cortico-cortical short or long-range interconnections. Such models postulate impaired long-range connectivity in autism and in HFA which is also hallmarked by decreased diffusion anisotropy that predominantly affects frontal, fronto-parietal and fronto-occipital pathways [29]. Conversely, a minority of studies revealed increased connectivity values or diffusion anisotropy [30], [31] particularly coinciding with the spatial patterns of short-range association fibers [19], [32]. The reported inconsistencies in neuroimaging studies reflect the puzzling nature of autism in which the core anatomical features and the common pathomechanisms of disrupted connectivity remains unknown.

Our investigation was designed to elucidate the putative effects of the autistic traits on the patterns of functional brain connectivity and diffusion anisotropy. Autism spectrum disorders are commonly understood to comprise traits with dimensionality; such traits can plausibly be quantified by clinical batteries or self-reported screening tests [33]. Accordingly, we assume that normally functioning adults with a certain degree of social reciprocity impairment and inflexible behavior present altered structural and functional brain connectivity. This supposition prompted us to employ a graph theoretical analysis that is presumably feasible and sensitive for characterizing changes in the functional topology of resting state networks (RSNs). We aim to relate the RSN topological correlates of autistic traits to the commonly reported observations from the broader spectrum of the disease.

Methods

Study population

Phenotypic information and imaging data of 207 subjects were taken from the public repository of the International Neuroimaging Data-sharing Initiative (INDI); we used the most recent release of the Nathan Kline Institute's Rockland Sample [34]. It is a freely available, large-scale, extensively phenotyped dataset for the discovery research and it contains healthy subjects from nearly all age groups. During sampling, inclusion criteria were: (i) age between 18 and 65, (ii) availability of the Wechsler Abbreviated Scale of Intelligence (WASI), Social Responsiveness Scale (SRS) and the autism spectrum screening questionnaire (ASSQ) scores, (iii) availability of DTI and resting-state fMRI scans, (iv) IQ above 70. As a result of this selection, we performed our analysis on a cohort of 127 neurotypical adults. The original NKI/Rockland data collection was preceded by the approval of the relevant ethical committee and subjects gave informed consent to the imaging studies.

Psychological evaluations

Our hypothesis is that inter-individual variability in dimensional autistic traits is correlated with characteristic changes in the functional network topology and diffusion anisotropy in neurotypical adults. The study was consequently designed to utilize multiple psychological tests as explanatory variables for functional topological measurements. First, we accessed data of the Autism Spectrum Screening Questionnaire (ASSQ), which is a 27 item checklist for completion by lay informants. This battery was reported to be plausible in assessing symptoms characteristic of Asperger's syndrome and other high functioning autism spectrum disorders in children and adolescents with normal intelligence [35]. The ASSQ was self-administered by the adults participating in this study. Data of self-reported autistic traits were completed by using the Social Responsiveness Scale (SRS) Adult version. The SRS is a 65 item questionnaire originally intended for an informant, but it was completed by the subjects here. Total SRS scores range from 0, corresponding to high social competence, to 195, corresponding to severe social impairment as observed in individuals with severe ASD. Scores between 60 and 80 are associated with mild forms of ASD [36]. SRS measures autistic traits that are continuously distributed in the normal population and has been used in similar investigations that correlate such traits with functional neuroimaging based data [37]. HFA is occasionally characterized by discrepant IQ profile, in other words, performance IQ scores significantly exceeding verbal IQ [38], but this observation is not ubiquitous [39]. To evaluate the intelligence profile, the authors of the NKI/Rockland dataset utilized the WASI, which consists of four subtests [40]. The performance IQ is composed of the scores of two subsets: the Block Design and Matrix reasoning, while the verbal IQ comprises the Vocabulary and Similarities tests. The distribution of the above mentioned demographic variables were tested for normality with one sample Kolmogorov-Smirnov (K-S) procedure. WASI verbal and full IQ scores were normally distributed (K-S; p = 0.57, 0.81). Age was not found to be normally distributed due to the lower and upper thresholds applied (K-S; asymptotic significance; p = 0.031). SRS total score was normally distributed (K-S; p = 0.23) while ASSQ showed significant deviation from the assumed distribution (K-S; p<0.001). Table 1 summarizes the basic demographic data and trait-based measures for the study group. Figure 1 illustrates the distribution of age and the psychological assessments.

Table 1

Demographic data of the study population and summary of psychological evaluations.

	Total population	Females	Males
n	127	49	78
Age	35.8±12.8 (18–65)	36.3±13.2 (19–64)	35.5±12.6 (18–65)
Full Scale IQ	108.3±13.1 (74–137)	107.9±13.4 (80–137)	108.5±13.1 (74–136)
Verbal IQ	106.9±13 (67–134)	106.4±12.9 (79–129)	107.1±13.2 (67–134)
ASSQ Score	5.2±7.9 (0–39)	4±5.9 (0–22)	6±8.9 (0–39)
SRS Total Score	27.2±19.4 (0–109)	22±14.6 (1–71)	30.4±21.4 (0–109)

Data are given in mean ± standard deviation (range: minimum – maximum) format. ASSQ: Autism Spectrum Screening Questionnaire, SRS: Social Responsiveness Scale.

Figure 1

Demographic data of the study population.

Distribution of subject age, full scale IQ, verbal IQ, ASSQ and SRS scores. Data are depicted as histograms (left panels) and Q-Q (quantile-quantile) probability plots (right panels) in which reference lines of the normal distribution are given (n = 127).

Imaging protocols

DTI and fMRI scans were acquired on a 3.0 T MRI system (Magnetom Trio Tim, Siemens, Erlangen, Germany) using a spin echo echo-planar imaging (EPI) sequence (TR = 10000 ms, TE = 91 ms) with the GRAPPA parallel imaging technique (acceleration factor: 3). Diffusion-weighting gradients were applied in 64 different directions, b-value: 1000 s/mm2. Volumes consisted of 58 transverse slices, slice thickness: 2 mm, voxel size: 2 mm * 2 mm, matrix size: 128 * 128 (field of view: 256 mm). Resting-state functional magnetic resonance imaging protocol was adapted from the Brain Genomics Superstruct Project common acquisition protocol, provided courtesy of Randy Buckner. Over an experimental period of approximately 10 minutes, 260 volumes were acquired with a spin echo EPI sequence (TR = 2500 ms, TE = 30 ms). 38 axial slices were acquired with a matrix resolution of 72 * 72 (field of view: 216 mm), isotropic voxel size was 3 mm.

Image processing – spatial standardization

Diffusion tensor and fMRI data were standardized to the MNI152 neuroimaging space. For guiding the transformation of the DTI space into the standard space, the fractional anisotropy (FA) images were used. Initially, FA images were registered to the FMRIB58 fractional anisotropy template linearly (FLIRT algorithm, 12 degrees of freedom registration using normalized correlation cost function, included in the FSL software library [41]). This affine transformation served as an initial stage for the non-linear registration (FNIRT algorithm).

The rs-fMRI data of each participant were realigned to the middle scan of the series using rigid body registration (MCFLIRT algorithm in FSL) [42]–[44] to adjust for motion and then co-registered to the T1-weighted MPRAGE image using rigid-body algorithm (FLIRT, FSL) [45]–[47]. Data were then spatially standardized by applying the warp field produced by the nonlinear co-registration of the skull-stripped [48] MPRAGE scan to the MNI152 non-linear template [49], [50] using the FNIRT tool in FSL. The diffusion-weighted and T1-weighted scans were skull-stripped using the BET command in FSL.

Image processing – brain region system

A customized brain region system was defined to comprise the entire cortex and subcortical areas including the cerebellum. When determining the nodes of the network, Smith et al. showed that it is of crucial importance to ensure that the ROIs provide whole gray matter coverage and thus eliminate the deleterious effect of shared extra-network inputs [51].

For this purpose, we constructed a custom region system containing 149 regions, based on multiple brain atlases. Atlas data were available in MNI152 standard space. Regions were sampled from the Harvard-Oxford Cortical and Subcortical Structural [52]–[55] atlas while cerebellar labels were taken from the UCL (University College, London) probabilistic cerebellar atlas [56]. The thalamus region of the Harvard-Oxford atlas was parceled into seven sub-regions based upon the Oxford Thalamic Connectivity Atlas [57], [58], each subdivision representing thalamic connectivity domains to major cortical areas. Each atlas comprised voxel-wise probabilities of the volumetric label maps. Accessing this data, we delineated each ROI at p>0.25 probability threshold, using the in-house developed BrainLOC (Software Access: www.minipetct.hu/brainloc) software package. To eliminate regional overlaps, voxels belonging simultaneously to two or more ROIs were assigned to the one with the highest probability. We provide a complete list of label names, abbreviations and data sources in Table S1.

Image processing - diffusion tensor estimation

Before the non-linear standardization of DTI acquisitions, tensors were fitted to the observed DWI signal using the DTIFIT algorithm (FSL), this step included the calculation of fractional anisotropy (FA) values. For the equation describing the relationship between the tensorial eigenvalues and FA, we refer to the commonly used ways in the literature [59]. Atlas-based labels were used to extract the regional values of FA. To determine the mean anisotropy value for each region, anisotropy measurements in each image voxel were weighted with the corresponding atlas-based probabilities.

Image processing and network construction - resting state fMRI

We characterized the strength of temporal coherence (functional connectivity) between each pair of brain regions with the wavelet correlations of the averaged regional BOLD signals [60] on a low (<0.1 Hz) frequency band, which is most commonly used in rs-fMRI studies [61], [62]. For estimating wavelet correlations we applied four level maximum overlap discrete wavelet transformation [63], [64] with a Daubechies least asymmetric wavelets of length 8 (so called LA(8) wavelet) to calculate four sets of wavelet coefficients for each BOLD signal. The four level wavelet decomposition refers to four frequency bands (0.1–0.2 Hz, 0.05–0.1 Hz, 0.025–0.05 Hz, 0.0125–0.025 Hz) depending on the [21/(k+1)/TR – 21/k/TR, k = 1…4] rule and the TR = 2500 ms repetition time, we chose the coefficients related to the 0.025–0.05 Hz frequency range for the further analysis (Band 3). The wavelet correlation for this band was calculated from the wavelet coefficients as described in [63], [64].

Wavelet correlations (rij) between each possible pair of brain regions were computed and stored in 149*149 symmetric correlation matrices. In order to preserve the continuous nature of the correlation information we chose the weighted undirected network model for characterizing network properties [65]–[67]. To emphasize strong correlations and punish weak correlations we defined the connectivity between two regions as a power 2 of the absolute value of the wavelet correlation coefficients [68]: wij = rij 2. Using this soft threshold approach we produced 149*149 weighted connectivity matrices for each subject. Figure 2 summarizes the image processing pipeline that was utilized for the rs-fMRI data.

Figure 2

Processing pipeline of resting state functional MRI data.

Panel A, far left: two atlas regions, Inferior Frontal Gyrus, triangular part (red, R23) and Middle Temporal Gyrus temporooccipital part (blue, R58) were overlaid on an atlas-space aligned fMRI image. The corresponding regional BOLD curves and their wavelet coefficients are shown on next images. Accordingly, the TR = 2500 ms repetition time and the four level wavelet decomposition the evaluated wavelet coefficients of four frequency bands (0.1–0.2 Hz, 0.05–0.1 Hz, 0.025–0.05 Hz, 0.0125–0.025 Hz) were generated (middle bottom image). In the 4th column the evaluated wavelet correlation matrix (top) and the distribution of correlation coefficients (bottom) are shown. We used these wavelet correlation values between each pair of brain regions to construct the matrix in which yellow color represents high correlation coefficients while red means low values. Panel B: In the bottom row, three adjacency matrices are shown with wired costs 0.1, 0.4 and 0.9. These matrices were generated from the weighted connectivity matrix (4th panel, right bottom image) by different weight thresholds. Vertical gray arrows represent the calculation procedure of nodal and global network parameters at different cost levels. The horizontal gray arrows demonstrate the final step of Monte-Carlo based cost-integration procedure in which the summed parameters are divided by the integration steps (MC).

Graph theoretical analysis of functional connectivities

In graph theory, the mathematical abstraction of a network is a G graph containing a set of N nodes and M edges. This graph can be represented as an N*N square “A” adjacency matrix with elements aij = 1 or 0, depending on whether an edge does (1) or does not exists (0) between nodes i and j. In the case of weighted networks the W connectivity/weight matrix containing wij elements represents the strength of connections between nodes i and j. In MRI-based functional brain network analysis the adjacency and weight matrices are zero-diagonal, symmetric matrices (wij = aij = 0, wij = wji, aij = aji) since modeling regional self-connections (i = j) and directed regional connections (wij≠wji, aij≠aji) are meaningless. These matrix properties correspond to undirected binary (so called unweighted) or undirected weighted graphs, depending on the used edge model (binary, weighted). In the case of weighted functional brain network without losing generality we assume that the weights of edges (i.e. the strength of connections) lie in the unit interval, [0,1].

Measurement of graph sparsity

The number of maximum edges of an undirected graph MC = N(N−1)/2, i.e. the number of elements of the upper-triangle of adjacency or weight matrix. The C index stands for “complete graph” in which every pair of nodes are connected. In the case of binary undirected graphs theratio defines the edge density or wiring cost of G. It is obvious that this parameter lies in the unit interval, [0,1], while the number of edges of empty and complete graphs is 0 and MC, respectively.

Substitution M and MC into formula of K(G):and soresults in the fact that the wiring cost is the mean of values of the off-diagonal elements of the adjacency matrix. The adaptation of wiring cost definition of binary graphs to weighted graphs is simple [69], because the K(G) can be calculated as the mean value of the off-diagonal elements of weight matrix, if the weight values lie in the unit interval, with maximal value 1:

Network topology: global graph features

Using the described graph theoretical notation we can formulate a large set of parameters that characterize the global, nodal or modular properties of the investigated networks, as reviewed by Rubinov and Sporns [67]. It is known from the literature that the functional human brain network shows small-world characteristics, therefore in our study we focused on quantifying only those nodal and global parameters that may have an influence on this property. Latora and Marchiori [70] introduced the efficiency-based characterization of small-world networks, which, due to its computational benefits, was proved to be more effective than the classical characteristic path length and clustering coefficient based calculation methods [71], [72]. The efficiency was introduced as the measure of the effectiveness of information exchange between nodes, while the average efficiency of nodes of a G graph was defined as global efficiency (Eg), which gives a normalized measure (lies in [0,1]) of the information transfer efficiency of parallel systems. The normalized local efficiency (El) of the network measures how efficient the local communication between first neighbors of a node is if this node is faulted or removed (fault tolerance). Latora and Marchiori also showed that the small-world behavioral network has high global and local efficiency [70]. According to this assumption, we calculated two nodal (degree and efficiency) and two global (global and local efficiency) parameters for all subjects networks using the following formulas.

Nodal parameters of unweighted graphs

The k [67] is defined as the sum of edges of node i:The efficiency coefficient of the path (path efficiencies) between node i and j is defined as the inverse of the l between these nodes [67]. The shortest path length l between nodes i and j means the minimum of all possible path lengths between nodes i and j, and can be calculated from the adjacency matrix by Dijksra's [73] or Floyd's [74] algorithms. Using these definitions the efficiency of node i is defined as the average of the path efficiencies of the given node:

Global graph parameters of unweighted graphs

Global efficiency Eg is given as the mean of nodal efficiencies and gives a normalized measure (lies in [0,1]) of the information transfer efficiency of the network:Local efficiency of an unweighted graph is defined as the average global efficiency of sub-graphs of nodes, because it measures how efficient the local communication is between the first neighbors of a node if this node is faulted or removed:Where Gi denotes the subgraph composed from the neighbors of node i, and by definition does not contain node i.

Cost-integrated measurement of topological metrics

One of the critical points of functional brain network analysis is the connection (edge) selection or the thresholding procedure. The fundamental aim of this process is to eliminate weak connections from the graph representation of network by an arbitrarily chosen weight threshold. The result can be a weighted graph containing only edges with higher weight than the applied threshold or it can be an unweighted graph containing only thresholded edges. Thresholding will change the wiring cost of graphs, in other words, it changes the topology of the investigated network. In the functional brain network assays we have an illusory solution to avoid this problem if we use the evaluated, unthresholded primary networks for any population level analysis. However, the primary networks differ in their average correlation coefficients and they also differ in their mean weights and wired-costs [75], which makes difficult to compare their topology parameters. Although, the adaptive thresholding technique [76] could guarantee the same weighted or unweighted wired cost for all networks in the population by uniquely selected thresholds, but it generates new questions: which wired cost is the best for the analysis and how the selected cost affects the results? Since there is currently no consensus in the literature regarding selection of threshold for weighted graph analysis we applied a cost-integrated technique as introduced by Ginestet et.al. in [69]. Following the definition of this paper we calculate the cost-integrated values of any X topological parameter of the graph G by this formula:where k is a cost value from the Ω set (set of all possible cost values), and p(k) the denotes the probability of the occurrence of k cost. The G k is an unweighted graph having k wired cost and it is generated from G by an appropriate threshold and X(G k) is the value of the X topological parameter calculated on G k. We assumed that k has uniform distribution, which means that the p(k) = 1/Mc, so the cost-integrated version of any X can be calculated as follows:In our case, the number of nodes is N = 149 and so the MC = 11026, we had to apply the Monte-Carlo (MC) based estimation to calculate cost integrated values (described also in [69]). Our preliminary analysis showed that in the case of 5 randomly selected subjects 200 MC-sample gave correct Eg and El values, so this sample rate was used in this study. The cost integration range was set as [0.05, 0.95].

In our special case the cost-integrated nodal degree (k, nodal efficiency (e and local efficiency (El*) were calculated for a subject as described in the Document S1 and illustrated in Panel B of Figure 2.

Software environment

For the network generation and for the cost-integrated evaluation of topological metrics we used FSL and the following in-house developed software:

BrainLOC (www.minipetct.com/brainloc) atlas database management software was used to generate 149-regional database from different atlases end evaluate averaged regional BOLD curves.

BrainNetTools parallelized utilities running on the high-performance cluster (HPC) of University of Debrecen were used to generate wavelet correlation matrices, connectivity matrices and to calculate cost-integrated topological metrics.

BrainCON (www.minipetct.com/braincon) was used for network visualization.

The components of BrainNetTool were validated by the following external software.

Brain Connectivity Toolbox (BCT, http://www.brain-connectivity-toolbox.net/), release 29/03/2012.

igraph 0.6 R-package (http://igraph.sourceforge.net/)

wavelets 0.2–7 R-package (http://cran.r-project.org/web/packages/wavelets/index.html)

NetworkAnalysis 0.3–1 R-package (http://cran.r-project.org/package=NetworkAnalysis)

Utilizing the statistical parametric network approach

The statistical parametric network (SPN) terminology was introduced by Ginestet and Simmons [77]. We used the population level mean SPN to demonstrate the subgraph of regional connectivity system. The generated mean SPN provides a method to statistically infer the mean wavelet correlation matrix of the population, and it contains z-scores of correlation coefficients:Where r− and sd(r) stand for the grand sample mean and grand sample standard deviation. The re

denotes the mean correlation coefficients of the edge e:

Correlation model, covariates and contrasts

A general linear model (GLM) based variance analysis was performed. The effects of the dimensional autistic traits (ASSQ and SRS scores) on the relevant graph theoretical and diffusion anisotropy were tested, age and gender terms were added as covariates. Random effects were modeled by using the randomly assigned subject identifier in the model. The general model was formulated according to:where Xij is the ith graph-theory based global or local topological descriptor or the regional diffusion anisotropy for the jth subject. Intercept and β terms are fixed effects (β1: explanatory variables – psychological scores, covariates: age, gender or scanner dependency), dij stands for random effects (within-person dependence) and eij is the residual error. In our study design, we used two different sets of covariates. Here we describe these models by listing the relevant dependent variables or explanatory variables, factors and covariates. The hypothesized interaction is defined by a contrast matrix in which the columns correspond to the elements of the parameter vector.

The first model assumes an inverse relationship between the psychological measurements and the predicted variable (i.e. functional topology or diffusion anisotropy). Age, gender and subject identifier are taken as covariates. The relevant parameter vector and contrast matrix are described below.We hypothesize that the intelligence profile (i.e. performance IQ and verbal IQ) interacts with regional graph theoretical parameters and also correlates with the autistic traits. Therefore we aimed to utilize a model which is adjusted for this putative effect. The intelligence profile is commonly characterized by discrepant values of performance and verbal IQ scores [38] although the direct relationship between HFA and IQ discrepancy has been questioned [39]. In our study population, verbal IQ significantly and positively correlated with performance IQ (Pearson's r = 0.55; p<0.001) while negative correlation was discovered between the verbal IQ and the ASSQ score (r = −0.393; p<0.001) and the SRS score (r = −0.343; p<0.001). Hence in the second model, we controlled the effects of the WASI verbal IQ profile:Custom scripts in IBM SPSS Statistics for Windows 20 were used to perform multiple univariate GLM tests for the 149 brain regions. The effects of the autistic trait were tested with variance statistics, for each test we provide the F values and uncorrected p-values. The contribution of each variable to the overall variance is characterized by reporting the slope (B) and the relevant significance values. The effect size is estimated by the partial eta squared value, which is calculated using the following equation:where F is the test statistic and dfh and dfe are its degrees of freedom and degrees of freedom for error. For the interpretation of effect size, we refer to the rule of thumb by Cohen, where similarly to r2 and R2 values, 0.01 denotes small, 0.06 moderate and >0.14 marks large correlation [78]. However, it is also important to note that partial eta measures tend to overestimate the effect size.

The relevant software syntax is described in the Document S2.

False discovery rate correction

Due to the exploratory design of our study and the large number of brain regions evaluated, it was necessary to account for multiple comparisons during hypothesis testing. Here we report results as if the 149 measurements (i.e. in each evaluated brain region) were independent. After obtaining the p-value vectors for the network degree, efficiency and regional white matter anisostropy, we utilized the Benjamini-Hochberg (BH) procedure for false discovery rate (FDR) control [79], with α level set at 0.05. The relevant SPSS syntax of the correction procedure is described in the Document S2.

Results

Validation of network measurements

We analyzed the population level cost-dependent distribution of Eg and El parameters, which is illustrated in Figure 3. Figure 4 demonstrates that the Eg and El of the brain networks monotonically increase by the cost with relatively low standard deviation. It means that all brain networks have high Eg and El values in the [0.34–0.5] cost-range, a conservatively defined small-world regime [80], which verifies that these networks have small-world properties [70]. We also calculated the cost-dependent efficiencies for random (red) and regular (blue) networks. In the case of random network containing 149 nodes, similarly to the brain networks, the edges were randomly generated in proportion to the given wired cost. The calculation was repeated 30 times for evaluating averaged cost-dependent Eg and El values. Efficiencies of regular networks were also calculated from 149-node networks, but in this case the required edge density was guaranteed by regularly distributed edges. It seems that the brain networks are better suited for efficiency-based small-world criteria (high E) than the regular and random network networks, i.e. at any cost level; the brain networks have higher Eg than the regular and higher El values than the random networks.

Figure 3

Distribution of global functional network properties.

Global efficiency (Eg) and average local efficiency (El). Data are depicted as histograms (left panels) and Q-Q (quantile-quantile) probability plots (right panels) in which reference lines of the normal distribution are given (n = 127).

Figure 4

Small-world properties of functional and synthetic networks.

Global (Eg) and local (El) efficiencies are depicted as a function of wired-cost for random (red), a regular (blue) and the investigated human brain networks. In latter case, the averages (black line) and the standard deviation (gray band) of efficiency values are shown. The Eg and El of brain networks monotonically increase by the cost with relatively low standard deviation. This means that all brain networks have simultaneously high Eg and El values in the [0.34–0.5] cost-range, which verify that these networks have small-world properties.

Global functional connectivity

Global brain functional connectivity was portrayed by two graph theoretical parameters (Eg, El); here we describe their observed characteristics and interactions with the phenotypic descriptors. All parameters were normally distributed in the study cohort (K-S test; asymptotic significances; Eg: p = 0.63; El: p = 0.99). The factor gender was not associated with any of the global graph theoretical parameters, age showed significant positive correlation with El (ηp 2: 0.046; F = 6.54; p = 0.012). WASI verbal IQ scores did not correlate with the connectivity parameters (controlling for age and gender). ASSQ values were found to significantly predict the averaged local efficiency parameter El when controlling for the effects of age and gender (ηp 2 = 0.046; F = 6.33; p = 0.013); the SRS Total Score alone did not show such effect. Similarly to the analysis of regional parameters, a combined model of the ASSQ and SRS, controlling for age and gender effects was utilized. This revealed small, significant association between the psychological trait and average local efficiency El (ηp 2 = 0.047; F = 5.99; p = 0.016). When controlling the model for the effects of the verbal IQ profile, the correlation with average local efficiency was slightly reduced (ηp 2 = 0.038; F = 4.39; p = 0.038). Numeric data of the described interactions are given in Table 2.

Table 2

Global graph theoretical network correlates of autistic traits, age, gender and IQ using GLM.

Dependent variables in the GLM design	Covariates in GLM (variables to control)	Association with Global efficiency Eg	Association with Average local efficiency El
Age	-	F: 0.25	F: 6.54
		p: 0.62	p: 0.012*
		ηp 2: 0.002	ηp 2: 0.046
		Bage: 6.07*10−5	Bage: 2.59*10−4
Gender	-	F: 0.55	F: 2.02
		p: 0.46	p: 0.16
		ηp 2: 0.004	ηp 2: 0.015
		BGender: 2.67*10−3	BGender: 3.9*10−3
WASI Verbal IQ	Age, Gender	F: 1.97	F: 0.68
		p: 0.16	p: 0.41
		ηp2: 0.016	ηp2: 0.0055
		BWASI: 1.94*10−4 (ηp 2: 0.016; p: 0.16)	BWASI: 8.79*10−5 (ηp 2: 0.0055; p: 0.41)
		BAge: 4.59*10−5 (ηp 2: 8.6*10−4; p: 0.75)	BAge: 2.45*10−4 (ηp 2: 0.04; p: 0.026)*
		BGender: 8.8*10−4 (ηp 2: 4.57*10−4; p: 0.81)	BGender: 0.0035 (ηp 2: 0.012; p: 0.22)
ASSQ Score	Age, Gender	F: 0.051	F: 6.33
		p: 0.82	p: 0.013 *
		ηp2: 3.8 * 10−4	ηp2: 0.046
		BASSQ: 5.24*10−5 (ηp 2: 3.8*10−4; p: 0.82)	BASSQ: −4.28*10−4 (ηp 2: 0.046; p: 0.013)*
		BAge: 6.21*10−5 (ηp 2: 0.0016; p: 0.65)	BAge: 2.45*10−3 (ηp 2: 0.044; p: 0.015)*
		BGender: 0.003 (ηp 2: 0.0037; p: 0.49)	BGender: 0.0048 (ηp 2: 0.024; p: 0.075)
SRS Total Score	Age, Gender	F: 1.09	F: 1.36
		p: 0.29	p: 0.25
		ηp2: 0.0088	ηp2: 0.011
		BSRS: −1.04*10−4 (ηp 2: 0.0088; p: 0.29)	BSRS: −8.72*10−5 (ηp 2: 0.011; p: 0.25)
		BAge: 4.92*10−5 (ηp 2: 9.2*10−3; p: 0.74)	BAge: 2.56*10−4 (ηp 2: 0.042; p: 0.022)*
		BGender: 0.0046 (ηp 2: 0.011; p: 0.25)	BGender: 0.005 (ηp 2: 0.023; p: 0.091)
ASSQ*SRS Total Score	Age, Gender	F: 0.025	F: 5.99
		p: 0.88	p: 0.016 *
		ηp2: 2.06 * 10−3	ηp2: 0.047
		BSRS: −1.23*10−4 (ηp 2: 8.6*10−3; p: 0.31)	BSRS: 1.95*10−5 (ηp 2: 4*10−4; p: 0.83)
		BASSQ: 8.46*10−5 (ηp 2: 7.02*10−3; p: 0.77)	BASSQ: −4.62*10−4 (ηp 2: 0.037; p: 0.033)*
		BAge: 5.65*10−5 (ηp 2: 1.2*10−3; p: 0.71)	BAge: 2.16*10−4 (ηp 2: 0.031; p: 0.053)
		BGender: 4.58*10−3 (ηp 2: 0.011 p: 0.25)	BGender: 0.0049 (ηp 2: 0.024; p: 0.089)
ASSQ*SRS Total Score	Age, Gender, WASI Verbal IQ	F: 0.77	F: 4.39
		p: 0.39	p: 0.038 *
		ηp2: 0.007	ηp2: 0.038
		BSRS: −9.14*10−5 (ηp 2: 0.0047; p: 0.47)	BSRS: 5.14*10−6 (ηp 2: 2.5*10−5; p: 0.96)
		BASSQ: 3.2*10−4 (ηp 2: 0.01; p: 0.28)	BASSQ: −4.26*10−4 (ηp 2: 0.03; p: 0.065)
		BAge: 3.73*10−5 (ηp 2: 5.2*10−3; p: 0.81)	BAge: 1.88*10−4 (ηp 2: 0.022; p: 0.12)
		BVIQ: 3*10−3 (ηp 2: 0.028; p: 0.078)	BVIQ: −7.34 *10−6 (ηp 2: 2.9*10−5; p: 0.96)
		BGender: 0.0023 (ηp 2: 0.0028; p: 0.58)	BGender: 4.8*10−3 (ηp 2: 0.021 p: 0.13)

GLM: General Linear Model, ASSQ: Autism Spectrum Screening Questionnaire, SRS: Social Responsiveness Scale.

: p<0.05

Regional functional connectivities and diffusion anisotropy

The regional degree and efficiency values were predominantly normally distributed in the study population (K-S tests). In the first part of the regional evaluations, we tested the joint effects of the ASSQ and SRS on graph theoretical parameters. This model was adjusted to control for the subject's age and gender. Significant inverse, moderate correlation was discovered between the autistic traits and the network degree and efficiency of the left and right posterior cingulate gyrus, although after FDR adjustment, only that of the right posterior cingulate gyrus remained statistically significant (nodal degree: ηp 2 = 0.11; F = 14.98; p = 0.00018; nodal efficiency: ηp 2 = 0.11; F = 14.43; p = 0.00024). The diffusion anisotropy of these two areas did not present significant association with the examined autistic traits. After accounting for contributions of age, gender and random effects, we found that higher scores on the SRS and ASSQ – indicating greater autistic traits – were associated with lower diffusion anisotropy in the right parahippocampal gyrus, posterior part (ηp 2 = 0.069; F = 8.81; p = 0.0036) and right temporal fusiform cortex, posterior part (ηp 2 = 0.054; F = 6.82; p = 0.01). However, we emphasize that this interaction was not proven to be significant after FDR adjustment by the BH procedure. We summarized the results of this analysis in Table 3.

Table 3

Regional neural correlates of self-reported autistic traits, GLM design 1.

Region name	Nodal Degree	Nodal Efficiency	Regional Diffusion Anisotropy
Right cingulate gyrus, posterior part	F: 14.98	F: 14.34	F: 0.27
	p: 0.00018 *** †	p: 0.00024 *** †	p: 0.61
	ηp2: 0.11	ηp2: 0.11	ηp2: 0.002
	BASSQ: −0.51 (ηp 2: 0.11; p: 1.64*10−4)***	BASSQ: −1.9*10−3 (ηp 2: 0.11; p: 1.8*10−4)**	BASSQ: −7.15*10−5 (ηp 2: 0.002; p: 0.65)
	BSRS: 0.083 (ηp 2: 0.019; p: 0.13)	BSRS: 3.5*10−4 (ηp 2: 0.022; p: 0.1)	BSRS: 2.37*10−6 (ηp 2: 1.1*10−5; p: 0.97)
	BAge: −0.018 (ηp 2: 0.001; p: 0.79)	BAge: −8.8*10−5 (ηp 2: 0.001; p: 0.74)	BAge: −9.47*10−7 (ηp 2: 0.012; p: 0.23)
	BGender: −0.41 (ηp 2: 4.33*10−4; p: 0.82)	BGender: −1.8*10−3 (ηp 2: 5.8*10−4; p: 0.79)	BGender: 4*10−3 (ηp 2: 0.034; p: 0.044)*
Left cingulate gyrus, posterior part	F: 8.05	F: 7.18	F: 0.21
	p: 0.0053 **	p: 0.0084 **	p: 0.65
	ηp2: 0.062	ηp2: 0.056	ηp2: 0.002
	BASSQ: −0.37 (ηp 2: 0.068; p: 0.004)**	BASSQ: −1.4*10−3 (ηp 2: 0.065; p: 0.004)**	BASSQ: 3.15*10−5 (ηp 2: 2.61*10−4; p: 0.86)
	BSRS: 0.072 (ηp 2: 0.016; p: 0.16)	BSRS: 3.08*10−4 (ηp 2: 0.02; p: 0.12)	BSRS: 3.78*10−5 (ηp 2: 0.002; p: 0.61)
	BAge: 0.062 (ηp 2: 0.008; p: 0.33)	BAge: 2.3*10−4 (ηp 2: 0.007; p: 0.34)	BAge: −1.11*10−4 (ηp 2: 0.012; p: 0.23)
	BGender: −0.27 (ηp 2: 2.1*10−4; p: 0.88)	BGender: −1.2*10−3 (ηp 2: 2.9*10−4; p: 0.85)	BGender: 2.2*10−3 (ηp 2: 0.007; p: 0.36)
Right parahippocampal gyrus, posterior part	F: 0.3	F: 0.34	F: 8.81
	p: 0.58	p: 0.56	p: 0.0036 **
	ηp2: 0.002	ηp2: 0.003	ηp2: 0.069
	BASSQ: 0.12 (ηp 2: 0.002; p: 0.66)	BASSQ: 6.7*10−4 (ηp 2: 1.4*10−3; p: 0.69)	BASSQ: −1.1*10−3 (ηp 2: 0.066; p: 0.005)**
	BSRS: 0.007 (ηp 2: 2.8*10−5; p: 0.95)	BSRS: 1.4*10−4 (ηp 2: 3.5*10−4; p: 0.84)	BSRS: 1.5*10−4 (ηp 2: 0.008; p: 0.34)
	BAge: 0.085 (ηp 2: 0.003; p: 0.54)	BAge: 1.52*10−6 (ηp 2: ∼0 p: 0.99)	BAge: 5.32*10−5 (ηp 2: 0.001; p: 0.78)
	BGender: 6.34 (ηp 2: 0.024; p: 0.089)	BGender: 0.051 (ηp 2: 0.04; p: 0.026)*	BGender: 3.8*10−3 (ηp 2: 0.005; p: 0.45)
Right temporal fusiform cortex, posterior part	F: 0.34	F: 0.51	F: 6.83
	p: 0.56	p: 0.48	p: 0.01 *
	ηp2: 0.003	ηp2: 0.004	ηp2: 0.054
	BASSQ: 0.14 (ηp 2: 0.003; p: 0.53)	BASSQ: 7.5*10−4 (ηp 2: 0.004; p: 0.48)	BASSQ: −6.7*10−4 (ηp 2: 0.051; p: 0.013)*
	BSRS: −0.029 (ηp 2: 0.001; p: 0.75)	BSRS: −1.11*10−4 (ηp 2: 0.001; p: 0.8)	BSRS: 8.95*10−5 (ηp 2: 0.005; p: 0.42)
	BAge: 0.084 (ηp 2: 0.005; p: 0.45)	BAge: 3.5*10−4 (ηp 2: 0.003; p: 0.53)	BAge: 9.46*10−5 (ηp 2: 0.004; p: 0.49)
	BGender: 3.7 (ηp 2: 0.013; p: 0.21)	BGender: 0.015 (ηp 2: 0.009; p: 0.29)	BGender: −1.4*10−3 (ηp 2: 0.001; p: 0.71)

We demonstrate the joint effects of Social Responsiveness Scale (SRS Total Score) and Autism Spectrum Screening Questionnaire (ASSQ), controlling for age and gender. GLM: General Linear Model.

: 0.01≤p<0.05.

: 0.001≤p<0.01.

: p<0.001.

: Significant interaction after false discovery rate correction with the Benjamini-Hochberg procedure.

In the next series of evaluations, we adjusted the GLM to account for the confounding effects of WASI verbal IQ. In this test, the putative interaction between regional graph theoretical variables and the verbal intelligence profile is controlled. The significance level of the interaction between autistic traits and right posterior cingulate gyrus degree and efficiency was mildly reduced but small correlations were found (degree: ηp 2 = 0.094; F = 11.46; p = 0.00098; efficiency: ηp 2 = 0.078; F = 9.45; p = 0.0026). After our conservative FDR adjustment procedure neither the left nor the right posterior cingulate gyrus showed significant association with the physiological markers, although the interaction effect in these regions still remained prominent when comparing to other areas. In this model, the diffusion anisotropy of the left anterior parahippocampal gyrus was found to be significantly correlating with the autistic trait (ηp 2 = 0.094; F = 11.46; p = 0.00098). Similar effect was noted for the right posterior parahippocampal gyrus (ηp 2 = 0.101; F = 12.29; p = 0.00066), left posterior temporal fusiform gyrus (ηp 2 = 0.093; F = 11.28; p = 0.0011) and the right posterior temporal fusiform gyrus (ηp 2 = 0.109; F = 13.46; p = 0.00038). The relationship of the predicted and observed values of graph measurements and diffusion anisotropy is depicted in Figure 5 and Figure 6. Supporting data are given in Table 4. Unadjusted p values for each brain region are listed in Table S2.

Figure 5

Regional graph theoretical correlates of self-reported autistic traits.

Left panels: observed vs. predicted values of the dependent variables in the general linear model analysis. Middle panels: predicted values vs. ASSQ score. Right panels: predicted values vs. SRS total scores.

Figure 6

Regional diffusion anisotropy correlates of self-reported autistic traits.

Left panels: observed vs. predicted values of the dependent variables in the general linear model analysis. Middle panels: predicted values vs. ASSQ score. Right panels: predicted values vs. SRS total scores.

Table 4

Regional neural correlates of self-reported autistic traits, controlled for verbal IQ: GLM design 2.

Region name	Nodal Degree	Nodal Efficiency	Regional Diffusion Anisotropy
Right cingulate gyrus, posterior part	F: 11.46	F: 9.45	F: 1.01
	p: 0.00098 ***	p: 0.0026 **	p: 0.32
	ηp2: 0.094	ηp2: 0.078	ηp2: 0.0022
	BASSQ: −0.49 (ηp 2: 0.102; p: 5.6*10−4)***	BASSQ: −1.8*10−3 (ηp 2: 0.092; p: 0.001)**	BASSQ: −1.2*10−4 (ηp 2: 0.0047; p: 0.47)
	BSRS: 0.081 (ηp 2: 0.017; p: 0.17)	BSRS: 3.5*10−4 (ηp 2: 0.02; p: 0.13)	BSRS: −2.7*10−5 (ηp 2:1.3*10−3; p: 0.71)
	BAge: −0.017 (ηp 2: 4.7*10−4; p: 0.82)	BAge: −6.9*10−5 (ηp 2: 5.4*10−4; p: 0.81)	BAge: −7.6*10−5 (ηp 2: 0.0069; p: 0.38)
	Bviq: 0.04 (ηp 2: 2.3*10−3; p: 0.61)	Bviq: 2.9*10−4 (ηp 2: 7.9*10−3; p: 0.35)	Bviq: −1.6*10−4 (ηp 2: 0.026; p: 0.088)
	BGender: −0.25 (ηp 2: 1.6*10−4; p: 0.89)	BGender: −1.9*10−3 (ηp 2: 6.2*10−4; p: 0.79)	BGender: 0.0053(ηp 2: 0.049; p: 0.019)*
Left cingulate gyrus, posterior part	F: 5.35	F: 3.82	F: 0.00
	p: 0.023 *	p: 0.53	p: 1.00
	ηp2: 0.046	ηp2: 0.033	ηp2: 0.0018
	BASSQ: −0.36 (ηp 2: 0.064; p: 0.007)**	BASSQ: −1.3*10−3 (ηp 2: 0.053; p: 0.014)*	BASSQ: −2.2*10−5 (ηp 2: 1.1*10−4; p: 0.91)
	BSRS: 0.095 (ηp 2: 0.025; p: 0.092)	BSRS: 3.9*10−4 (ηp 2: 0.03; p: 0.067)	BSRS: 1.4*10−5 (ηp 2: 2.8*10−4; p: 0.86)
	BAge: 0.074 (ηp 2: 0.01; p: 0.28)	BAge: 2.9*10−4 (ηp 2: 0.011; p: 0.28)	BAge: −1.5*10−4 (ηp 2: 0.021; p: 0.13)
	Bviq: 0.071 (ηp 2: 0.008; p: 0.34)	Bviq: 3.8*10−4 (ηp 2: 0.016; p: 0.19)	Bviq: −1.3*10−4 (ηp 2: 0.014; p: 0.22)
	BGender: −0.056 (ηp 2: 9*10−6; p: 0.97)	BGender: −9.6*10−4 (ηp 2: 1.7*10−4; p: 0.89)	BGender: 0.0028 (ηp 2: 0.01; p: 0.29)
Left parahippocampal gyrus, anterior part	F: 0.37	F: 0.83	F: 11.46
	p: 0.544	p: 0.36	p: 0.00098 *** †
	ηp2: 0.003	ηp2: 0.007	ηp2: 0.094
	BASSQ: 0.33 (ηp 2: 0.016; p: 0.18)	BASSQ: 1.8*10−3 (ηp 2: 0.019; p: 0.14)	BASSQ: −7.5*10−4 (ηp 2: 0.073; p: 0.004)**
	BSRS: −0.2 (ηp 2: 0.031; p: 0.062)	BSRS: −8.4*10−4 (ηp 2: 0.023; p: 0.11)	BSRS: −5.7*10−6 (ηp 2: 2.5*10−4; p: 0.96)
	BAge: 0.11 (ηp 2: 6.1*10−3; p: 0.41)	BAge: 4.7*10−4 (ηp 2: 4.7*10−3; p: 0.47)	BAge: −1.5*10−4 (ηp 2: 0.012; p: 0.25)
	Bviq: 0.075 (ηp 2: 3.2*10−3; p: 0.59)	Bviq: 5.6*10−4 (ηp 2: 5.8*10−3; p: 0.42)	Bviq: −5.6*10−4 (ηp 2: 0.12; p: 1.7*10−4)***
	BGender: 0.94 (ηp 2: 1.4*10−3; p: 0.78)	BGender: 3.3*10−3 (ηp 2: 3.6*10−4; p: 0.84)	BGender: 8.7*10−4 (ηp 2: 5.7*10−4; p: 0.8)
Right parahippocampal gyrus, posterior part	F: 0.19	F: 0.56	F: 12.29
	p: 0.66	p: 0.46	p: 0.00066 *** †
	ηp2: 0.002	ηp2: 0.005	ηp2: 0.101
	BASSQ: 0.11 (ηp 2: 1.3*10−3; p: 0.69)	BASSQ: 9.8*10−4 (ηp 2: 2.8*10−3; p: 0.58)	BASSQ: −1.2*10−3 (ηp 2: 0.088; p: 0.001)**
	BSRS: −4.3*10−4 (ηp 2: ∼0; p: 0.99)	BSRS: 1.8*10−4 (ηp 2: 5.1*10−4; p: 0.81)	BSRS: 7.1*10−5 (ηp 2: 1.7*10−3; p: 0.66)
	BAge: 0.14 (ηp 2: 7.8*10−3; p: 0.35)	BAge: 2.8*10−4 (ηp 2: 8.2*10−4; p: 0.76)	BAge: −1.2*10−4 (ηp 2: 0.003; p: 0.54)
	Bviq: 0.025 (ηp 2: 2.2*10−4; p: 0.87)	Bviq: 7.3*10−4 (ηp 2: 4.7*10−3; p: 0.47)	Bviq: −3.6*10−4 (ηp 2: 0.025; p: 0.098)
	BGender: 6.2 (ηp 2: 0.022; p: 0.11)	BGender: 0.048 (ηp 2: 0.035; p: 0.046)*	BGender: 0.0047 (ηp 2: 0.007; p: 0.37)
Left temporal fusiform cortex, posterior part	F: 0.01	F: 0.11	F: 11.28
	p: 0.93	p: 0.74	p: 0.0011 ** †
	ηp2: 6.9 * 10−5	ηp2: 9.7 * 10−4	ηp2: 0.093
	BASSQ: −0.012 (ηp 2: 3.1*10−4; p: 0.95)	BASSQ: 3.1*10−4 (ηp 2: 1.1*10−3; p: 0.73)	BASSQ: −7.4*10−4 (ηp 2: 0.073; p: 0.004)**
	BSRS: −3.6*10−3 (ηp 2: 1.7*10−5; p: 0.97)	BSRS: −4.8*10−5 (ηp 2: 1.4*10−4; p: 0.9)	BSRS: −1.8*10−6 (ηp 2: 3.1*10−6; p: 0.99)
	BAge: 0.11 (ηp 2: 0.01; p: 0.29)	BAge: 4.3*10−4 (ηp 2: 7.6*10−3; p: 0.36)	BAge: −1.1*10−4 (ηp 2: 0.006; p: 0.41)
	Bviq: −0.094 (ηp 2: 6.4*10−3; p: 0.39)	Bviq: −2.3*10−4 (ηp 2: 2*10−3; p: 0.64)	Bviq: −4.3*10−4 (ηp 2: 0.076; p: 0.003)**
	BGender: −8.2*10−3 (ηp 2: ∼0; p: 0.99)	BGender: 1.7*10−4 (ηp 2: 2*10−7; p: 0.99)	BGender: 1.3*10−3 (ηp 2: 1.2*10−3; p: 0.71)
Right temporal fusiform cortex, posterior part	F: 0.1	F: 0.32	F: 13.46
	p: 0.75	p: 0.57	p: 0.00038 *** †
	ηp2: 9.1 * 10−4	ηp2: 0.003	ηp2: 0.109
	BASSQ: 0.14 (ηp 2: 3.3*10−3; p: 0.54)	BASSQ: 8.9*10−4 (ηp 2: 5.6*10−3; p: 0.43)	BASSQ: −8.7*10−4 (ηp 2: 0.085; p: 0.002)**
	BSRS: −0.076 (ηp 2: 5.3*10−3; p: 0.45)	BSRS: −3.4*10−4 (ηp 2: 4.5*10−3; p: 0.48)	BSRS: −6.8*10−6 (ηp 2: 3.1*10−5; p: 0.95)
	BAge: 0.043 (ηp 2: 1.1*10−3; p: 0.73)	BAge: 1.8*10−4 (ηp 2: 8.8*10−4; p: 0.76)	BAge: 5.9*10−5 (ηp 2: 1.6*10−3; p: 0.67)
	Bviq: −0.13 (ηp 2: 8.3*10−3; p: 0.34)	Bviq: −4.6*10−4 (ηp 2: 4.7*10−3; p: 0.47)	Bviq: −4.9*10−4 (ηp 2: 0.083; p: 0.002)**
	BGender: 4.2 (ηp 2: 0.016; p: 0.18)	BGender: 0.016 (ηp 2: 0.01; p: 0.29)	BGender: 4*10−4 (ηp 2: 1.1*10−4; p: 0.91)

We demonstrate the joint effects of Social Responsiveness Scale (Total Score) and Autism Spectrum Screening Questionnaire, controlling for age and gender and WASI Verbal IQ ( = viq).

: 0.01≤p<0.05.

: 0.001≤p<0.01.

: p<0.001.

: Significant interaction after false discovery rate correction with the Benjamini-Hochberg procedure.

An important part of the analysis is to interpret graph theoretical results in a broader context: to relate regional impairments to the brain's overall functional network organization. To portray the sub-network to which the posterior cingulate gyrus belongs, we provide quantitative values for the strengths of first degree connections in this area in Table 5. In subjects with high ASSQ scores, the right posterior cingulate is most strongly connected to its contralateral equivalent and to anterior cingulate gyri in both hemispheres. Further strong connections were found with the precuneus and cuneus cortex, to the temporal connectivity domains of the thalamus and to the right frontal pole. A notable difference is seen in the network structure in low and high scorers of the ASSQ test (the later marking greater traits of autism). For illustrating these networks, we used the correlation coefficient threshold of 0.634 (95th percentile strongest connection), this selected the 9 regions that were most strongly interconnected to the right posterior cingulate gyrus. The same threshold in low scorers revealed 34 connected regions. As the network density is defined by the wiring cost [69], correlation coefficients may not be suitable for inter-subject comparisons – we propose to use the statistical parametric network approach. This would help tackle this problem in a way that graph edges (connections) significantly deviating from the group mean can be pinpointed. Using a correlation threshold of p = 0.005 during SPN analysis, a similar difference was noted between high and low ASSQ scorer groups. Using the SPN approach, 6 regions were connected to the right posterior cingulate in high ASSQ scorers, while the same cut-off threshold depicted 75 connections in low ASSQ scorers. This approach reveals an overall reduction of first degree connections of the right posterior cingulate, more specifically, reduced connections to the frontal lobe, superior temporal regions and the occipital lobe. We illustrate the network structure of the right posterior cingulate gyrus in Figure 7, Figure 8 and Figure 9.

Table 5

Connections of the right posterior cingulate gyrus in subjects with low and high degree of autistic traits.

Group	Connected brain regions and averaged functional connectivity values (correlation coefficient)
Subjects with the highest ASSQ Total score (n = 52, cut-off value = 2)	l-CingGyP (0.91), r-PreCC (0.82), l-CunC (0.79), r-CingGyAnt (0.79), l-PreCC (0.78), r-CunC (0.76), l-CingGyAnt (0.65), r-ThTemp (0.63), r-FrP (0.63)
Subjects with the lowest ASSQ Total score (n = 7, cut-off value = 24)	l-CingGyP (0.94), r-PreCC (0.85), l-PreCC (0.83), r-CingGyAnt (0.78), l-CunC (0.75), l-CingGyAnt (0.74), r-CunC (0.74), r-ThPreFr (0.72), r-ThTemp (0.72), l-ThPreFr (0.71), l-LingGy (0.7), r-SupFrGy (0.69), l-ThPostPar (0.69), r-LingGy (0.69), r-FrP (0.68), l-IntCalC (0.68), l-ThTemp (0.67), r-IntCalC (0.67), l-LatOccCS (0.66), r-PreCentGy (0.66), r-MiFrGy (0.66), l-CerVI (0.66), l-CerCrI (0.66), r-LatOccCS (0.66), l-ThOcc (0.65), l-LatOccCInf (0.65), l-Hyp (0.65), r-ParCingGy (0.65), r-ThPostPar (0.64), r-SupTempGyA (0.64), r-CerVI (0.64), r-SupCalcC (0.64), r-Hyp (0.64), r-PlaT (0.69)

Connection strengths are given as the wavelet correlation coefficients of the BOLD time curves. Cut-off threshold for lowest and highest ASSQ scorers were the 5th and 95th percentile values. Regions were listed in descending order of the connection strengths. Key for the abbreviations is given in Table S1.

Figure 7

Functional network of the right posterior cingulate cortex, calculated for sub-populations with the lowest and highest ASSQ scores.

Graph edges were depicted based on the strongest functional connectivity (threshold criteria for averaged networks: wavelet correlation coefficient >0.634; 95th percentile strongest connection in high scorers; threshold criterion for the statistical parametric network (SPN) method: p<0.005). The nomenclature for brain region abbreviations is given in Table S1. Raw connectivity data are provided in Table 5. First row: averaged functional network (first degree connections) in subjects with the lowest ASSQ score (n = 52; 5th percentile, cut-off threshold: 2). Second row: lowest ASSQ scorer group, SPN based determination of connectivity strengths. Third row: averaged functional network in subjects with the highest ASSQ score (n = 7; 95th percentile, cut-off threshold: 24). Fourth row: high ASSQ scorer group, connections are visualized using SPN analysis.

Figure 8

Demonstrating the functional connections of the right posterior cingulate cortex with a circular connectivity profile.

Connectivity data were calculated for sub-populations with the lowest and highest ASSQ scores. Key for the abbreviations is given in Table S1.

Figure 9

Circular representation of the functional connections of the right posterior cingulate cortex, SPN method.

Connectivity data were calculated for sub-populations with the lowest and highest ASSQ scores. Key for the abbreviations is given in Table S1.

Discussion

Global functional network correlates of autistic traits

Network-based approaches are valuable for autism research: they allow interpreting brain connectivity information as complex interactions between remote neuronal groups and therefore provide feasible means to study the effects of autism globally. In the graph theoretical model each network node is interconnected with a large number of brain regions, such connections are weighted by the degree of temporal coherence of low-frequency BOLD signals. Hence if a pathological process selectively damages functionally associated groups of neurons, the largest effect will presumably be observed in the region or subgraph having the highest density of connections to the affected areas (i.e. connector hubs). Similarly, generalized impairments in connections by pathologic processes can change the nature of information transfer globally, leaving their mark on graph theoretical descriptors such as global efficiency or degree [81]. In our study the association of regional network efficiency was only significant in the right posterior cingulate gyrus after adjusting for the false discovery rate. In contrast to this observation, reduced average local efficiency (El) in the entire network (averaged over 149 regions) was weakly associated with the autistic trait (ηp 2 = 0.047; F = 5.99; p = 0.016). This result can putatively be ascribed to the summation of impairments in distributed regions that are not significant after the relatively conservative FDR adjustment (Table S2). Recently, Dennis and colleagues provided whole-brain measures for the association of structural network structure and the presence of rs2710102 CNTNAP2 autism risk gene [82]. A significant reduction of global efficiency, small-worldness and elevated characteristic path length was discovered in carriers of the gene. Using EEG data, Barttfeld et al. showed that the functional network of subjects with autism spectrum disorder have less clustering coefficient and increased characteristic path length than controls [83].

Regional functional network correlates of autistic traits

According to our functional network analysis, the posterior cingulate gyrus was affected in both hemispheres while we revealed significantly decreased network degree and efficiency for the right posterior cingulate cortex. The cingulate region was already proved to play a role in the decreased inhibitory control in high-functioning autists [84]. Additional support comes from a study describing neuronal migrational disorder in the cingulate cortex [85], although the case report by Korkmaz and co-authors does not imply that such a migration deficit could represent common pathogenic pathway for autism spectrum disorders. A histological analysis of autistic brains by Oblak and colleagues showed that the PCC had altered cytoarchitecture with irregularly distributed neurons, poorly demarcated layers IV and V, and increased presence of white matter neurons [86]. The possible role of the posterior cingulate in the manifestation of adverse socio-emotional behaviors was postulated: reduced GABA(A) receptors were found in this region and also in the fusiform gyrus [87]. By explicating the neural correlates of moral reasoning in ASD, Schneider et al. revealed decreased activation in the amygdala while an increase was seen in the cingulate gyri [88]. When we controlled the extended functional network of the posterior cingulate gyrus, first degree connections were found with the cuneus and precuneus cortex, dorsolateral prefrontal cortex (DLPFC) and superior temporal areas (Table 5, Figure 7). This pattern coincides with the common descriptions of the task negative network ( = TNN) [89]. The TNN is known to activate during performance of emotional, social and introspective tasks, such as the theory of mind [90], social perception, emotional processing, viewing familiar faces [91] and its impairments have been frequently linked with the autism spectrum disorder. Even though there is a considerable uncertainty about the function of the posterior cingulate cortex, functional neuroimaging studies postulate it as a cortical hub and assume a role in regulating the balance between internally and externally directed cognition [92]. Intrinsically defined maps of the TNN were found to be altered in autists [93]. Kennedy and co-authors showed that the posterior cingulate cortex and precuneus are reproducibly included in TNN in autistic subjects while other components are diminished, although this effect was only significant for the medial prefrontal cortex and the left angular gyrus. The lowered temporal synchrony of such interconnected areas of the TNN, including the posterior cingulate, can result in reduced network degree or nodal efficiency which is coherent with our findings. The limbic system including the cingulate gyrus is closely tied to emotion and social behaviors and it is an almost trivial statement that the disrupted limbic circuitry could provoke some of the behavioral deficits seen in autists [94].

A review by Courchesne et al. postulates that frontal lobe abnormality is to be held responsible for the social, emotional and cognitive impairments in autism [29]. This is supported by observations of early localized enlargement of the frontal lobe in young children [95] and impaired frontal activation in tasks evaluating selective attention [96] or language skills [97]. In contrast to the commonly described frontal abnormality linked with the autistic trait, we were unable to show such a pattern. Interestingly, the temporal connectivity domain of the thalamus (as defined by the Oxford Thalamus connectivity atlas) was found to be among the first degree connections of the right posterior cingulate gyrus and hence its influence on the network efficiency can be assumed. The thalamus plays a role in cognitive processes by relaying cortico-striato-thalamo-cortical loops: three of such circuits originate from the prefrontal (dorsolateral prefrontal and medial orbitofrontal cortex) and limbic (anterior cingulate cortex) cortical areas and pass through the mediodorsal thalamic nucleus [98]. This model was further stratified by in vivo neuroimaging studies mapping such neuronal pathways [99] and revealing correlations between the thalamic representation of such circuits and executive functions of human subjects [100]. In the present study we have included the thalamus and various other subcortical territories in the construction of the functional network. Although we did not reveal localized disruption in network topology of the thalamus, it occurs that the autism spectrum disorder as a distributed disease could affect the thalamus [101], [102] and the basal ganglia as well [103].

People with autism spectrum disorder exhibit heterogeneous social-communication characteristics, a phenotype that presumably extends into the neurotypical population [104]. So far, we reported our findings in comparison to the studies using clinical (low or high functioning autistic) populations. A different design of investigations was necessary to find evidence about such autistic traits in normally developed subjects. In a study by Hagen et al. structural and functional impairment of superior temporal sulcus (pSTS) was implicated to correlate with self-reported traits autism spectrum characteristics [105]. Similarly to our study, Di Martino et al. described the neural correlates of the autistic traits in neurotypical populations and revealed reduced cingulo-insular functional connectivity in resting state functional imaging data [37]. Significant negative relationship of this specific connection was found with the SRS score which is partially coherent with our findings, however, currently we do not report the effects on edges ( = individual connections) but provide descriptors of network topology.

White matter impairments linked with autistic traits

Measuring functional network topology and regional diffusion anisotropy in the same autistic subject group can provide complementary information on a more pervasive structural deficit that is characterized by the disruption of white matter integrity, altered connections or conversely, increased anisotropy could indirectly indicate “more ordered” axons compared to the control group [106]. In both hemispheric temporal fusiform gyri, we have revealed diffusion anisotropy deficits. A mixed pattern of elevated and reduced FA values was reported by Cheng and colleagues [19]. Although individuals affected by autistic spectrum disorders do not usually show severe behavioral consequences of prosopagnosia, abnormal response to face stimuli was reported [107], [108], reduced functional connectivity from the fusiform face area (FFA) and its associated networks was shown [109], but conversely, the role of FFA impairment in autism was doubted by others [110]. Moreover, neuropathological evaluations discovered microscopic structural deficits in this territory, smaller neuronal volumes and densities were found [111]. Besides the diffusion abnormalities of the fusiform gyrus, we showed the bilateral association of the FA values of the parahippocampal gyrus and the self-reported autistic traits. Among distributed pattern of cortical regions, the cortical thickness of the left parahippocampal gyrus was found to be a predictor of autism spectrum disorders [112]. In contrast to this observation, Ke and colleagues linked the reduction of the right parahippocampal gray matter volume with HFA [113]. The significance of the posterior cingulate cortex in autism spectrum disorders was also noted by Uddin et al. who used multivariate pattern analysis to find the most discriminative features between normal subjects and ASD. The largest discriminative power was achieved by using gray matter thickness measurements in the posterior cingulate cortex, the parahippocampal gyri and the anterior temporal lobe [114].

Limitations and future directions

A major limitation in our study protocol is that neither the ASSQ nor the SRS provides explicit diagnostic data on the autism spectrum disorder including HFA. Furthermore, the diagnostic border between Asperger's syndrome and HFA is not well-defined; it is acknowledged that the key to the diagnosis is the presence of early language developmental impairment. We did not have such historic information about the adult subjects involved. In conclusion, our results can only be interpreted as the effects of the self-reported autistic trait on the functional brain connectome and diffusion anisotropy patterns, observed in normally functioning adults. We raised the specificity of our study by employing a conservative FDR adjustment method, although many of our measurements are highly dependent from each other; for instance, hemispheric equivalent regions are highly synchronized (r>0.9). Therefore the decision threshold in terms of accepted p-values could be higher and impairments in a larger set of regions would be linked to the autistic traits. Furthermore, by building and analyzing larger and clinically controlled datasets such as the ABIDE ( = Autism Brain Imaging Data Exchange), it will become possible to explicate the neural correlates of this puzzling mental condition while more plausible biomarkers could be developed for the diagnostics and prognostics of ASD.

Conclusions

Recent large-scale neuroimaging data-sharing initiatives such as the INDI provided valuable normative functional imaging resources for discovery research. Such complex datasets allow the observation of the relationship between brain structure and dimensional psychological scales in neurotypicals. One such challenging task is exploring traits with relatively low prevalence like autism. We conclude that high functioning autistic adults carry peculiar brain abnormalities that are linked with regionally disrupted functional networks and altered regional diffusion anisotropy of the brain, especially in the limbic cortex while reduced diffusion anisotropy is present in temporal fusiform and parahippocampal gyri. A general tendency was observed towards lower averaged local efficiency in subjects with higher scores of self-reported ASSQ and SRS. In contrast to studies revealing frontal and temporal abnormalities in autism, we did not reveal such effects in neurotypical adults with different degrees of autistic traits.

Description of the custom brain region system used for creating network nodes and regions of the diffusion anisotropy analysis.

(DOCX)

Click here for additional data file.

Regional neural correlates of self-reported autistic traits, unadjusted F and p values for all brain regions.

(XLS)

Click here for additional data file.

Pseudo code of cost integration method used for calculating the global and regional graph theoretical network descriptors.

(PDF)

Click here for additional data file.

Scripts used in SPSS to perform GLM analysis and false discovery rate adjustment.

(DOCX)

Click here for additional data file.