Functional connectome fingerprinting: Identifying individuals and predicting cognitive functions via autoencoder.

Functional network connectivity has been widely acknowledged to characterize brain functions, which can be regarded as "brain fingerprinting" to identify an individual from a pool of subjects. Both common and unique information has been shown to exist in the connectomes across individuals. However, very little is known about whether and how this information can be used to predict the individual variability of the brain. In this paper, we propose to enhance the uniqueness of individual connectome based on an autoencoder network. Specifically, we hypothesize that the common neural activities shared across individuals may reduce the individual identification. By removing contributions from shared activities, inter-subject variability can be enhanced. Our experimental results on HCP data show that the refined connectomes obtained by utilizing autoencoder with sparse dictionary learning can distinguish an individual from the remaining participants with high accuracy (up to 99.5% for the rest-rest pair). Furthermore, high-level cognitive behaviors (e.g., fluid intelligence, executive function, and language comprehension) can also be better predicted with the obtained refined connectomes. We also find that high-order association cortices contribute more to both individual discrimination and behavior prediction. In summary, our proposed framework provides a promising way to leverage functional connectivity networks for cognition and behavior study, in addition to a better understanding of brain functions.

INTRODUCTION

Functional magnetic resonance imaging (fMRI) allows for noninvasive interrogation of brain functions based on the blood‐oxygenation‐level‐dependent‐signal (BOLD signal) (Biswal, Zerrin Yetkin, Haughton, & Hyde, 1995; Damoiseaux et al., 2006; Greicius, Krasnow, Reiss, & Menon, 2003). To comprehensively map the networks in the brain for cognitive and behavioral study, the functional connectome has been widely used by the fMRI community. Intriguingly, a functional connectome based on functional connectivity (FC) extracted from the fMRI time series provides a promising tool to investigate individual differences in human cognitive and behavioral performance from a network perspective. Recently, numerous studies have reported individual variability in functional connectivity. For instance, Airan et al. provided an in‐depth study to evaluate the degree of influence that standard fMRI acquisition and analysis schemes have on individual subject variability (Airan et al., 2016). This variability is assumed to be associated with both genetic and environmental factors, and thereby influencing neural development. Meanwhile, such variability also may partially affect individual cognition and behavior (Baldassarre et al., 2012; Gerraty, Davidow, Wimmer, Kahn, & Shohamy, 2014).

More importantly, the functional connectome of the human brain constitutes individualized patterns that enable us to identify one from a pool of individuals just like a fingerprint (Finn et al., 2015). Specifically, Finn et al. demonstrated that such connectivity profiles could be used to distinguish individuals among adult participants across rest/task modalities. In their work, they showed that the discriminative subnetworks of individuals contributed most to the prediction of fluid intelligence scores (Finn et al., 2015; Xiao, Stephen, Wilson, Calhoun, & Wang, 2019). Kaufmann et al. reported that the functional profile developed into a stable, individual wiring pattern during adolescence, and they demonstrated that reduced mental health induced a delay and an overall reduction of such wiring (Kaufmann et al., 2017). It has been pointed out that the individual differences of the functional connectome were stable across months or even years (Horien, Shen, Scheinost, & Constable, 2019; Jalbrzikowski et al., 2020). The studies mentioned above used a standard procedure to extract the functional profiles, but overlooked the influence from the group‐wise contribution. Inspired by this, we refined the measures of capturing individual connectomes by increasing inter‐subject variability across FC values in a population to improve the power of predicting individuality with FC fingerprints (Cai et al., 2019). In addition, Amico et al. showed that the individual fingerprint of a human functional connectome could be improved from a reconstruction procedure based on a group‐wise decomposition (Amico & Goñi, 2018). Recently, the limitation of static connectivity has been widely realized, and the concept of dynamic connectivity has emerged to emphasize the time‐varying characteristics of the FC (Allen et al., 2014; Cai et al., 2017, 2018; Calhoun & Adali, 2016; Calhoun, Miller, Pearlson, & Adalı, 2014; Hutchison et al., 2013; Hutchison, Womelsdorf, Gati, Everling, & Menon, 2013; Zhang et al., 2019; Zhang et al., 2019; Zhang, Fang, Liang, Calhoun, & Wang, 2018). Incorporating the information from the time‐varying FC, Liu et al. studied whether and how the dynamic properties of the chronnectome acted as a fingerprint of the brain to identify individuals (Liu, Liao, Xia, & He, 2018). Mounting evidence indicates that during the dynamic FC analysis, a state (termed stable state) that resembles the static FC reoccurs more frequently than others (Allen et al., 2014; Hutchison, Womelsdorf, Allen, et al., 2013). These stable states are similar across subjects and share a basic configuration representing all dynamic patterns. Hence, we hypothesize that the common neural activities, without the individual‐specific characteristics in the connectomes, can be represented by the static FC. This may impede the revealing of individualized characterization which, to our knowledge, has not yet been considered when performing the individualized pattern analysis.

To this end, we use a dimensionality reduction technique to extract the basic configuration. More specifically, we project fMRI data onto its underlying subspace with a common structure. A simple and commonly used method is principal component analysis (PCA), which finds the direction of the greatest variance in the data set and represents each data point by its coordinates along each of these directions (Wold, Esbensen, & Geladi, 1987). However, PCA cannot extract nonlinear structures modeled by higher than second‐order statistics. Various methods have been proposed for nonlinear dimension analysis, such as the auto‐associative networks, generalized PCA, and kernel PCA (Karhunen & Joutsensalo, 1995; Kramer, 1991; Malthouse, 1998; Schölkopf, Smola, & Müller, 1998). The sophisticated autoencoders (AEs) (Hinton & Salakhutdinov, 2006) are typically a class of neural network‐based nonlinear dimensionality reduction methods. Through multilayer neural networks, the autoencoder and its extensions demonstrate powerful performance in learning key features from data (Lee, Ekanadham, & Ng, 2008; Rifai, Vincent, Muller, Glorot, & Bengio, 2011; Vincent, Larochelle, Bengio, & Manzagol, 2008), which found successful applications described previously (Hinton & Salakhutdinov, 2006; Le, 2013; Vincent, Larochelle, Lajoie, Bengio, & Manzagol, 2010). Thus, in this work, we use an autoencoder to estimate the common neural activity from rest/task fMRI data.

