Keywords

1 Introduction

Research has indicated that the brain’s Functional Connectivity (FC) does not remain stationary over time, and that it has multiple states distinguished by unique connectivity patterns [1, 3, 7, 18, 19]. Efforts are being made to characterize how neurological ailments affect brain connectivity and its dynamic properties [10].

Several studies have shown that Autism Spectrum Disorder (ASD) alters brain connectivity [8, 9]. Accuracy has been improved by taking the dynamic property of FC into account in identifying subjects with ASD [12]. This suggests that brain dynamics possess noteworthy information about the disorder. In this work, our three primary goals are:

  1. 1.

    To model brain dynamics in autistic children.

  2. 2.

    To examine if there is any deviation in brain dynamics of autistic subjects in comparison to typically developing ones of similar demographics.

  3. 3.

    To investigate how autism affects the temporal properties of the brain if any anomalies caused by autism are discovered.

The triple network model consists of three Intrinsic Connectivity Networks (ICNs), namely Default Mode Network (DMN), Salience Network (SN) and Central Executive Network (CEN) [15]. It has been proven to be of primary importance in understanding higher-order cognitive function and dysfunction [11]. Aberrations in the engagement and disengagement of these three core networks have been shown to play a significant role in neurological disorders such as schizophrenia, depression, anxiety, dementia and autism [5, 6, 11, 16]. Supported by previous research claims indicating that neurological disorders have an impact on the functioning of the triple network model, we focus our work on SN, CEN, and DMN.

Despite the existence of multiple approaches to model brain dynamics, majority of them have not focused on computing useful statistical measurements of the temporal properties of these time-varying states. Moreover, many of the approaches rely on the usage of sliding windows, which might lead to erroneous assumptions of the number of brain states. A recent framework leveraging the power of Hidden Markov Models (HMMs) with Variational Bayesian Inference has been proposed by Ryali et al. [15] that addresses the aforementioned setbacks by accurately estimating several useful parameters. It also eliminates the need for arbitrary assumptions on the number of brain states. In their work, they have used the VB-HMM framework to show how interactions between DMN, SN and CEN mature with age by conducting a study on adult and child cohorts [15].

As a natural extension, we thought of characterizing how autism influences the brain dynamics in children aged between 6.5 and 14. It has been found that children diagnosed with ASD exhibit a deviation in brain connectivity in comparison to typically developing ones of similar demographics. In this paper, several insightful properties of the brain states have been quantified and discussed. Such detailed analysis on the effects of autism is believed to be contributory to better comprehending of the disorder.

2 Methods

2.1 Variational Bayes Hidden Markov Models (VB-HMM)

VB-HMM is a framework that employs Bayesian Inference approximation to handle large scale problems that would otherwise be intractable using traditional ML-HMMs [13], and was proposed by Ryali et al. [15]. Notations and equations used in this section to describe VB-HMM are adapted from their paper [15].

Notations: Matrices are represented using uppercase, while scalars and vectors are represented using lowercase. Let T be the number of time samples, and S be the number of subjects. \(y_{t}^{s}\) is an M dimensional time sample for subject s at time t, where M is number of ROIs (Region of Interests), then \(Y = \left\{ \left\{ y_{t}^{s}\right\} _{t=1}^{T}\right\} _{s=1}^{S}\) represents the observed voxel time series data. Let \(Z = \left\{ \left\{ z_{t}^{s}\right\} _{t=1}^{T}\right\} _{s=1}^{S}\) be the underlying hidden discrete states, where \(z_{t}^{s}\) is the state label for subject s at time t.

Z is a first order Markov chain with transition matrix A and initial distribution \(\pi \). The probability of the observation \(y_{t}\) given its state is assumed to be a multivariate Gaussian distribution with mean \(\mu _{k}\) and covariance \(\varSigma _{k}\).

Let the unknown parameters of the HMM model be \(\varPhi =\{\pi , A, \varTheta \}\), where \(\varTheta =\left\{ \mu _{k}, \varSigma _{k}\right\} _{k=1}^{K}\). Then, \(p(Y, Z, \varPhi )\) will be the joint probability distribution of the observations, hidden states, and the required HMM parameters. The objective is to model the HMM parameters \(\varPhi \) accurately. The traditional maximum likelihood approach may result in inaccurate characterization of the parameters since it requires that the number of hidden states to be specified a priori. With the help of Bayesian inference, high precision approximations can be attained.

Let us consider \(p(Z, \varPhi |Y)\) to be the true posterior and \(q(Z, \varPhi | Y)\) to be any arbitrary probability distribution. The Kullback-Leibler (KL) divergence is used to measure of how different q is from the true posterior distribution:

$$\begin{aligned} KL(q|p)=-\int dZd\varPhi q(Z, \varPhi |Y)\log \frac{p(Z, \varPhi |Y)}{q(Z, \varPhi |Y)} \end{aligned}$$
(1)

