Functional connectivity is abnormal in autism, but the nature of these abnormalities remains elusive. Different studies, mostly using functional magnetic resonance imaging, have found increased, decreased, or even mixed pattern functional connectivity abnormalities in autism, but no unifying framework has emerged to date. We measured functional connectivity in individuals with autism and in controls using magnetoencephalography, which allowed us to resolve both the directionality (feedforward versus feedback) and spatial scale (local or long-range) of functional connectivity. Specifically, we measured the cortical response and functional connectivity during a passive 25-Hz vibrotactile stimulation in the somatosensory cortex of 20 typically developing individuals and 15 individuals with autism, all males and right-handed, aged 8–18, and the mu-rhythm during resting state in a subset of these participants (12 per group, same age range). Two major significant group differences emerged in the response to the vibrotactile stimulus. First, the 50-Hz phase locking component of the cortical response, generated locally in the primary (S1) and secondary (S2) somatosensory cortex, was reduced in the autism group (P < 0.003, corrected). Second, feedforward functional connectivity between S1 and S2 was increased in the autism group (P < 0.004, corrected). During resting state, there was no group difference in the mu-α rhythm. In contrast, the mu-β rhythm, which has been associated with feedback connectivity, was significantly reduced in the autism group (P < 0.04, corrected). Furthermore, the strength of the mu-β was correlated to the relative strength of 50 Hz component of the response to the vibrotactile stimulus (r = 0.78, P < 0.00005), indicating a shared aetiology for these seemingly unrelated abnormalities. These magnetoencephalography-derived measures were correlated with two different behavioural sensory processing scores (P < 0.01 and P < 0.02 for the autism group, P < 0.01 and P < 0.0001 for the typical group), with autism severity (P < 0.03), and with diagnosis (89% accuracy). A biophysically realistic computational model using data driven feedforward and feedback parameters replicated the magnetoencephalography data faithfully. The direct observation of both abnormally increased and abnormally decreased functional connectivity in autism occurring simultaneously in different functional connectivity streams, offers a potential unifying framework for the unexplained discrepancies in current findings. Given that cortical feedback, whether local or long-range, is intrinsically non-linear, while cortical feedforward is generally linear relative to the stimulus, the present results suggest decreased non-linearity alongside an increased veridical component of the cortical response in autism.
Autism spectrum disorder (ASD) is diagnosed behaviourally, and has a complex genetic basis (Berg and Geschwind, 2012; Jeste and Geschwind, 2014) and no clear identifying biomarkers. Yet, identifying objective ASD biomarkers is absolutely crucial in the quest to better understand the aetiology of ASD. One emergent direction for the identification of objective biomarkers is based on functional connectivity in the brains of individuals with ASD. Abundant evidence suggests that ASD is associated with functional abnormalities in cortical processing, including abnormalities in functional connectivity. However, there are many seemingly disparate and often contradictory findings in the field, which show functional over-connectivity, functional under-connectivity, and even normal functional connectivity in ASD (Kana et al., 2011; Muller et al., 2011; Wass, 2011; Vissers et al., 2012). While there are many potential explanations for these disparities, ranging from methodological concerns (Jones et al., 2010b; Van Dijk et al., 2012; Nair et al., 2014) to the inherent heterogeneity of ASD (Betancur, 2010; Pelphrey et al., 2011; Di Martino et al., 2014), an often overlooked possible source is the diversity in subtypes of functional connectivity itself. Functional connectivity can be mediated, altered and constrained by different cortical rhythms (Lopes da Silva, 1991; Fries, 2005), by interactions between these rhythms, by changes in anatomical pathways, by the nature of the mediating neurons (e.g. inhibitory or excitatory), by the directionality of the connectivity (feedforward versus feedback), etc. Delineating the nature and mechanisms of the diverse functional connectivity abnormalities observed in ASD is essential for gaining insight into the fundamental aetiology and neural mechanisms of ASD.
We recently found, using graph theory, that network efficiency during resting state is increased in ASD in the gamma frequency band, but reduced in the beta frequency band (Kitzbichler et al., 2015). This finding raises the intriguing possibility that these two distinct cortical rhythms might mediate abnormal functional connectivity in ASD, but in opposite directions. Specifically, these findings suggest that functional connectivity in ASD might be abnormally increased when mediated by the gamma band, but abnormally decreased when mediated by the beta band.
This hypothesis is particularly significant in light of the fact that the gamma band has been associated with feedforward connectivity, while the beta band has been associated with feedback connectivity (Wang, 2010). Feedforward connectivity is defined as connections arriving into layer 4, the granular layer in the cortex, and progressing upwards along the cortical hierarchy. Feedforward connectivity is equivalent to bottom-up connectivity, but it is defined independently of context (e.g. emotional load, attention state, memory, or learning). Feedback connectivity refers to inputs arriving into the supragranular layers of the cortex and usually progresses downwards or in parallel across the cortical hierarchy. Feedback connectivity is equivalent to top-down connectivity, but like feedforward connectivity it is defined independently of context. If long-range functional connectivity in ASD were decreased in the feedback direction, but increased in the feedforward direction, that might explain some of the apparent discrepancies in the field, as well as clarify the aetiology of functional connectivity abnormalities in ASD.
The division of functional connectivity into feedforward and feedback components applies only to long-range connections. It has long been hypothesized that local functional connectivity is also abnormal in ASD, and in particular that it is increased in individuals with the disorder (Belmonte et al., 2004; Minshew and Williams, 2007). The term ‘local’ has, however, been used loosely in the ASD field. Henceforth we use the term ‘local’ to mean connectivity within a single functionally defined cortical region. This spatial scale aligns well with the spatial resolution achievable with non-invasive methods. More importantly, it is consistent with the well-known and anatomically traceable functional connections that spread laterally across the non-granular layers within a functionally defined cortical region. These local connections are usually referred to as recurrent connections (Lamme and Roelfsema, 2000).
In a previous study (Khan et al., 2013) we found evidence of reduced local functional connectivity in ASD that was localized to the fusiform face area and manifested as reduced interactions across different cortical rhythms, during a face processing task. This finding contradicts the hypothesis that local functional connectivity is increased in ASD, as least while performing this paradigm, and perhaps during other paradigms as well. If reduced local functional connectivity in a region driven by a functionally appropriate task is indeed characteristic of ASD, as we speculated, it would suggest that we should be able to observe reduced local functional connectivity in other cortical areas and using different experimental paradigms.
To address these open questions surrounding both long-range and local functional connectivity abnormalities in ASD, we sought a study that would allow us to distinguish between long-range feedforward and feedback interactions, while also distinguishing both from local functional connectivity mediated via recurrent local interactions. Results from a previous study of behavioural responses to vibrotactile stimuli (rapid vibrations applied to the skin) on fingertips in ASD suggested that studying the cortical responses to vibrotactile stimulation in ASD would achieve this goal. In healthy individuals, vibrotactile stimulations on the fingertips dramatically increase the behaviourally measured temporal separation threshold needed for distinguishing between two consecutive tactile pulses, but no such effect is observed in ASD individuals (Tommerdahl et al., 2008). These results suggest that the vibrotactile stimulus elicited a very different cortical response in ASD individuals. The interpretation proposed was that local functional connectivity is reduced in the somatosensory cortex in ASD.
The cortical response to a vibrotactile stimulus has been studied extensively in neurotypical subjects and localizes well to early somatosensory cortex (Jamali and Ross, 2012). It consists of phase locking at the stimulus frequency, along with additional components of the response that manifest as phase locking at integer multiples of the stimulus frequency generated locally in the cortex (Langdon et al., 2011). The focus on the somatosensory cortex made this paradigm compelling for testing feedforward and feedback connections in ASD, as it is possible to measure directed (i.e. feedforward and feedback) connectivity between S1 and S2, given sufficient temporal resolution.
In parallel to the vibrotactile paradigm, we also examined long-range feedback and feedforward connectivity during the resting state by measuring the sensorimotor mu rhythm. The mu rhythm is observable at rest in the sensorimotor cortex and consists of two components: mu-α, in the alpha frequency range (8–13 Hz), and mu-β, in the beta frequency range (20–24 Hz) (Pineda, 2005). Recent computational modelling work suggests that the mu-β component of the somatosensory mu rhythm is generated via balancing between long-range feedforward and feedback connections (Jones et al., 2009), making it relevant to the current study of balance between these two streams of functional connectivity in ASD.
Behavioural and diagnostic data also support the importance of investigating cortical neurophysiology in somatosensory cortex in ASD. ASD has long been associated with tactile processing abnormalities (Tommerdahl et al., 2008; Wiggins et al., 2009; Green et al., 2013). Yet the neural signatures of such abnormalities, and the neurophysiological mechanisms that may underlie them, have not been mapped to date.
Based on our previous findings, we hypothesized that we would find reduced local functional connectivity (reduced recurrence) and increased feedforward connectivity in ASD individuals during a passive vibrotactile paradigm. We also hypothesized that long-range feedback connectivity, assessed indirectly via the mu rhythm, would be reduced in ASD. Lastly, we hypothesized that the magnitude of functional connectivity abnormalities in ASD would correlate with behaviourally measured tactile abnormalities and with ASD severity, thus linking them directly to behavioural manifestations of ASD.
To test our hypotheses, we used magnetoencephalography (MEG), a technique with millisecond time resolution and good spatial resolution, to study the cortical responses to 25 Hz pneumatic vibrotactile stimulus applied passively to the fingertips in 20 typically developing and 15 ASD individuals, aged 8–18, and the mu rhythm during resting state in a subset (12 per group) of these participants.
Materials and methods
Participants were 15 males diagnosed with ASD and 20 age-matched typically developing males, ages 8–18 (11.6 mean age). ASD participants had a prior clinically verified ASD diagnosis, met a cut-off of >15 on the Social Communication Questionnaire, Lifetime Version, and were assessed with either Module 3 (n = 3) or 4 (n = 12) of the Autism Diagnostic Observation Schedule (ADOS, Lord, 1999), administered by trained research personnel who had previously established interrater reliability. Individuals with autism-related medical conditions (e.g. Fragile-X syndrome, tuberous sclerosis) and other known risk factors (e.g. premature birth) were excluded from the study. All typically developing participants were below threshold on the Social Communication Questionnaire and were confirmed to be free of any neurological or psychiatric conditions, and of substance use for the past 6 months, via parent and self-reports. The ASD and typically developing groups did not differ in verbal or non-verbal IQ, as measured with the Kaufman Brief Intelligence Test – II (Kaufman and Kaufman, 2004). Handedness information was collected using the Dean Questionnaire (Piro, 1998). Only right-handed participants were included in the study. Additional details on the participants are provided in Supplementary Tables 1 and 2.
Vibrotactile stimulation in the MEG consisted of pulses applied to the index and middle right fingers at 25 Hz using a custom made pneumatic tactile stimulator with latex tactor tips, based on a published design (Briggs et al., 2004). The duration of each stimulus train was 500 ms with an interstimulus interval of 3 s and a 500 ms jitter. The stimuli were presented while participants were watching a movie. Participants were instructed to not pay attention to the stimulation and not move their hands. Hands were kept still using an armrest, and a blanket positioned over the arm. The sequence of stimuli was presented using the psychophysics toolbox (www.psychtoolbox.org).
The resting state paradigm consisted of a red fixation cross at the centre of the screen, presented for 5 min continuously, while participants were seated and instructed to fixate on the cross.
Structural MRI data acquisition and processing
T1-weighted high-resolution magnetization-prepared rapid gradient echo (MPRAGE) structural images were acquired using a 3.0 T Siemens Trio whole body magnetic resonance scanner (Siemens Medical Systems) and a 32-channel head coil. The in-plane resolution was 1 × 1 mm2, slice thickness 1.3 mm with no gaps, and a repetition time/inversion time/echo time/flip angle 2530 ms/1100 ms/3.39 ms/7°. Cortical reconstructions and parcellations for each subject were generated using FreeSurfer (Dale et al., 1999; Fischl et al., 1999a). After correcting for topological defects, cortical surfaces were triangulated with dense meshes with ∼130 000 vertices in each hemisphere. For visualization, the surfaces were inflated, thereby exposing the sulci (Dale et al., 1999).
MEG data acquisition and preprocessing
MEG data were acquired inside a magnetically shielded room (IMEDCO) using a whole-head VectorView MEG system (Elekta-Neuromag), comprised of 306 sensors arranged in 102 triplets of two orthogonal planar gradiometers and one magnetometer. The signals were filtered between 0.1 and 200 Hz and sampled at 600 Hz. The position and orientation of the head with respect to the MEG sensor array was recorded continuously with help of four head position indicator (Zaidel et al., 2009) coils (Uutela et al., 2001). To allow co-registration of the MEG and MRI data, the locations of three fiduciary points (nasion and auricular points) that define a head-based coordinate system, a set of points from the head surface, and the sites of the four head position indicator coils were digitized using a Fastrak digitizer (Polhemus) integrated with the Vectorview system. The ECG and electrooculography signals were recorded simultaneously to identify epochs containing heartbeats as well as vertical and horizontal eye movement and blink artefacts. During data acquisition, online averages were computed from artefact-free trials to monitor data quality in real time. All off-line analysis was based on the saved raw data. In addition, 5 min of data from the room void of a subject were recorded before each experimental session for noise estimation purposes.
Cleaning and motion correction
The data were spatially filtered using the signal space separation method (Elekta-Neuromag Maxfilter software) to suppress noise generated by sources outside the brain (Taulu et al., 2004; Taulu and Simola, 2006). Signal space separation also corrects for head motion between and within runs (Taulu et al., 2004). Cardiac and ocular artefacts were removed by signal space projection (Gramfort et al., 2013). The data were low-pass filtered at 145 Hz to remove the head position indicator coil excitation signals.
Epoching: tactile paradigm
The data were epoched into single trials lasting 2.5 s, from 1000 ms prior to stimulus onset to 1500 ms after it. A total of 100 trials were collected. Epochs were rejected if the peak-to-peak amplitude during the epoch exceeded 1000 fT and 3000 fT/cm in any of the magnetometer and gradiometer channels, respectively. This resulted in the loss of 2–20 trials per participant. To maintain a constant signal-to-noise ratio across conditions and participants, the number of trials per condition per participant was fixed at 80, the minimum number of accepted trials that we had for each condition and participant. For participants that had more than 80 good trials, we selected 80 trials randomly from the available trials.
Epoching: resting state
The data were epoched into single trials lasting non-overlapping 4 s. A total of 75 trials were collected. Epochs were rejected if the peak-to-peak amplitude during the epoch exceeded 1000 fT and 3000 fT/cm in any of the magnetometer and gradiometer channels, respectively. This resulted in the loss of 5–25 trials per participant. As explained before, we selected, randomly when needed, 50 good trials per participant.
For both tactile and resting state MEG recordings, there were no group differences in overall quality of the data, and the number of good (unrejected) trials per condition was similar between groups and across conditions. For each participant, the same set of trials was used for all analyses.
Mapping MEG data onto cortical space
The cortical source space consisted of 10 242 dipoles per hemisphere, corresponding to a spacing of ∼3 mm between adjacent source locations. The forward solution was computed using a single-compartment boundary-element model (Hämäläinen and Sarvas, 1989). The individual inner skull surface triangulations for this model were generated with the watershed algorithm in FreeSurfer. The current distribution was estimated using the minimum-norm estimate by fixing the source orientation to be perpendicular to the cortex (Gramfort et al., 2014). The noise covariance matrix was estimated from data acquired in the absence of a subject prior to each session. We used depth weighting to reduce the bias of the minimum norm estimates towards superficial currents (Lin et al., 2006).
Intersubject cortical surface registration for group analysis
A morphing map to optimally align the cortical surface of each participant to an average cortical representation (FsAverage in FreeSurfer) was computed in FreeSurfer (Fischl et al., 1999a).
Delineating the primary and secondary somatosensory areas
To identify S1 and S2, we first mapped averaged MEG sensor space data to the individual cortical manifold through MNE. The activation during the steady state period (250 ms to 550 ms) was first averaged on the cortical manifold in time. Regions of interest were then delineated individually by setting a global threshold on the cortical activations at Z-score > 8 (Dale et al., 2000; Gramfort et al., 2014). The Z-score represented the statistical distance of cortical activation during the steady state period from the empty room MEG recordings collected before the subject arrival. The average sizes of the delineated S1 and S2 were 321.2 mm2 ± 25.7 mm2 (38.7 ± 2.9 dipoles) and 158.9 mm2 ± 9.6 mm2 (22.3 ± 1.4 dipoles), respectively. The same procedures were used for both S1 and S2. Since steady state stimulation produces very strong cortical responses, averaging across time during this period allowed us to reliably identify S1 and S2 in each participant individually. The MNI coordinates of each region of interest are shown in Supplementary Table 3. Lastly, for group analysis, we mapped S1 and S2 to FsAverage space by using morphing maps, computed earlier (Fig. 1B).
The exact location of S1 and S2 is not uniform across subjects. This is why we presented S1 and S2 as contours (Fig. 1). Supplementary Fig. 1 shows examples of S1 and S2 locations from four participants (two typically developing, two ASD), chosen randomly, on both the uninflated and the inflated brain.
Point spread and source localization
Due to the spatial resolution of ∼1 cm, the MNE for a focal source can extend across sulcal walls separated by only a few millimetres (Liu et al., 2002; Lin et al., 2006; Hauk et al., 2011). In our case, the problem is somewhat worsened by the intersubject variability in S1 localization on the central sulcus, as illustrated in Supplementary Fig. 1. This issue is prominent in Fig. 4, where the spread of activation around S1 appears quite large. However, this spread is due mostly to variability in S1 location across participants, and to ‘bleeding’ across sulci walls, as illustrated in Supplementary Figs. 1 and 2, on the uninflated brain. In a visualization on the inflated cortex the spread appears accentuated because adjacent sulcal walls appear to be spatially separated. Due to these point spread effects, we cannot ascertain with absolute certainty that the 25 Hz and 50 Hz phase locking is localized exclusively to S1 or S2. The extent of the spread is much higher in S1 than in S2 due to the larger activation observed there, but as is obvious from Supplementary Fig. 2, even if the actual source of the activity is spread beyond S1, the spread is limited, and the activation is still local, and more importantly, continuous, i.e. showing no signs of recruitment of long-range functional connectivity.
Intertrial phase locking (PL) is a method to quantify phase synchrony across multiple trials. To compute phase locking, we convolved the epoched time series with a dictionary of complex Morlet wavelets (each spanning seven cycles). We then normalized the resulting complex coefficients by dividing by their absolute magnitude and averaging the unit-norm phasors across trials for each time-frequency bin. We then took their absolute value so that each number ranged between 0 and 1, with 0 representing a uniform distribution of phase angles and 1 representing perfectly synchronized phase angles, across trials (Tallon-Baudry et al., 1996; Makeig et al., 2002). Mathematically phase locking is defined as:
Where ØK represent instantaneous phase resulting from convolution of the trial with the complex Morlet wavelet, and N is the numbers of trials.
Normalized phase locking
To compute Z-PL, we compared each phase locking value to a set of surrogate null distributions, to correct for statistical biases proportional to the number of epochs. This approach is non-parametric, and makes no a priori assumptions besides the independence across the trials in the experimental data. The independence across trials was motivated by the fact that there was an average 3 s time interval between trials, and anticipation effects were eliminated because our experimental paradigm had a 500 ms jitter in Stimulus-Onset Asynchrony. Z-PL was computed as follows: each trial was first circularly shifted by a random lag [, where T = period (1/f) in samples] and phase locking was computed on the shifted epoched data. This process was repeated 500 times. Z-PL was then computed by subtracting the mean and dividing by the standard deviation of the null distributions from the actual phase locking values (Supplementary Fig. 4).
Phase locking and Z-PL in S1/S2
Epoched time series of every vertex inside S1 and S2 were decomposed in complex time frequency domain by convolving them with seven-cycle Morlet Wavelets on a linear frequency scale (10–130 Hz). Phase locking was then computed as described earlier and normalized by the null distribution, with each realization of the null distribution computed by first circularly shifting each trial by a random lag , and then taking the maximum value across times for each frequency. To form a null distribution, this process was repeated 500 times. By taking the maximum value across time we also corrected for the multiple comparisons using maximum statistics. The process was repeated for every vertex inside S1 and S2 and the results were then averaged. All the multiple comparison corrections are done in both time and frequency using tempo-frequential maximum cluster statistics (Maris and Oostenveld, 2007). More specifically, we used 1000 permutation, and Wilcoxon Rank-Sum test as the test statistics.
Phase locking and Z-PL on the full cortical surface
We estimated phase locking at 25 Hz and 50 Hz, as previously described, at each cortical vertex, during the steady state component of the response (250 ms to 550 ms). We then normalized phase locking with respect to the null distribution as above for . Finally, to correct for the multiple comparison as phase locking was computed across every vertex on the cortex, we again used maximum statistics, and took the maximum value across the vertices. This process was repeated 500 times, forming a null distribution. Z-PL was then computed as above. Z-PL thus represents a Z-score of the phase locking corrected for multiple comparison in space.
Granger causality between S1 and S2
We used a non-parametric approach to compute Granger causality between S1 and S2 (Dhamala et al., 2008a). We first computed time frequency decomposition of both the S1 and S2 epoched time series using seven-cycle Morlet Wavelets on a linear frequency scale (20–60 Hz). Spectral matrix factorization and subsequently Granger causality scores were computed as described by Dhamala et al. (2008b). Again, all the multiple comparison corrections are done in both time and frequency using tempo-frequential maximum cluster statistics (Maris and Oostenveld, 2007), with 1000 permutation, and Wilcoxon Rank-Sum test as the test statistics.
Resting state power
Spectral power (2–55 Hz) of resting state MEG cortical data was estimated using multi-taper Fourier transform with dpss taper at spectral smoothing of 2 Hz (Thompson, 1982). We then identified individual mu-α () and mu-β power by averaging spectral power over motor and somatosensory cortices and identifying relative peaks in the spectrum. Both peaks were identified independently by visual inspection of the individual Fourier spectrum.
Local functional connectivity index and the mu-β index
Selection of behavioural measures
We used a hypothesis-driven approach to select the behavioural measures for the correlations. The ADOS scores were correlated with both the local functional connectivity index (LFCi) and the Granger causality score for the ASD group only. For sensory processing, we tested only the two scores from the Sensory Processing Questionnaire that were relevant for somatosensory processing: the touch score and the multisensory processing score. Each of these was correlated with both the LFCi and the Granger causality score, for each group separately.
All correlation coefficients and the corresponding P-values were computed using Pearson correlation. Correlations resulting in significant P-values were then tested using Robust Correlation (Pernet et al., 2012), which strictly checks for false positive correlations using bootstrap resampling and six correlation tests (bootstrap Pearson correlation, bootstrap Spearman correlation, bootstrap Bend correlation, bootstrap Pearson skipped correlation and bootstrap Spearman skipped correlation). Significant correlations were further tested for survival of multiple comparison correction by controlling for family-wise error rate using maximum statistics through permutation testing (Groppe et al., 2011). Only correlations that remained significant are shown in Figs 6 and 7A–C. Lastly, we also performed hierarchical regression analyses to test for group main effects. This analysis further confirmed the significance of each tested correlation, and also further confirmed that the LFCi and Granger causality values were indeed independent.
Linear discriminant analysis
Subjects were randomly partitioned into 10 equal-size subsamples. Of the 10 subsamples, a single subsample was retained as the validation data for testing the model, and the remaining nine subsamples were used as training data. The cross-validation process was then repeated 10 times, with each of the subsamples used once as the validation data. The standard error associated was <0.05% for all values (sensitivity, specificity, accuracy), confirming the validity of the classifier.
Statistical analyses of MEG data
Our statistical analyses were based on cluster-based statistics and a non-parametric method (Maris and Oostenveld, 2007) that also corrects for multiple comparisons. We used 1000 permutations and the test statistics used were Wilcoxon Rank Sum test. For time-frequency plots (Figs 2 and 3), correction was done in both time and frequency. For the whole cortex analyses (Figs 4 and 5), the corrections for multiple comparisons was performed across the vertices on the cortical manifold.
The model is described in detail elsewhere (Jones et al., 2009). To adapt the model to the vibrotactile paradigm, we first made sure to keep all parameters in their neurophysiological range, as detailed in the original publication. To simulate typically developing data, the 25 Hz component was reproduced with a 25 Hz stochastic feedforward burst of two spikes with an interstimulus interval of 7 ms into the lower (feedforward) layer of the model. Feedback was simulated as a stochastic burst of three spikes at 25 Hz into the upper (feedback) layer in the model, delayed by 14 ms with respect to the feedforward drive. To best replicate the averaged response, a small temporal variability (σ2 = 2 ms) was introduced to both the feedforward and feedback drives. For modelling the ASD results, we used a data-driven approach to determine the per cent increase and decrease of the feedforward and feedback components, respectively. To that end, we first calculated the (EDR) of the 25 Hz and 50 Hz components (EDR25 and EDR50) in the S1 region of interest as EDR = (ETD – EASD) / ETD, where ETD (EASD) was the average energy density of the typically developing (ASD) signal at 300 ms between 23 and 27 Hz and 48 and 52 Hz, respectively. EDR50 and EDR25 were found to be +94% and −21%, respectively. Starting from the simulated typically developing signal, feedback conductance gain was decreased from its initial value by 32% to get a simulated EDR50 of 94%, and feedforward conductance gain was then increased to 50% of its initial (typically developing) value, to get a simulated energy density ratio value of −21% (Supplementary Fig. 6). The simulations ran for 1 s and the first 550 ms were discarded to allow for initialization of the model’s internal dynamics.
Evoked response to tactile vibrations
We began by analysing the evoked responses to the vibrotactile stimulus and using these responses, recorded in the MEG sensors, to map the data onto cortical space. As expected, the evoked responses to the stimulus (Fig. 1A) localized to the contralateral (left) S1 and S2 (Fig. 1B). Although the evoked responses for the ASD and typically developing groups did not fully overlap, group differences were not statistically significant.
Local functional connectivity via phase locking within S1 and S2
The evoked response is the result of averaging over multiple frequency bands. Therefore, it may mask more subtle group differences, for instance differences that are frequency band specific. We therefore proceeded to investigate the cortical response using phase locking, and a region of interest-based analysis, focusing only on the cortical response inside S1 and S2 (Fig. 1B). For a steady-state vibrotactile stimulus, we expected to see strong phase locking at the frequency of the stimulus itself (25 Hz) as well as integer multiples of the stimulus frequency (Langdon et al., 2011). To test this hypothesis, we computed phase locking inside S1 and S2 for every frequency between 8 Hz and 125 Hz.
In S1, as expected, we observed phase locking at 25 Hz in both groups. The typically developing group also exhibited the expected strong phase locking at 50 Hz. In contrast, the ASD group showed weak phase locking at 50 Hz at the onset but nearly no phase locking subsequently, during the steady-state period of the response (Fig. 2A and B). The onset period of the 50 Hz phase locking was marginally stronger in the ASD group than in the typically developing group (P < 0.04 corrected). The steady-state component of the phase locking at 50 Hz was significantly stronger in the typically developing group (P < 0.003 corrected).
Similar results were observed in S2. Both groups showed phase locking at 25 Hz, and the typically developing group also showed significant phase locking at 50 Hz (Fig. 2C and D). Unlike in S1, phase locking at 25 Hz during the steady state period of the response was significantly stronger in the ASD group (P < 0.046 corrected), while the phase locking at 50 Hz during the steady state period was stronger in the typically developing group (P < 0.006 corrected). In S2 there were no significant group differences in the onset component of the response.
There was also evidence of phase locking at 75 Hz and 100 Hz in the typically developing group. Those signals were very weak, however, so we did not analyse them further.
S1-S2 long range functional connectivity
We next tested whether increased feedforward connectivity from S1 to S2 could account for the increased strength of the 25 Hz phase locking in the ASD group during the steady state period. To that end, we computed the directed (i.e. both feedforward and feedback) coherence between S1 and S2 using non-parametric Granger causality. We found significantly increased feedforward S1 to S2 connectivity (P < 0.004 corrected) at 25 Hz in the ASD group relative to the typically developing group (Fig. 3). No other significant group differences were observed for this analysis, including no differences in S2 to S1 feedback, as anticipated, since this no-task paradigm is not expected to drive feedback connectivity.
Spatial distribution of the phase locking response
To extend our region of interest-based analysis spatially and verify the localization of the observed responses, we computed phase locking at 25 Hz and 50 Hz at every vertex on the cortical surface, collapsing over the time dimension and considering only the steady state period of the response. Phase locking at 25 Hz and 50 Hz in the typically developing group was only statistically significant in areas contiguous with S1 and S2 (Fig. 4A and B, top), and no spatially or functionally independent cortical areas showed significant phase locking. In the ASD group, phase locking was statistically significant at 25 Hz in S1 and S2 and their contiguous surroundings, but there was no significant phase locking anywhere on the cortex at 50 Hz (Fig. 4A and B, bottom). Thus, two contiguous areas surrounding each functionally derived region of interest, one for S1 and one for S2, showed significantly reduced 50 Hz phase locking in the ASD group (P < 0.001 and P < 0.04 corrected, respectively).
Quantification of local functional connectivity during vibrotactile stimulation
It has been shown that the higher frequency component of the response, the 50 Hz phase locking in this case, is generated locally in the cortex (Langdon et al., 2011). To quantify the extent to which local functional connectivity, i.e. the 50 Hz component of the response, was reduced in ASD, we defined the vibrotactile LFCi as the strength of phase locking at 50 Hz normalized by the total strength of phase locking at 50 Hz plus 25 Hz. We computed the LFCi for every participant and for every vertex on the left cortical hemisphere. As expected from the phase locking results, the LFCi was significantly lower in the ASD group in the contoured areas shown in Fig. 4C, spanning S1 and S2 (P < 0.001, P < 0.04 corrected, respectively).
Functional connectivity and the mu-alpha and mu-beta resting state rhythms
This difference between the groups in the LFCi measure was so stark that it led us to search for additional phenomena that demonstrate frequency doubling, i.e. two related cortical responses at integer multiples of one another. The mu-β rhythm presents a particularly intriguing instance of cortical frequency doubling because it is generated during resting state, has been modelled as the result of a balance between feedforward and feedback inputs (Jones et al., 2009), and is known to be not just a harmonic double frequency component of mu-α but instead far more complex in nature (Pineda, 2005). In parallel to the evoked vibrotactile data, we obtained resting state data from a subset of the same participants (12 typically developing, 12 ASD) and analysed the mu-α and mu-β rhythms in these data. As expected, the general alpha rhythm was bilaterally strong in all the participants, and mu-α was prominent in multiple cortical sensorimotor regions (Fig. 5A). We also found, as expected, strong mu-β in sensorimotor areas in the typically developing group. However, the mu-β component was significantly weaker (P < 0.04, corrected) in the ASD group in the sensorimotor cortex (Fig. 5B).
Relationship between mu-beta and the response to vibrotactile stimulation
To probe whether mu-β and the 50 Hz phase locking component of the response to the vibrotactile stimulus may share underlying neural mechanisms, we defined the mu-β index, mu-βi, to mirror our definition of LFCi, as the strength of the mu-β component of the mu-rhythm normalized by the total strength of mu-α plus mu-β . Because the vibrotactile response localized only to the left hemisphere, we focused on the left mu-βi. As shown in Fig. 5C, mu-βi was significantly reduced in the ASD group in the sensorimotor cortex, including the area that spanned S1 (P < 0.04 corrected).
We then correlated the individual values for the 24 participants for whom both measurements were available. We found that the correlation between LFCi and mu-βi was 0.78 (P < 0.00005) when both groups were combined, and 0.72 (P < 0.008) and 0.74 (P < 0.005) within the typically developing and ASD groups, respectively (Fig. 6).
Correlation with behavioural measures
The neurophysiologically-derived LFCi was negatively correlated with the severity of ASD as measured on the ADOS (Fig. 7A). The LFCi was also positively correlated with the touch perception score (Fig. 7B), while the strength of S1 to S2 connectivity—as measured using Granger causality—was correlated with the multisensory processing score (Fig. 7C), both derived from the Sensory Profile Questionnaire. Finally, a linear discriminant analysis classifier based on the neurophysiological data had 89% accuracy, 87% sensitivity, and 90% specificity (Fig. 7D and Supplementary Fig. 3) for blindly separating the groups.
No age effects
Despite the wide age range of 8 to 18, age was not a significant factor in any of our analyses.
Modelling the origins of the 50 Hz component of the cortical response
To further test our hypothesis, a computational model that allowed the strengths of feedforward and feedback connectivity to be altered independently was examined. The model, which is biophysically realistic, shows that a reduced feedback-to-feedforward ratio eliminates the mu-β component of the cortical response. The details of the model are in the original and follow-up publications (Jones et al., 2009, 2010a), and our replication of the mu-α and mu-β results are shown in Supplementary Fig. 5.
To simulate our vibrotactile results, we began by altering the feedforward input to the model to a 25 Hz 500 ms signal, thus mimicking our vibrotactile stimulation. Feedback was simulated as 25 Hz arriving into the upper layers of the model with a delay and slightly increased strength relative to the feedforward input. We then compared model output and the typically developing group MEG data on two parameters, evoked responses and power. These parameters resulted in a reasonable reproduction of the evoked responses in the typically developing group (blue traces in Fig. 8D versus Fig. 8A) and a relatively faithful reproduction of both the 25 Hz and 50 Hz steady-state power of the MEG data obtained from the typically developing group (Fig. 8E versus Fig. 8B). We next tested whether the model could reproduce the results observed in the ASD group by altering the feedforward-to-feedback ratio, where the altered ratio estimates were derived from the MEG data of power in the ASD versus typically developing groups. This indeed resulted in a reasonably replicated evoked responses (orange traces in Fig. 8D versus Fig. 8A) and a relatively faithful reproduction of both the 25 Hz and 50 Hz steady-state power of the MEG data obtained from the ASD group (Fig. 8F versus Fig. 8C).
In this study, we found direct evidence of increased long-range feedforward functional connectivity between S1 and S2 in ASD, using passive vibrotactile stimulation. We also found reduced local functional connectivity in both S1 and S2 in ASD, using the same paradigm. In parallel, we found reduced mu-β, which is likely mediated via long-range feedback functional connectivity, in somatosensory cortex in ASD during resting state. All of the neurophysiological measures in this study were correlated with behaviourally measured traits, ASD severity, and diagnosis.
Increased feedforward connectivity in ASD
The ASD group showed a much stronger 25 Hz (i.e. stimulus frequency) phase locking in S2, alongside a non-significant trend towards increased 25 Hz phase locking during steady state in S1. These results suggest a stronger veridical representation of the stimulus in the ASD group, as the 25 Hz phase locking faithfully reproduced the stimulus. Our Granger causality analysis of the data showed greatly increased S1 to S2 feedforward connectivity in the ASD group, which likely explains why the 25 Hz phase locking component of the response was stronger in S2 than in S1 in ASD. Although no other studies to date directly measured feedforward functional connectivity in ASD, this interpretation of increased feedforward connectivity in ASD is supported by many indirect studies that show increased bottom-up processing tendencies in ASD (Neumann et al., 2006; Cook et al., 2012; Amso et al., 2014). Evidence of increased local connectivity during resting state (Berman et al., 2015) might also be the result of increased feedforward activity during this non-task state.
Reduced local functional connectivity in ASD
The most striking finding of this study was the near complete absence of 50 Hz phase locking in the ASD group during the steady state component of the response to the 25 Hz vibrotactile stimulus. To interpret this finding, it is necessary to first understand the source of this component of the response. We considered several possible scenarios. An obvious question is whether this signal is a result of an artefact or simple aliasing. In that case, however, it would appear equally in both groups, so this hypothesis was dismissed. We also verified that the 50 Hz response was not phase locked to the 25 Hz response (Supplementary Fig. 7). It is therefore not simply a harmonic component, but instead a more complex process, as in the case of mu-β (Pineda, 2005) and as suggested by previous studies of this stimulus (Langdon et al., 2011). The fact that the envelopes of the 25 Hz and 50 Hz components were not identical in time, as can easily be observed in Fig. 2, further confirms that the two components of the response are not simple harmonic multiples. Another scenario in which a ‘doubled’ response to a stimulus may occur is when neurons respond to both the ‘on’ and the ‘off’ components of the stimulus, as occurs in the magnocelullar neurons in the visual system and biphasic firing (Lee et al., 1979). That hypothesis can also be ruled out, for two primary reasons. First, we also observed faint traces of responses at 75 Hz and 100 Hz, which are incompatible with this hypothesis. Perhaps most germane, we also observed this frequency doubling (the mu-β rhythm) in the typically developing group during resting state, in the absence of any stimulation, and the lack of such doubling in the ASD group.
In agreement with previous studies of these frequency doubling harmonic-like responses in the cortex (Pineda, 2005; Langdon et al., 2011), we therefore concluded that the double frequency phase locking is generated locally via recurrent connectivity in the cortex, as indicated by the source localization of the response shown in Fig. 4 (see also Supplementary Fig. 2).
The 50-Hz response to vibrotactile stimulation versus mu-β
The highly significant positive correlation we observed between the strength of the 50 Hz component of the vibrotactile response and the strength of the mu-β rhythm across individual participants suggests that the two phenomena are related. The ease with which we were able to adapt the computational model for mu-β proposed by Jones et al. (2009) to our own findings was striking. In both cases, reduced feedback in the model ‘breaks’ the doubled frequency component of the cortical signal, just as was observed in the ASD group.
It is important to stress, however, that this computational model does not discriminate between possible different origins of what is defined as the feedback component in the model. In theory, what is referred to as feedback in the model could be generated via long-range feedback, hierarchically parallel inputs, or local functional connectivity (recurrent connections). The common theme among all of these scenarios is that the inputs would arrive at the upper layers in the cortex. The results from source localization (Fig. 4 and Supplementary Fig. 2), showing the response is localized around S1/S2, suggest it is more likely that the 50 Hz component of the vibrotactile response is driven by local recurrent connections than by long-range feedback connectivity. In contrast, given the widespread cortical distribution of the mu-β component, it is more likely to be generated at least in part via long-range functional connectivity.
The convergence and correlation between the 50 Hz component of the vibrotactile response and the mu-β component during resting state likely stem from the fact that both phenomena require additional processing in the supragranular layers of the cortex, whether via local functional connectivity (i.e. recurrent connections) or long-range feedback. These findings suggest that cortical processing in the supragranular layers is abnormal in ASD.
Studies of functional connectivity using functional MRI
The vast majority of functional connectivity studies in ASD to date were carried out using functional MRI, a technique whose temporal resolution is limited by the slow time course of the haemodynamic response. This coarse temporal resolution means that feedforward and feedback connections mediated by oscillations at different frequency bands are necessarily averaged together in functional MRI analyses. If our hypothesis, that feedforward connections are increased in ASD while feedback connections are reduced, is indeed true, then such averaging could result in reduced, increased, or seemingly unchanged averaged connectivity in ASD, depending on whether the connections being measured are dominated by feedforward or feedback components. The paradigms being studied, along with the analysis methods, would determine which processing stream prevails. Thus, our results might explain some of the inconsistencies in functional MRI studies of functional connectivity in ASD.
Models of ASD
An extensive body of evidence, spanning both psychological and neurophysiological data, suggests increased bottom-up and reduced top-down processing in ASD (Cook et al., 2012; Takarae et al., 2014). Previously, however, no specific neural mechanism had been proposed for this hypothesis. This theory is consistent with our findings of reduced feedback and increased feedforward processing in ASD, as feedforward pathways mediate bottom-up processing while feedback pathways are the primary mediators of top-down processing. Our findings also support the enhanced perceptual abilities model of ASD (Mottron et al., 2006), which suggests enhanced abilities in individuals with ASD in some tasks relying on perceptual abilities. The abilities described in the model are likely to rely on veridical feedforward processing (Mottron et al., 2013). Lastly, our findings are consistent with the reduced signal-to-noise hypothesis of ASD (Rubenstein and Merzenich, 2003). Reduced local recurrence and reduced feedback would both diminish the signal and allow increased noise from spontaneous fluctuations in brain activity that would be left unregulated, while the increased feedforward signal prevails. Such a scenario is consistent with recent findings (Edgar et al., 2015), including findings of increased trial-to-trial variability in the cortical response, as recently observed (Dinstein et al., 2012), which is also a measure of signal-to-noise ratio in the brain, in ASD.
Somatosensory abnormalities in ASD
ASD has long been associated with sensory abnormalities in general and tactile abnormalities in particular. Our findings align well with other studies finding abnormalities even in low-level cortical areas in ASD (Bertone et al., 2005; Dakin and Frith, 2005; Tommerdahl et al., 2008; Roberts et al., 2011; Dinstein et al., 2012). Together, this body of evidence suggests that some core features of ASD might result from the downstream consequences of abnormal processing of basic sensory information during development. In this study, the neurophysiological measures of the cortical response were strongly correlated with behavioural scores. Specifically, reduced local functional connectivity in S1, the primary somatosensory area, was correlated with a more severe ADOS score and with increased tactile processing abnormalities. In parallel, increased feedforward connectivity from S1 to S2 was correlated with an increase in multisensory processing abnormalities. This latter result is probably due to the fact that S2 is not only a somatosensory area, but is in close proximity to sensory association areas and has an important role in multisensory processing (Zhou and Fuster, 2000). The combination of both MEG measures, i.e. abnormalities in both local functional connectivity and long-range feedforward connectivity, was sufficient to blindly classify participants carrying an ASD diagnosis with 89% accuracy. Given that the neurophysiological abnormalities documented here are highly correlated with the ASD phenotype, our results offer a potential underlying neural mechanism for the behavioural findings of abnormal processing of tactile stimuli in ASD individuals.
The reduced non-linearity/increased veridical representation hypothesis of ASD
It is well known that both long-range feedback and local functional connectivity operate in intrinsically non-linear ways on inputs arriving to the cortex. The observations in the typically developing group of highly non-linear components in the response to a 25 Hz oscillatory input signal, presented as phase locking at 50 Hz or a strong mu-β, which has complex non-linear origins, are therefore to be expected.
The striking and novel aspect of our results is that such basic processes, stemming from the intrinsically and fundamentally non-linear processing properties of cortical function, appear severely disrupted in the ASD group. The lack of doubled frequency components in the cortical response to a vibrotactile stimulus in the ASD group, combined with the reduced mu-β, suggest a significant and substantial reduction of all inherently non-linear aspects of cortical processing in ASD. Coupled with the augmentation of the feedforward component of the cortical response, which is inherently linear (i.e. veridical), these findings suggest that ASD is characterized by reduced non-linearity in cortical processing, which can also be formulated as an abnormally increased veridical representation of the external world.
It is worthwhile noting that feedforward and feedback projections follow different developmental time trajectories, with feedback projections developing later than feedforward projections (Batardiere et al., 2002). As ASD is a developmental disorder, it is possible that due to these different developmental time courses, feedback projections would be more impacted in ASD than feedforward projections, with feedforward projections then compensating by becoming abnormally enhanced. This hypothesis would also be interesting to test in the context of minimally verbal autism relative to high functioning autism, a possibility given the simple and passive nature of the paradigm.
On the level of synaptic mechanisms, there is ample evidence that neuronal inhibition, which is central to the non-linear components of cortical processing, is abnormally reduced in ASD (Buxbaum et al., 2002; Casanova et al., 2003; Fatemi et al., 2009; Oblak et al., 2009; Rojas et al., 2014). That said, breakdowns in other synaptic mechanisms could definitely contribute to these observed abnormalities. Invasive studies in animal models would be required to discriminate between the effects of different potential contributors to the observed reduced non-linearity in ASD at the receptor and layer specific levels. The observed increase in veridical processing (i.e. feedforward connectivity) is likely also driven by reduced neuronal inhibition, or reduced gain control, of the input to the granular layer. Here too, invasive studies would be required to verify this hypothesis.
In this study we investigated long-range feedforward and feedback functional connectivity, alongside local functional connectivity, in the somatosensory cortex in typically developing and ASD individuals. We documented decreased local functional connectivity and increased long-range feedforward connectivity in ASD. We also found indirect evidence for reduced long-range feedback connectivity in ASD.
In an earlier study (Khan et al., 2013) we documented reduced local functional connectivity in ASD by measuring phase-amplitude cross frequency coupling in the fusiform face gyrus during viewing of faces. Thus we have now mapped local functional connectivity during a task in ASD in two distinct brain regions (fusiform face area previously and S1/S2 here) using two completely different analytic approaches (interactions across different cortical rhythms previously, higher order phase locking here) and two distinct paradigms (active and visual previously, passive and tactile here). In both of these otherwise independent studies, we found strong evidence of reduced recurrent interactions, i.e. reduced local functional connectivity, in ASD. In combination, our findings suggest that reduced local functional connectivity in an area driven by its functionally relevant task is likely characteristic of ASD.
Our long-range functional connectivity results are 2-fold. The mu-β rhythm results indicate that long-range feedback connectivity is likely reduced in ASD, whereas our feedforward connectivity results from S1 to S2 clearly show an increase in such connectivity in the ASD group. Both of these findings are in-line with our previous results on network efficiency in ASD (Kitzbichler et al., 2015). Along with many studies documenting what is most likely reduced feedback connectivity in ASD (Just et al., 2004; Kashino et al., 2005; Dinstein et al., 2011; Kenet et al., 2012), our findings suggest that increased long-range feedforward connectivity and reduced long-range feedback connectivity are likely characteristic of ASD.
Based on these and previous findings, we propose that a unifying feature of cortical processing in ASD is a reduction in the inherent non-linearity of the cortical response, which stems from a combination of increased weight of the veridical feedforward functional connectivity components and reduced weights of the non-linear long-range feedback and local functional connectivity components. It is important to note that even minor changes in the feedforward/feedback ratios in ASD could significantly alter perception and cognitive processing, due to the inherent non-linearity in the effects of feedback and local recurrent connections. Such altered ratios would likely be at the core of the observed sensory processing abnormalities in ASD, as well as at the core of many or most of the observed behavioural abnormalities associated with the disorder.
We would like to also thank Drs Stephanie Jones and Oleg Favorov for illuminating discussions about these findings.
This work was supported by grants from the Nancy Lurie Marks Family Foundation (T.K., S.K.), Autism Speaks (T.K.), The Simons Foundation (SFARI 239395, T.K.), The National Institute of Child Health and Development (R01HD073254, T.K.), The National Centre for Research Resources (P41EB015896, M.S.H.), National Institute for Biomedical Imaging and Bioengineering (5R01EB009048, M.S.H.), and the Cognitive Rhythms Collaborative: A Discovery Network (NFS 1042134, M.S.H.). Article needs to be deposited in PMC.
Supplementary material is available at Brain online.
Autism Diagnostic Observation Schedule
autism spectrum disorder
local functional connectivity index