The remainder of this paper is organized as follows. In Section 2, we first describe the data set used in this work. Then, we introduce our proposed framework for estimating FCs step by step, followed by a series of experiments. In particular, we analyze whether refined connectomes extracted by our proposed method can better distinguish each individual from a pool of participants, and predict high‐level cognitive behaviors. The corresponding results are illustrated in Section 3. Some discussions and concluding remarks are given in Sections 4 and 5, respectively.

MATERIALS AND METHODS

Data acquisition

We used the publicly available S1200 Data Release of the human connectome project (HCP) (Van Essen et al., 2013). The S1200 release contains behavioral and 3 T MR imaging data from 1,206 healthy young adult participants collected from August 2012 to October 2015. A total of 889 subjects have complete data for all four 3 T MRI modalities in the HCP protocol: structural images (T1w and T2w), resting‐state fMRI (rsfMRI), task fMRI (tfMRI), and high angular resolution diffusion imaging (dMRI). A written informed consent was obtained for each subject. All HCP subjects were scanned on a customized Siemens 3 T “Connectome Skyra” housed at Washington University in St. Louis, using a standard 32‐channel Siemens receiver head coil and a “body” transmission coil designed by Siemens specifically for the smaller space available using the special gradients of the WU‐Minn and MGH‐UCLA Connectome scanners. To address head motion, dynamic head position information was acquired using an optical motion tracking camera system (Moire Phase Tracker Kineticor). fMRI was acquired using a whole‐brain multiband gradient‐echo (GE) echoplanar (EPI) sequence with the following parameters: TR/TE = 720/33.1 ms, flip angle = 90°, FOV = 208 × 180 mm, matrix = 104 × 90 (RO × PE), multiband factor = 8, echo spacing = 0.58 ms, slice thickness = 2 mm. The resulting normal voxel size was 2.0 × 2.0 × 2.0 mm.

The resting‐state runs (R1 and R2) were acquired in separate sessions on two different days. Task runs included the following: working memory (Wm), motor (Mt), language (Lg), and emotion (Em). The working memory task and motor task were acquired on the first day, while the language and emotion tasks were acquired on the second day. Note that not all participants have these six modalities of fMRI data. Thus, we filtered out subjects lacking one or two modalities of the fMRI scanning session. Following data selection, a cohort of 862 participants (aged 22–35 years, 409 male and 453 female) was included in our analyses. Within each session, oblique axial acquisitions alternated between phase encoding in a left‐to‐right (LR) direction in one run and phase encoding in a right‐to‐left (RL) direction in another run. Note that not all the participants have both left‐to‐right and right‐to‐left data. Here, we included only the left‐to‐right encoding runs to avoid potential effects of different phase encoding directions on our findings. More details about S1200 Data Release of the HCP can be found in the reference manual (WU‐Minn, 2017).

Data preprocessing

Our study used the fMRI data set from HCP with the minimal preprocessing pipeline, which included gradient distortion correction, head motion correction, image distortion correction, spatial normalization to standard Montreal Neurological Institute (MNI), and intensity normalization (Glasser et al., 2013). Further, we applied the standard preprocessing procedures to reduce biophysical and other noise sources in the minimally processed fMRI data. These procedures contained the removing linear trend and performing band‐pass filtering (0.01–0.25 Hz). Notably, Finn et al. has reported that the smoothing level had essentially no effect on identification accuracy (Finn et al., 2015). Thus, we investigated the analysis based on the data without applying spatial smoothing.

To facilitate the understanding of behaviors associated with different brain regions, we applied a 268‐node functional atlas provided by Finn et al. (2015)), which was defined using a group‐wise spectral clustering algorithm (Shen, Tokoglu, Papademetris, & Constable, 2013). More specifically, we extracted the time series of each node by averaging the time courses of all voxels that belonged to that node. Then, we assigned these nodes into eight functional networks, including medial frontal (Med F), frontoparietal (FP), default mode (DMN), subcortical‐cerebellum (Sub‐Cer), motor (Mt), visual I (Vis I), visual II (Vis II) and visual association (Vis Assn) regions. Axial, sagittal, and coronal views of these functional networks were displayed in Figure 1.

FIGURE 1

Axial, sagittal, and coronal views (from left to right) of 8 functional networks provided by Finn et al. 1. Med F: medial frontal; 2. FP: frontoparietal; 3. DMN: default mode; 4. Sur‐Cer: subcortical‐cerebellum; 5. Mt: motor; 6. Vis I: visual I; 7. Vis II: visual II; 8. Vis Assn: visual association. The figure was generated by the BioImage Suite Web

Autoencoder network construction

An autoencoder was used to extract common neural activity from the BOLD time series. For each participant of a modality, BOLD time courses with p ROIs and n time points (p, ) are available. These signals are marked as the original time series and set as inputs to the AE network. The AE network can be defined as: where the deterministic mapping f is called the encoder, which transforms an input vector into the hidden representation . Its parameter set is θ = {W, b}, where W is a d × d weight matrix with d  ≤ d and b is an offset vector of dimensionality d .