The Kullback-Leibler measure is greater than or equal to zero at all times and becomes zero only when both of the distributions in comparison are identical. Now, the log of marginal distribution of observations Y can be expressed as

$$\begin{aligned} \log P(Y)=F(q)+KL(q||p) \end{aligned}$$
(2)

where F(q) is the negative free energy, which is computed as

$$\begin{aligned} F(q)=\int dZd\varPhi q(Z, \varPhi |Y)\log \frac{p(Y, Z, \varPhi )}{q(Z, \varPhi |Y)} \end{aligned}$$
(3)

A strict lower bound is defined by F(q) on \(\log P(Y)\), and the objective of Bayesian approximation is to model q such that the lower bound F(q) is maximized. The posterior distribution is estimated using Expectation-Maximization (EM) algorithm similar to the Baum-Welch algorithm. Further details on the workings of VB-HMM can be found in Ryali’s paper [15].

Viterbi decoding algorithm is used for obtaining optimal hidden state sequence. For each distinct state, mean lifetime and occupancy rate are calculated from the optimal state sequence. The transition probabilities are estimated by the VB-HMM algorithm during the maximization step.

2.2 Flow of VB-HMM Analysis

The steps used in conducting the study are listed below in a sequential manner:

  1. 1.

    Build an HMM model for the observed fMRI time series data using the VB-HMM algorithm (\(\varPhi =\{\pi , A, \varTheta \}\) is estimated).

  2. 2.

    Use Viterbi decoding for estimating the underlying hidden state sequence with help of the parameters learned.

  3. 3.

    Compute mean lifetimes and occupancy rates of the dynamic brain states.

  4. 4.

    Compute FC matrices for each state from the estimated covariance matrices.

  5. 5.

    Apply Louvain Community Detection algorithm on the FC matrices to detect their underlying community structures.

  6. 6.

    Combine states with similar underlying community structure and compute temporal properties of the resultant states.

Finally, we analyze the impacts of ASD by comparing the brain dynamics of autistic and typically developing cohorts.

3 Data Description