The resulting hidden representation y is then mapped back to a reconstructed matrix z with a dimension of d × p in the input space, that is, z = g (y). This mapping g is called the decoder. s 1 and s 2 are activation functions for encoding and decoding layers, respectively. Here, we used rectified linear units (ReLUs) in all encoder/decoder pairs, except for the first pair (Because the time series have both positive and negative values. Here, hyperbolic tangent was used) (Nair & Hinton, 2010). The training was performed by minimizing the least‐square error . After training of one layer, we used its output f (x) as the input to train the next layer. Note that the AE network is usually implemented to minimize the loss between the reconstructed output and input. In our case, we apply the residual to exclude the common patterns shared by subjects. Hence, we need to avoid over‐parameterized scenarios by selecting the parameters in the DNNs (deep neural networks). More details about the parameters selection can be found in (Van Der Maaten, 2009). More importantly, inspired by the work provided by van der Maaten (Van Der Maaten, 2009), we set layers' dimensions in AE as d‐500‐500‐2000‐10‐2000‐500‐500‐d for all fMRI modalities, where d is the number of time points (i.e., d = n ) varying with respect to different fMRI modalities. All layers are fully connected. Then, we calculated the residual time courses, which are the differences between the original time series and reconstructed ones generated by the AE network. Next, the residual time series were set as the inputs to the subsequent sparse dictionary learning (SDL) model. An illustration of the workflow is displayed in Figure 2a.

FIGURE 2

An illustration of the workflow to refine the brain connectivity. (a) Extraction and reduction of the effect of common neural activity using the autoencoder network. The dimensions of the AE network are set to be d‐500‐500‐2000‐10‐2000‐500‐500‐d for each participant. The difference between the original and reconstructed time series (residual time series) is set as the input to a sparse dictionary learning model (SDL). (b) Decomposition of the FC into both the group‐wise and subject‐wise patterns using the SDL model. Note the assumption that the subject‐specific FC may carry most of the identification information for a participant, which is tested and implemented here. The figure was generated by the NN‐SVG and Matlab

Increased individual identifiability using SDL model

In our previous study, we indicated that individual connectivity analysis benefits from group‐wise inferences and the refined connectomes are indeed useful for brain mapping (Cai et al., 2019). Thus, to further improve the inter‐subject variability across FCs, we implemented the same pipeline to reduce group‐wise contribution. Assume that we have subjects. For the residual time courses (p ROIs and n time points, ), we first calculate a correlation matrix, , for each subject. C (b 1, b 2) is the Pearson correlation between ROIs b 1 and b 2 for the i th subject across the entire residual time series. In consideration of the symmetry of the correlation matrix, we discard the upper triangular part of C . This leads to the edge weight vector for each subject. Next, we concatenate edge weight vectors from all subjects to form a matrix Y = [e 1, e 2, …, e ] with the size of m × n, where m = p(p − 1)/2. Identifying the sparse representation of the functional connectivity across subjects (Y) can be modeled as an SDL problem. By solving the following formulation, we can approximate the given data Y:where L is a nonnegative model parameter to control the sparse level of representations. denotes the dictionaries, and K is the size of dictionaries. is the representation matrix and ∥·∥0, ∥·∥ denote the l 0 and Frobenius norms, respectively. Here, the ksvdbox13 was implemented. More details about the SDL model can be found in our previous work (Cai et al., 2019).

Since we want to improve the inter‐subject variability, group‐wise contributions can be excluded from each correlation matrix C to obtain a new refined functional connectome . The refined functional pattern is defined as follows:where is the correlation matrix reconstructed from the lower triangular information Dx . The framework for the SDL model is illustrated in Figure 2b. Note that during analysis, the SDL model is performed on each subject of ROI network from different fMRI modalities individually.

Individual identifiability analysis

To explore the use of functional connectomes as fingerprints using our pipeline, we investigated individual identification ability proposed by Finn et al. (2015). Identification is performed across pairs of scans consisting of one target and one session from the HCP database, with the requirement that the target and database sessions are from different days to avoid interference as much as possible. For the target session with a given subject (e.g., resting‐state1, R1), we would like to identify that the connectivity pattern from the session in the database (e.g., language, Lg) belongs to the same subject. More specifically, for each participant, we first compare the correlation matrix of this subject from session 1 to each of the matrices of all the participants from session 2 (). For each comparison, the similarity scores between the connectivity patterns from sessions 1 and 2 are simply estimated using the Pearson correlation coefficient. Then, we assign this participant the same label with the subject in session 2 who has the maximal similarity score with this participant. If the FCs with the same label are indeed from the same participant, the identification accuracy is considered to be 100%. Otherwise, it is designated as 0%. By calculating the proportion of subjects with the correct identification, we determine the identification accuracy of all the participants. Finally, the sessions 1 and 2 are reversed, and the procedures discussed above are repeated (). Because we have three fMRI modalities (R1, Wm, Mt) for 1 ay and three modalities of fMRI (R2, Lg, Em) for another day, this results in nine possible combinations for (likewise, nine possible combinations for ).

After obtaining the identification accuracy for all the participants, we performed 10,000 nonparametric permutation tests (two‐sided) to assess whether the observed accuracies were significantly above chance. For each permutation, we randomize the identities of the subjects in the data set session, perform the identification procedures, and record the accuracies. A significant level of p‐value = .05 is used as the threshold for the 10, 000 permutation tests.

We then investigated the identification accuracy on the basis of each specific functional network to figure out which brain network contributes more to individual discriminability. These functional networks are defined in the section of data preprocessing. During this process, a single network or a combination of networks are used to estimate the individual identification. Note that, if we denote the set of nodes belonging to network j as V = v , k = 1, 2, …, K , where K is the total number of nodes in network j, only connections within the selected network are included.

Common neural activities contributing to individual identification

To check the hypothesis that common neural activities may weaken the individual variability, we compare the performance of identifiability with and without the AE network processing. For this purpose, we exclude the SDL model in this experiment to avoid its influence. As a first pass, we calculate correlations between connectivity matrices of all participants across nine possible combinations () under these two scenarios (with and without the AE network). For each correlation matrix, the row and column are symmetric. Thus, diagonal elements are similarity scores from the matched subjects, while off‐diagonal elements are the ones from the unmatched participants. By observing the difference between the mean values of diagonal and off‐diagonal factors, the individual identifiability can be evaluated. The larger the difference, the stronger the discriminative power.

We then estimate the identification accuracy for these two scenarios, respectively. The procedures have already been provided in the section of individual identifiability analysis. If the identification rates generated by the scheme with the AE network are much higher than those without the AE network, we assume that using AE can therefore increase the subject‐specific identifiability by reducing the effect of the common neural activities. Finally, we reverse sessions 1 and 2 and repeat the above proceedings ().

To further explore the effects of applying the AE network, we analyze the reconstructed FCs from the AE's output. First, we investigate the similarity between original‐residual FCs and original‐reconstructed patterns. Here, we use the Frobenius norm to calculate the distance between two matrices. We then calculate the identification accuracy by using the output of the AE networks. If the common neural activities do not carry the subject‐specific information, the identification accuracies should be much lower than that based on applying the original static FCs. In this way, we can figure out whether the residual time series possess meaningful individual information or are pure noise.

Afterwards, we investigate the regions in the FC, to which the signals removed by the AE network belong. To filter out the influence induced by activities in the task runs, we restrict the analysis to resting‐state fMRI (R1 and R2). Group differences between the analyses with and without AE are considered here. First, we calculate the correlation matrix for each subject. Next, we transform correlation matrices into the edge‐weight vectors (e ), and concatenate them into the data Y as mentioned in the previous section. Finally, a two‐sample t test is applied to the data Y with a significant level of q = 0.01 to examine the group differences.

Examination of AE network's effectiveness in capturing common neural activities

To assess the AE network's effectiveness in capturing the common neural activities, we implement the group‐average time series to make a comparison. More specifically, we first calculate the group‐average time series for each fMRI modality and hypothesize that these signals also represent the shared FC patterns. We remove these group‐average signals from the original time series of each subject. Next, we use the remaining time courses to repeat the identification experiments across all nine pairs of database and target sessions. We expect to see the result that the identification rates provided by the remaining time courses are lower than those obtained using the AE networks. Otherwise, the AE network in our pipeline may be unreasonable, since the common neural activities can be estimated using the group‐average signals. Meanwhile, the effectiveness of the AE network in extracting the common neural activities can be validated.

Edgewise contribution to identification

To investigate which connections of the FC contribute more to subject identification, we estimated the modified differential power (DP) provided by Liu et al.(Liu et al., 2018). In this part, we also restrict the study to R1 and R2 sessions. The modified differential power is defined as follows:where P (i, j) is an empirical probability to quantify the differential power of an edge for the purpose of subject identification; l and k (l ≠ k) represent the labels of two different participants; i and j (i ≠ j) denote two different nodes within the functional connectivity; n is the total number of participants in the analysis (n = 862). |ϕ (i, j) > ϕ (i, j)| indicates the probability that |ϕ | between two different subjects is higher than |ϕ | of the same participant. Given two sets of connectivity matrices obtained from the R1 and R2 sessions after z‐score normalization, the corresponding edge‐wise product vector ϕ (i, j) can be calculated as follows:and |ϕ | can be obtained in the same way. DP reflects each edge's ability to distinguish an individual subject. For a given functional connectivity, a higher DP value means a greater contribution to individual identification. Furthermore, to investigate the network‐dependent contribution to subject‐specific identification, we also count the number of the highest DP values (top 1%) within or between functional networks. In this manner, we examine whether specific brain networks play a significant role in discriminating individuals.

Prediction analysis for individual cognitive behavior