Preprocessed Data. Preprocessed fMRI time series data has been obtained from the “ABIDE preprocessed” project [4]. Data preprocessed with Configurable Pipeline for the Analysis of Connectomes (CPAC) using the Automated Anatomical Labeling (AAL) parcellation has been used in conducting the study. A detailed report on the preprocessing steps can be obtained from the ABIDE website (http://preprocessed-connectomes-project.org/abide/).

For every subject, 176 frames were acquired with a repetition time (TR) of 2s and the scan duration used in this study is approximately six minutes (352s).

Regions of Interest (ROIs). While each of the three ICNs consist of several regions, few representative regions have been selected based on their importance in higher-order cognition. The regions used to represent SN are Right Insula (Insula_R) and Right Midcingulate Area (Cingulum_Mid_R). For CEN, the regions used are Right Middle Frontal Gyrus (Frontal_Mid_R) and Right Inferior Parietal Lobule (Parietal_Inf_R). And for DMN, Right and Left Precuneus (Precuneus_L and Precuneus_R), and Right middle frontal gyrus-orbital part (Frontal_Med_Orb_R) have been selected. These regions have been identified as crucial for representing their respective networks by a previously published study [17].

Description of Cohorts. To study the impacts of ASD in children, 20 typically developing individuals and 20 individuals clinically classified as autistic were selected. These two cohorts were matched in age, gender and IQ demographics. ADOS (Autism Diagnostic Observation Schedule) total scores for ASD subjects ranged from 8 to 19.

Cohort

Subjects

Age (Mean ± Std, Range)

FIQ Score (Mean ± Std, Range)

NYU ASD

20, Male

10.41 ± 1.68, 6.5–14

108.15 ± 8.74, 95–130

NYU TD

20, Male

9.94 ± 1.91, 6.5–14

111.09 ± 7.19, 95–130

Fig. 1.
figure 1

Temporal properties of dynamic brain states in ASD & TD Cohorts: (a, e) Time evolution of dynamic states in each subject. (b, f) Occupancy Rates (c, g) Mean Lifetimes (d, h) Transition Probabilities of states with occupancies above chance level

Fig. 2.
figure 2

(a) Dynamic FCs and (b) DFNs of ASD Cohort. States with occupancies above chance level are shown here. In (a), for ROI names of corresponding ROI numbers, see Fig. 3. In (b), edges of the same color belong to the same community.

4 Experiments and Results

To address the sensitivity of HMMs to initialization, the experiment is run 100 times with random initializations and the model with maximum lower bound F(q) is selected. The FC matrices are computed from estimated covariance matrices of each state and the interdependence of every pair of ROIs is represented by Pearson Correlation values. To investigate the properties of underlying connectivity patterns of each brain state, the Louvain Community Detection algorithm [2] present in the Brain Connectivity Toolbox [14] has been employed and this revealed interactions among the ICNs, and also between the regions within them. The BrainNet Viewer tool was used for visualizing the communities discovered [20].

Fig. 3.
figure 3

Segregated DFN (DFN-S) and Cross-Connected DFN (DFN-C)

Brain states possessing similar connectivity patterns have been merged into two distinct dynamic functional network configurations. States in which SN, CEN and DMN formed fully segregated communities are grouped together and referred to as segregated DFNs (DFN-S) and the remaining DFNs that possess cross-network interactions are referred to as cross-connected DFNs (DFN-C) (see Fig. 3).

Fig. 4.
figure 4

(a) Dynamic FCs and (b) DFNs of TD Cohort. States with occupancies above chance level are shown here. In (a), for ROI names of corresponding ROI numbers, see Fig. 3. In (b), edges of the same color belong to the same community.

4.1 Temporal Properties of Dynamic Brain States and Their Underlying Functional Network Configurations

For allowing the VB-HMM algorithm to discover number of brain states directly from the fMRI time-series data without making arbitrary assumptions, we initialize the number of hidden states to a large upper-bound number. As a result, all of the excess states would be assigned zero occupancy leaving us with only relevant states. In this work, we set an upper-bound of 25 states (leading to an initial chance level of 4% for each state) and only seven non-zero occupancy states were found in both of the cohorts.

Dynamic Brain States. In both ASD and TD cohorts, four states with an occupancy rate above the initial chance level were found while remaining states had negligible occupancy. We can observe that the top two states in both of the cohorts have similar occupancy rates (close to 30%). The third and fourth states of each cohort had dissimilar occupancy rates, suggesting a difference in patterns of brain connectivity between ASD and TD cohorts. The mean lifetime of the top four states in ASD Cohort ranged from 24 to 30 s, while in TD Cohort it ranged from 10 to 42 s indicating a varying level of persistence in brain states. High values along the principal diagonal of the markov transition matrices indicate the stability of brain states (see Fig. 1).

Fig. 5.
figure 5

Temporal properties of DFN-S & DFN-C in ASD & TD Cohorts: (a, d) Time evolution of the two configurations in each subject. (b, e) Occupancy Rates & Mean Lifetimes (c, f) Transition Probabilities

Dynamic Network Configurations. Deeper insights about the brain dynamics are revealed by studying the underlying functional connectivity configurations. It is interesting to note that the first two states in both ASD and TD cohorts have the same underlying connectivity pattern (see Figs. 2 and 4), suggesting that these states are universally present across all of the subjects unaffected by the disease. Moreover, these states formed communities within the ICNs and had no interactions across them. DFN-S configuration (corresponding to states 1, 2 and 5 in ASD cohort and 1, 2 and 4 in TD cohort) in ASD cohort had significantly less occupancy in comparison to TD cohort. It is also observed that DFN-S has higher mean lifetime in healthy individuals, demonstrating higher persistence and stability of DFN-S in typically developing subjects (see Fig. 5).

5 Discussion

Recent studies have proven that SN, CEN and DMN formed segregated communities only intermittently [15]. It has also been shown that children’s resting brain activity spent less time in segregated functional network configurations compared to adults [15].

In addition to these earlier observations, our findings now indicate that children with ASD have lesser occupancy and shorter mean lifetime of DFN-S, exhibiting a clear deviation in brain dynamics as a result of autism. Narrowing the focus to a set of regions with well-documented roles in higher-order cognition enabled us to highlight the existence of impairment in the working memory and cognitive control. All three ICNs of the triple network model have been investigated extensively for their involvement in behavioral attributes, information processing and cognitive functioning. SN plays a crucial role in identifying and responding to salient external stimuli and internal events, whereas the CEN is responsible for high-level cognitive functions such as planning, decision making, and the control of attention and working memory [11]. The Default Mode Network forms an integrated system for self-related cognitive activity, including autobiographical, self-monitoring and social functions. Dynamic interactions among these ICNs govern attention transition and access to domain-general and domain-specific cognitive resources [11]. Impaired functioning of these networks can be held responsible for difficulties in social interaction, communication, behavior and sensory sensitivities, all of which are the most common traits of autism.

Investigating the dynamic changes in brain connectivity due to the presence of ASD at a subject level is beneficial in understanding the nature of the breakdown in function of the subject being studied. Autism refers to a wide spectrum of mental and behavioral conditions, often many of these being unique to the subject. The subject level time evolution of the brain states and their underlying configurations generated in this study can be used to identify atypical characteristics in the functioning of these networks in the corresponding subject.

6 Conclusion

In summary, we have modeled the effects of autism on the temporal dynamics of the brain. Patterns and trends that seemed unusual when compared to typically developing subjects have been identified and the underlying network configurations associated with these aberrations have been discussed.