To determine whether our refined FC applying the AE network could benefit individual cognitive prediction, we apply the regression analysis on the cognitive behaviors. Here, we select items of high‐level cognition from the HCP protocol, including fluid intelligence (Penn progressive matrices, HCP: PMAT24_A_CR, mean ± SD: 17.04 ± 4.71, range: 4–24), cognitive flexibility/executive function (Dimensional Change Card Sort, HCP: CardSort_AgeAdj, mean ± SD: 102.54 ± 9.89, range: 57.79–122.65), inhibition/executive function (Flanker task, HCP: Flanker_AgeAdj, mean ± SD: 102.05 ± 9.94, range: 72.81–123.56) and language/vocabulary comprehension (Picture Vocabulary Test, HCP: PicVocab_AgeAdj, mean ± SD: 109.44 ± 15.07, range: 68.68–153.09). More details can be found at the HCP website (https://db.humanconnectome.org/). Then, the refined FCs from the R1 session are applied as the features to define these high‐cognition scores.

To determine whether the refined FC profiles can better predict the individual cognitive behavior relative to those without using the AE network, we use the leave‐one‐subject‐out cross‐validation (LOOCV) strategy to estimate the prediction accuracy (Cai et al., 2019). For instance, to assess the ability of refined FCs to describe fluid intelligence, in each LOOCV fold, one participant is used as the test sample, and the remaining n − 1 subjects are considered as the training samples. First, we concatenate all the connections within the FC profiles (i.e., 35,778 connections) to generate a feature vector for each subject. Second, we investigate a feature selection step, which calculates the correlation between each connection of FCs (Pearson correlation between two ROIs) and fluid intelligence scores on the training set. If the correlation is significant (p‐value <.001), the corresponding feature is retained. Third, a predictive model is built using a simple regression model to fit the selected features to the fluid intelligence coefficients in the training set. Finally, we adopt the model on the unseen test data to generate the behavioral score. During this procedure, each participant is used as the test sample once. After all the LOOCV folds are completed, we assess the predictive power through the correlation values between the predictive and observed fluid intelligence scores.

At the end, we perform the permutation test (10,000 times) to evaluate the statistical significance of the observed behavioral scores. For each permutation, the observed behavioral scores of the subjects are randomly shuffled before the regression analysis. In this way, we can examine whether the prediction performance is obtained by chance.

RESULTS

Refined FC based individual identification

As a first pass, we evaluated the identification accuracy applying the whole‐brain connectivity matrices (268 nodes, without prior network definitions) to validate that the refined FCs can highlight the subject‐specific variability. Identification rates are described in Figure 3a. Even with a large number of participants (862 subjects), refined FCs worked well for individual identification. The success rates were 99.3 and 99.5% based on a database‐target rest1–rest2 and the reverse rest2–rest1, respectively. Meanwhile, the identification rates ranged from 93.1 to 96.4% for rest–task pairs and 93.8 to 96.3% for task–task combinations. Relative to raw FCs (900 subjects, 89% accuracy between the two rest sessions, 65 ± 14% for other pairs in Finn et al.'s work [Finn et al., 2017]), refined FC profiles largely improved the performance of individual discrimination. To make the comparison directly, we repeated the identification experiment by applying the two approaches. More details can be found in Part 5 of the supplementary materials. We observed that our proposed pipeline achieved the highest identification accuracy for all the sessions' combinations among these three different strategies. Given that identification trials were not independent of each other, we performed 10,000 nonparametric permutation tests (two‐sided) to assess the significance level of these results. Across 10,000 iterations, the p‐value for each pair of sessions is below .0001. It indicates that the success rates of identification are significantly above chance. The overfitting is often a critical issue in practice, but we demonstrated that there is no such concern with our framework according to a cross‐validation experiment. More details and results are shown in Figure S3. By checking the identification rates in Figure S3, we found that the accuracy for all the scenarios had been significantly improved, and the performance was almost the same as that described above. Note that inappropriate setting of the parameters involved in the SDL, that is, dictionary size and sparse, will induce overfitting and lower identification rates for the proposed framework (Cai et al., 2019). Thus, we believe that our results are not primarily influenced by data leakage.

FIGURE 3

Identification accuracy across session pairs and networks. (a) Identification rates from the whole brain are highlighted in color‐coded matrices to compare the accuracies across rest–rest, rest–task, and task–task sessions, respectively. R1, Rest1; Wm, working memory task; Mt, motor task; R2, Rest2; Lg, language task; Em, emotion task. (b) Identification results based on all nine sessions in the database and target combinations. Each row shares the same database session and each column shares the same target session. The color of the bar (gray or white) indicates which fMRI modality was used as the database, and the other one was served as the target. Graphs display the identification rate based on each network as well as the combination network 1 and 2 and the whole brain (all). 1. Medial Frontal; 2. Frontoparietal; 3. Default Mode; 4. Subcortical‐cerebellum; 5. Motor; 6. Visual I; 7. Visual II; 8. Visual Association. The figure was generated by the Matlab

Next, we examined the identification accuracy based on each functional network to explore which brain network contributes more to the individual variability. These networks are defined by Finn et al. (2015) and depicted in Figure 1. The medial frontal network (network 1) and the frontoparietal network (network 2) achieve the highest success rates in individual discrimination, which comprises the higher‐order association cortices in the frontal, parietal and temporal lobes. In comparison with the medial frontal network, the frontoparietal network performs much better especially for the rest–task and task–task pairs. Furthermore, we also checked whether the combination of networks 1 and 2 can provide better performance than each individually. As shown in Figure 3b, in each scenario, the identification accuracies using the combination of networks 1 and 2 are higher than implementing network 1 or 2 independently, and are pretty close to those generated through applying the whole‐brain nodes. For other networks, subcortical‐cerebellum network (network 4) and motor network also contribute to subject‐specific variability.

Contribution of common neural activities to the individual identifiability

To investigate the contribution of common neural activities to individual identifiability, we estimated the correlations between connectivity matrices of all subjects across nine possible combinations for both time courses with and without the AE network processing. For each correlation matrix, the row and column are symmetric by subject. Thus, diagonal elements are correlation coefficients from the matched subjects, while off‐diagonal elements are those from the unmatched participants. The results for all nine possible pairs are displayed in Figure 4a. In comparison with raw cross‐subject correlation coefficients, scores generated after the AE network become significantly weaker for both diagonal and off‐diagonal factors. However, applying the AE network improves the difference between diagonal and off‐diagonal elements in the correlation matrix. It indicates that reducing the contribution from common neural activities mentioned above indeed helps individual discrimination.

FIGURE 4

Evaluation of the influence of common neural activities on individual identification. (a) Analysis of identifiability matrices based on all nine pairs of the combination of a session in the database and the target (). For the top line of each sub‐figure, from left–right: identifiability matrix (i.e., Pearson correlation coefficient between functional connectivity across subjects and modalities) of the original data; identifiability matrix of the AE residual data. The row and column subject order of identifiability is symmetric. Hence, diagonal elements are correlation scores from the matched subjects, while off‐diagonal coefficients are from the unmatched participants. Mean correlation coefficients for both diagonal (match) and off‐diagonal (unmatched) elements are also calculated (bottom, err bars indicates ± SD). ** means p‐value<10−5 for two‐tailed t test. Besides, the difference (mean value) between the diagonal and off‐diagonal elements for both the original and AE residual data are displayed. (b) Comparison of identification accuracies across all nine pairs of database session and target session between using original time series (top) and AE residual data (bottom). Note that only the situation with whole‐brain nodes (264 nodes) is considered here. (c) Static functional network connectivity for the original data (left) and AE residual signals (middle) estimated by the Pearson correlation. Meanwhile, group differences (right) between them is obtained by applying a two‐sample t test. These variations are visualized by plotting the log of p‐value with the sign of t statistics, − sign(t)log 10(p). Note that in this part, the results from the rest1 and rest2 are displayed individually. The figure was generated by the Matlab

Next, we repeated the identification experiments applying connectivity matrices from the two scenarios discussed above to further validate our hypothesis. By checking the identification results in Figure 4b, we observe that identification accuracies for all nine pairs have increased through reducing the common neural contribution. More specifically, for the pairs of rest–rest, the identification rates were improved around 4%. As to rest–task combinations, when using connectivity matrices of resting‐state fMRI in the database, the ability of discrimination was significantly enhanced (R1‐Lg: 36.2%, R1‐Em: 27.3%, R2‐Wm: 31.4%, R2‐Mt: 23.6%). By contrast, the rates applying matrices achieved from task‐based fMRI in the database gained around 4%. Meanwhile, with the AE network preprocessing, the identification accuracies increased for every condition of task–task combination, ranging from 11.2% (Em‐Mt) to 25.7% (Lg‐Wm). Thus, weakening the signals from common neural activities can significantly enhance inter‐subject differences, and the pairs of rest–rest possess the strongest identification power.

We then further analyzed the output of the AE network. By observing the results in Figure S1a,b, we find that the FCs generated by the residual time series are much closer to the original ones relative to those extracted from the AE's output. This means that the residual time series indeed carry subject‐specific information instead of noise. Thus, the feasibility of our pipeline is verified. We also checked the identification ability of those FCs constructed by the AE's output. As shown in Figure S1c, the highest accuracy is 20.8% (R2–R1). Most of the rates under other combinations cannot achieve 10%, which further proved our conjecture that the common neural activities may weaken the individual discrimination.

To explore which connections were removed by the AE networks, we tested the group average functional connectomes before and after applying the AE networks and examined the difference between them using a two‐sample t test. From Figure 4c, we obtain that the strength of links in the connectomes reduces overall. However, the significant difference is largely related to the frontoparietal and subcortical‐cerebellum networks.

The effectiveness of AE network in capturing common neural activities

To validate the AE network's effectiveness in extracting common neural activities, we first obtained the shared patterns across subjects applying the group‐average signals for each fMRI modality. Next, we estimated the remaining time series by subtracting the group‐average signals and repeating the identification experiments discussed above. Finally, we compared the results with those generated by the AE networks. The results are shown in Figure S2. We observe that the FC patterns constructed by the group‐average signals have more details relative to those from the AE network. In our opinion, some of the subject‐specific features may be weakened by implementing this simple approach. As to the identification experiment results, the rates are close to the accuracy using the original time series. It shows that the group‐average method cannot extract common neural activities appropriately. For specific pairs (e.g., Wm‐Lg, Mt‐Em), the accuracy is lower than the one from the original FC patterns. However, for the same fMRI combinations, removing the common neural activities estimated by the AE network has significant improvement in the accuracies (for instance, from 60.7 to 72.0% for Wm‐Lg, from 40.1 to 63.2% for Mt‐Em). It indicates that the AE network outperforms group‐average approach for estimating the common neural activities.

Evaluation of edgewise contributions to the individual identification

To determine which connections contribute more to subject‐specific identification, we calculated the modified differential power (DP). The modified DP reflects each edge's ability to distinguish an individual from a pool of participants. A connection with high DP tends to have a similar value within an individual across modalities, but possesses a different degree across individuals regardless of modalities.

By restricting the analysis to resting‐state fMRI (rest1 and rest2), we estimated the modified DP for all edges in the brain. We determined which connections were in the 99.9 percentile across all of the links (Figure 5). We observe that the majority of edges in the 99.9 percentile of the edges are in the frontal, parietal, and temporal lobes. Meanwhile, most of the nodes with high DP values were found in the frontoparietal and medial frontal networks. Some of them belong to the default mode network. By checking the results in Figure 5b, we find that for the connections possessing high DP values in connectivity, 27.3% of them link the medial frontal and frontoparietal networks. In addition, 59% are connections linking these networks to others (38% of them link the frontoparietal network to other networks). It indicates that connections related to high‐order association cortices are the most discriminative of individuals. On the other hand, medial frontal and frontoparietal networks play a significant role in individual identification.

FIGURE 5

Edgewise contributions to individual identification. (a) Connections that possess the highest DP scores in individual connectivity profiles (top, circle plot). Axial, sagittal, and coronal views of these links are also provided (bottom, from left to right). Note that connections with the highest 1% DP values are shown here. In the circle plots (top), the 268 nodes (the inner circle) are organized into a lobe scheme (the outer circle) roughly reflecting brain anatomy from anterior (top of the circle) to posterior (bottom of the circle) and split into left and right hemisphere. Lines indicate edges or connections. (b) The percentage of connections within and between each pair of networks (eight functional networks defined in Figure 1) using the same data as (a). The color depth of the grid in the matrix indicates the fraction of DP edges for each pair of networks. The figure was generated by the BioImage Suite Web and Matlab

Connectivity profiles predict high‐level cognitive behaviors

To test whether refined FC profiles improve the behavior prediction, we explored the prediction abilities across different high‐level cognitive scores of connectivity profiles under the scenarios with and without the AE network. Note that for these two conditions, the SDL model was included during the analysis. As demonstrated in the scatter plots in Figure 6, the predicted scores by the connectivity profiles with the AE network have higher correlation coefficients with the observed scores relative to those without the AE network. Besides, the range of predicted scores by the refined FCs is much narrower across all cognitive scores (especially for language/vocabulary comprehension). To validate the prediction power of our proposed framework applying the AE network, we performed 100 nonparametric permutations for each score. In addition, we indicated that the overfitting did not significantly affect the linear regression as depicted in Figure S4. The results illustrate that the prediction of each high‐level cognitive behavior (the correlation between observed and predicted scores) is above chance (fluid intelligence: p‐value<.01; cognitive/executive flexibility: p‐value<.08; inhibition/executive function: p‐value<.06; language comprehension: p‐value<.01). It means that reducing the common neural activities also helps predict cognitive behavior.

FIGURE 6

Connectivity profiles predict cognitive behavior. Scatter plots display prediction results from a leave‐one‐subject‐out cross‐validation (LOOCV) analysis comparing the predicted and the observed high‐level cognitive scores. Both the connectomes with and without applying the AE network processing are considered. Note that under these two scenarios, the sparse dictionary learning (SDL) model is included. In the scatter plot, each dot represents one subject, and the area between dashed lines reflects 95% confidence interval for the best‐fit line, which is used to assess the predictive power of the model. R‐values are the correlation coefficients between the predicted and observed high‐level cognitive scores. Furthermore, for each cognitive scores, edges retained from the feature selection step (p‐value<.001) are also depicted (with the AE network process). Axial, sagittal, and coronal views of these connections in the brain are provided. In the circle plots, the 268 nodes (the inner circle) are organized into a lobe scheme (the outer circle), roughly reflecting brain anatomy from anterior (top of the circle) to posterior (bottom of the circle), and divided into left and right hemisphere. Red and blue lines mean positive and negative connections, respectively. (a) fluid intelligence; (b) cognitive flexibility/executive function; (c) inhibition/executive function; (d) language/vocabulary comprehension. The figure was generated by the BioImage Suite Web and Matlab

Furthermore, by observing the selected features from refined FC profiles, we find that different brain regions exhibit distinct contributions to various high‐level cognitive parameters. More specifically, frontal and parietal lobes contribute most to fluid intelligence prediction according to this study. In detail, nodes located in the frontal and parietal regions provide positive connections highly related to fluid intelligence. Moreover, some of the positive links are relevant to the right cerebellum, while a large portion of negative links connects with the insula region. When using the refined FC profiles to predict the performance of cognitive flexibility, in comparison with positive edges, negative connections play a more critical role. Most of these negative edges are linked with frontal and parietal lobes. Notably, the region of the motor strip is involved in both the positive and negative connections regarding cognitive flexibility. Similarly, negative connections predominantly contribute to inhibitory function as well. Among these negative links, the majority of them interact with the parietal lobe. In addition, the negative edges from the motor strip and subcortical regions have equivalent influence on the prediction of inhibitory function. In the prediction of language comprehension, through both the positive and negative edges in the connectivity, the frontal and temporal lobes closely connect with language comprehension. As a consequence, the data‐driven framework applying the refined FC profiles can adequately enable us to find out which brain regions closely interact with high‐level cognitive behavior.

DISCUSSION

Recently, the study and use of dynamic functional network analysis has drawn more attention. Some basic configurations in brain connectivity appear across all time‐varying states, and may not contribute to subject‐specific discrimination. In our previous work, we pointed out that refining the measures of individual connectomes can help improve the fingerprinting power of individuality (Cai et al., 2019). In this study, we applied the autoencoder network to extract the signals of common neural activities, which may hinder the identification of individual differences. Therefore, we removed them from the raw time courses when constructing connectivity network between ROIs. We then applied dictionary learning on all the data. Based on the pipeline, we generated the refined FC connectomes (Cai et al., 2019). We showed that the refined FC profiles can successfully distinguish one individual from a pool of the population with a high identification accuracy. Moreover, the refined connectomes can significantly predict high‐level cognitive behavior, including fluid intelligence, cognitive flexibility, inhibition function, and language comprehension. Notably, reducing the signals of common neural activities benefited individual identification and cognition prediction, where frontal, parietal, and temporal lobes contributed significantly. Collectively, our findings supported our assumption that common neural activities may impede the individual identifiability and its removal can help enhance the uniqueness of the connectome for each individual.

When testing the identifiability of using refined FCs, we found that regardless of database‐target combinations, the connectomes applying the proposed approach successfully distinguished an individual from all the participants. The accuracy was from 99.3 to 99.5% for the rest–rest pairs. As to the rest–task and task–task combinations, the success rates ranged from 93.1 to 96.4% and from 93.8 to 96.3%, respectively. The strong individual identification power of the refined FCs from our results suggest that connectomes vary across participants and are unique to each subject. This is in agreement with previous findings that the connectomes can be used as fingerprints to identify an individual (Cai et al., 2019; Finn et al., 2015). Compared to results obtained by raw FCs, reducing the effect of common neural activities and group factors can significantly improve the identification rates (Finn et al., 2017). For the pairs of rest‐rest fMRI, we gained 10% accuracy improvement. Also, the success rates increased by around 30% for other combinations. The findings validate our assumption that some patterns caused by common neural activities may impede the individual variability, and thus reducing their impact helped capture unique characteristics of the brain.

The contributions of functional networks to individual identification were also examined. Although discrimination based on the whole connectivity matrix performed best, the medial frontal and frontoparietal networks achieved sufficiently high accuracy. The combination of these two networks provided better performance than only one of them alone. These networks are composed of higher‐order association cortices (frontal, parietal, and temporal lobes), which have been proven to show the highest inter‐subject variance (Cole et al., 2013; Finn et al., 2015). Relative to the medial frontal network, frontoparietal contributed more to the identification, especially for rest–task and task–task pairs. It is consistent with the function of the frontoparietal system, which is particularly active in tasks requiring a high degree of cognitive control. Even for the rest–rest combination, the frontoparietal network worked very well for identification. Thus, we believe that the frontoparietal system plays a significant role in determining each brain's uniqueness regardless of whether the mind is at rest or not. In addition, we detected that the subcortical‐cerebellum and motor network correlated positively with individual differentiation. It matches with the conclusion that there was a gradual increase of variability in primary regions of the visual and sensorimotor systems specific to subcortical and cerebellum structures as the brain developed (Li et al., 2017).

To further estimate the contribution of common neural activities to individual variability, we compared the correlation matrices with and without AE network, and repeated the identification experiments for these two cases. While applying the AE network weakened both diagonal and off‐diagonal factors of correlation matrices, reducing the signals from common connectivity patterns increased the difference between diagonal and off‐diagonal elements. Besides, without using the SDL model, reducing the common neural contribution increased the identification rates across all nine combinations. All these findings indicate that weakening the contribution of common neural activities helps strengthen the subject‐specific variability and enhances the uniqueness of each brain. Furthermore, the identification rates generated by the AE's output are at most 20.8%. Compared to the FCs constructed by the AE's output, those extracted by the residual signals are more similar to the original patterns. It demonstrates that the residual time series possess the subject‐specific information and validates the feasibility of our approach. Intriguingly, we observed that the time courses identified by the AE networks are mainly related to the frontoparietal and subcortical‐cerebellar systems, which have a significant influence on the identification of individuals. In light of this, the interaction between these networks with individual predictions will be studied in future work.

While applying the simple group‐average method to capture the common neural activities, we discovered that the individual identification ability was not significantly improved (for some fMRI combinations, the rates were even reduced). We believe that the common neural activities in time series cannot be estimated easily by the group‐average method. In contrast, the AE network performs much better in detecting common neural activities, facilitating the extraction of subject‐specific features.

By examining the modified differential power of each edge in the connectivity map, we determined that most of the connections with high DP values were related to the higher‐order association (frontal, parietal, and temporal) lobes. In addition, 27.3% of high DP edges were connected with the medial frontal and frontoparietal systems, and 57.3% of them linked these networks to others. These findings further validate the function of the medial and frontoparietal networks in individual identification.

When exploring the prediction abilities across different cognitive parameters with and without AE processing, we observed that the predicted scores by applying AE processing had higher correlation coefficients with the observed scores than without AE processing. This suggests that the reduction of common functional patterns can improve the power of cognitive behavior prediction. In addition, we analyzed the selected features that contributed more to the cognition prediction and found that different brain regions had various effects on a cognitive measurement. We demonstrated large contributions of the frontal and parietal lobes to individual fluid intelligence and execution function (e.g., flexibility and inhibition function). These findings are consistent with previous studies (Cole et al., 2013; Finn et al., 2015; Liu et al., 2018). Specifically, we detected that insula was closely associated with fluid intelligence. This agrees with the statement that fluid intelligence is correlated with a distributed network comprising regions of the frontal, insula, and parietal cortex (Tschentscher, Mitchell, & Duncan, 2017). Furthermore, we found that the motor strip region was highly related to executive function. This point is supported by the argument that motor and cognitive processes are functionally related and most likely share a similar evolutionary history. It is well established that multiple brain regions integrate both motor and cognitive functions (Leisman, Moustafa, & Shafir, 2016). For language comprehension, we demonstrated that the frontal and temporal lobes closely interacted with it. Broca's area (in the frontal lobe) and Wernicke's area (in the temporal lobe) are cortical areas that respond to human language. In sum, connectomes refined by our approach can improve our ability to characterize the relationship between brain regions and cognitive behaviors, which can in turn enhance our understanding of the human brain.

Several issues need further consideration. First, in this work, we still used static functional network connectivity. However, recent research has noted that dynamic functional connectivity could provide complementary individual information (Liu et al., 2018). A combination of both dynamic and static connectivity is a promising direction for analyzing individual variability. Second, several confounding factors such as parcellation schemes may affect the performance of our proposed approach. Follow‐up studies are needed to further study the effects of these factors. Third, we focus on group‐wise common and individualized aspects of the connectome. More work is needed to fully study their associations with brain functions.

CONCLUSION

In this work, we assumed that the removal of common neural activities can enhance the difference in brain connectivity across participants. Using autoencoder, we reduced these common patterns of connectivity to increase the uniqueness of connectome for each individual. We observed that refined FC profiles by our proposed pipeline can identify each individual with high accuracy (up to 99.5% for the rest–rest pair). Moreover, they can also be used to predict high‐level cognitive behaviors (e.g., fluid intelligence). We also showed that frontal, parietal, and temporal lobes contributed more significantly to brain functions. In summary, the findings in this study validate our hypothesis, and our proposed approach provides a promising way to study individualized brain networks.

Appendix S1 Supporting Information.

Click here for additional data file.