Hill Climbing Optimized Twin Classification Using Resting-State Functional MRI
Twin imaging studies are an important part of human brain research that can reveal the importance of genetic influences on different aspects of brain behavior and disorders. Accurate characterization of identical and fraternal twins allows inference to be performed on the genetic influence in a population. In this paper, we propose a novel pairwise feature representation to classify the zygosity of twin pairs using resting state functional magnetic resonance images (rs-fMRI). Specifically, we project an fMRI signal to a set of cosine series basis, and use the projection coefficients as the compact and discriminative feature representation of noisy fMRI data. The pairwise relation is encoded by a set of twin-wise correlations between the new feature representations across brain regions. We further employ Hill Climbing variable selection to identify the brain regions that are most genetically affected. The proposed framework has been applied to 208 twin pairs in the Human Connectome Project (HCP) and we achieved classification accuracy in determining the zygosity of paired images.
keywords:Twins, heritability, resting-state fMRI, sparse model, Hill Climbing, Human Connectome Project
The development and functioning of human brain are influenced by both genetic and environmental factors. However, the extent by which these factors shape brain structure and function is still unknown even though the topic has been of great research interest for decades. A substantial new advance in this area has been recently made and the Recently, a dataset containing brain images of twins has become available through the Human Connectome Project (https://www.humanconnectome.org, HCP). Twins provide a valuable source of information, as their unique relationship allows researchers to pull apart and examine genetic and environmental influences in-vivo. The power of twin studies arise from the fact that there are only two types of twins, identical (or monozygotic, MZ) and fraternal (or dizygotic, DZ), that share different amount of genetic information: MZ twins are expected to share 100% of genes (up to the amount of possible de novo somatic mutations (McConnelleaal1641)), and DZ twins on average are expected to share only 50% of genes (Twins). Comparing the similarity of MZ and DZ twins, we can draw inference on the genetic influence in a population.
The differences in genetic makeup between MZ and DZ twins provide the fundamental basis for studying phenotypical brain variations. The significance of genetic contribution in twin studies has been reliably shown for a wide range of structural brain measures, including brain size (Pennington2000), grey matter volume (Thompson2001) and neuroanatomical shape (Ge2016), cortical surface area (YOON2012169; Chen1634) and thickness (Schmitt6774; Strike2018), white matter microstructure (BROUWER20101085; JAHANSHAD2013455) and whole-brain white matter architecture (Yeh2016; KUMAR2017242) and development (CHIANG20112308), fractional anisotropy (blokland2012; KOCHUNOV2015300), as well as structural brain network connectivity (SHEN2014628; Bohlken2014).
An extensive review of the literature has shown that the vast majority of twin-related research in brain imaging is still primarily devoted to the study of heritability (measure of genetic influence) of brain microstructure and morphometry (TwinStudy2007; TwinStudy2008). However, in recent years an increasing number of studies have focused on measuring genetic influence on brain function (TwinStudy2015; TwinStudy2016), e.g., detecting heritability of neural activation during simple visuomotor task (PARK20121132), calculation (PINEL2013306), oral reading recognition and picture vocabulary comprehension (Feremi2017), estimating genetic contribution to brain activation in neural networks supporting working memory tasks (KARLSGODT2007191; BLOKLAND200870; Blokland10882; Koten1737). Though all of these studies are mainly focused on investigating the genetic influence of task-based brain activity, little research has been performed using resting-state brain activity. To the best of our knowledge, there is only a handful of twin studies using resting-state functional MRI that show significant genetic contributions to functional network connectivity architecture of the human brain (TwinStudy2010; Fornito3261; VANDENHEUVEL201319; Gao11288; Yang2016). However, these studies did not investigate the utility of the possible genetic effects at the voxel level.
The interest in using resting-state fMRI in contrast to task-related fMRI is explained by the fact, that it can be used to investigate the baseline functional connectivity of the brain (rsfMRI1; rsfMRI2). Therefore, it has a broader use for clinical applications, including the study of the human brain development and function (Cote2007; Vergun2013), and early detection and treatment of different neurological and mental pathologies (greicius2008resting; HUANG2010935; Jie2014; DAMARAJU2014298; rsfMRI3; RestingState; Wee2016). The identification of highly heritable brain regions can serve as the first step for uncovering genetic variants that contribute to neurological disorder or psychopathology. Findings revealed in this paper can also help to provide a better understanding of the recently identified physiological mechanism responsible for resting fMRI functional connectivity (MATEO2017936).
In the past, a wide range of machine learning methods, including Artificial Neural Networks (ANN), Support Vector Machines, -Nearest Neighbor, Gaussian Naïve Bayes and Fuzzy Classifiers (wang2004training; singh2007detection; PEREIRA2009S199; Peltier5332805; Vergun2013; Wang2014; honorio2015classification; wang2017depression) have been used to analyze fMRI data. However, none of them have been used in classifying paired twin images. It is unclear how to apply existing classifiiers to twin rs-fMRI. In this paper, we propose a unified classification framework that automatically determines the zygosity of a pair of twins using resting-state fMRI and identifies brain regions with significant genetic influence. The framework combines the algebraic representation of fMRI time series at the voxel-level and the sparse version of supervised machine learning algorithm that utilizes a multi-layer Artificial Neural Network (Figure 1).
For a given non-trivial classification problem setting, where each classification identity is represented by a pair of twin images, a special type of neural networks is required, specifically designed to process paired data through a comparative analysis. This type of neural networks is called Siamese Neural Network (SNN) (Siamese1; zagoruyko2015learning) and has recently gained interest from the medical imaging society (WANG2017137; kouw2017mr; Ktena2017; WANG201812; KTENA2018431). Whilst one of the main features of siamese neural networks is the ability to learn the suitable embedding space of the original data, this embedding space is usually constructed using highly abstract and semantic information from the original data, which are difficult to interpret biologically in a meaningful fashion. In order to have biological interpretability, we propose to adopt the idea behind SNN but instead of learning paired representation through a blackbox, we introduce a reliable parametric algebraic framework that provides an easy-to-comprehend twin representation.
In this paper, we propose a dimension reduction approach to construct a new compact feature representation for the fMRI signal at each voxel using the cosine series representation (CSR). The correlation between CSR is feed into in a two-layer ANN classification model. Using Hill Climbing, a variable selection procedure, along with sparse-ANN, we determine the significance of AAL regions and discover the most heritable regions.
The main contributions of the paper are: (1) we introduce a new feature representation to characterize twin-wise relationship in the whole brain that substantially improves the classification performance; (2) we attempt for the first time to identify zygosity using rs-fMRI data only, we use a sparse classification solution to classify MZ and DZ twins with accuracy; (3) we propose a principled way to determine the most genetically heritable regions.
The proposed framework has proved its robustness in one of the largest publicly available twin datasets – HCP that contains rs-fMRI of 131 MZ and 77 DZ twin pairs.
2 Material and Methods
We use resting-state fMRI scans from 208 healthy same-sex twin pairs. This data was collected as part of the Washington University-Minnesota Consortium Human Connectome Project (VANESSEN20121299; HCP; VANESSEN201362). All participants gave informed consent. The information regarding the zygosity of twin pairs is determined by genotyping. There are 131 MZ twin pairs (age , 56M/75F) and 77 DZ twin pairs (age , 30M/47F). Demographic details of the subjects are given in Table 1.
|Twin group||Sample size||Sex (M/F)||Age|
All subjects were scanned on a customized Siemens 3T Connectome Skyra scanner housed at Washington University in St. Louis, using a standard 32-channel Siemens receive head coil and a customized SC72 gradient insert and a customized body transmitter coil with 56 cm bore size. Resting-state functional MRI were collected in four runs of 14 minutes and 33 seconds using a gradient-echo-planar imaging (EPI) sequence with multiband factor 8, time repetition (TR) 720 ms, time echo (TE) 33.1 ms, flip angle , (ROPE) matrix size, 72 slices, 2 mm isotropic voxels, and 1200 time points. Runs were combined in two sessions, and within each session, oblique axial acquisitions alternated between phase encoding in a right-to-left (RL) direction in one run and phase encoding in a left-to-right (LR) direction in the other run. During each scanning participants were at rest with eyes open with relaxed fixation on a projected bright cross-hair on a dark background and presented in a darkened room (HCP1200). In this study, we use fMRI of different phase-encoding directions (LR and RL) from the first scanning session.
To benefit from the high-quality data provided by the Human Connectome Project, we used fMRI scans that undergone spatial (GLASSER2013105) and temporal (SMITH2013144) preprocessing, including: removing spatial distortions, realigning volumes to compensate for subject motion, registering the fMRI data to the structural MNI template, minimal highpass temporal frequency filtering, independent component analysis (ICA)-based artefact (non-neural spatiotemporal components, e.g. motion parameters) removal. The resulting volumetric data contains resting-state functional time series with 2-mm isotropic voxels at 1200 time points.
2.2 Cosine Series Representation
We start with normalizing fMRI signal at voxel and time by subtracting by its mean over time so that . Then normalized fMRI time series is modeled as
a function of the scaled time variable , where is a zero-mean noise and is unknown signal. We scaled the 1200 time points of fMRI into interval .
We have , the space of square integrable functions in . The orthonormal basis in is given by eigenfunctions of the Laplace operator that can be found as the solution of the Helmholtz equation
To reduce the boundary effects, we introduce the periodic constraint and the additional constraint of evenness (CSR2). The unique solution satisfying eq. 1 with these constraints is the basis
The eigenfunctions satisfy orthonormality
where is the Kronecker’s delta.
At each fixed , we estimate in the subspace spanned by up to the -th degree basis functions:
The least squares estimation (LSE) of in gives the -th degree cosine series representation (CSR):
where the coefficients are estimated as
Figure 1 illustrates the CSR at one particular voxel. The expansion degree is emiprically determined at 120 that gives 10% compression rate for fMRI time series with 1200 time points.
The CSR for multiple time series can be computed by setting up a matrix equation:
where is the matrix of CSR coefficients to be estimated. The least squares estimate of can be expressed as
2.3 Correlations between CSR
Consider the CSR of fMRI obtained from the first and second subject in a twin pair at the same voxel location . For instance, in Figure 1, voxel location is marked with red and blue circles for the first and second twin. With a slight abuse of notation, consider two CSR
are CSR coefficients. Since fMRI were normalized with respect to the mean, we have
be the variance of and , respectively. The correlation between and can be computed using the inner product on the CSR coefficients.
are vectors of cosine series representation coefficients normalized to the unit norm, i.e.,
Note is not the correlation of the original fMRI signals but their low frequency components. The expansion degree is truncated as follows. As , the correlation between the original rs-fMRI. In task-related twin fMRI studies, correlation between twins can be computed in a straightforward manner in the temporal domain because the timing of neuronal activation is comparable across subjects due to the exterior task. However, in the resting-state fMRI there is no external anchor that will ’lock’ brain activation of twin subjects across time. Because of this, we utilize cosine series representation of fMRI signals to compute correlation between twin subjects in the frequency domain instead of the temporal domain. Similar approaches were used in curtis.2005; ombao.2008, where frequency components of signals are correlated using coherence.
2.4 Heritability Index Filtration
The concept of heritability comes from genetics where it was used for decades to estimate the influence of genes on the phenotype of different traits. In the standard ACE model, the population-level variance of any trait is due to three components: genetic effects (heritability), common environment , events that happen to both twins, affecting them in the same way, and unique environment , events that occur to one twin but not the other, or events that affect either twin in a different way (Twins3). According to the ACE model, MZ twins are expected to share both 100% of their genes, and all of the shared environment. Any differences arising between them are random and due to mostly unique environmental effects. In the current problem setting we are more interested in comparing MZ and DZ twin pairs, rather than performing twin-wise comparison in MZ twin pairs. The correlation between MZ twins is then estimated as . DZ twins also share the common environment but share on average only of genetic information. The correlation between DZ twins is given by . More generally, we assume that not only twins in a twin pair share the same common environment, but all MZ and DZ twin pairs share it. Then, the heritability index (HI) is given by the Falconer’s formula (HI) as
and quantifies the influence of genetic effects.
We propose a framework aiming to identify a collection of the most genetically affected regions. We employ the AAL atlas with regions to compute pairwise twin correlation at region level by averaging across voxels in each parcellation. Since we cannot average correlations by taking the arithmetic mean without biasing, we transform correlations computed at the voxel level using the Fisher -transform first:
Then, Fisher transformed correlations are averaged and back projected via the inverse transform (Fisher2). At the -th region of AAL parcellation, we have modified Falconer’s formula
where is the average pairwise correlation between twins of a given zygosity across all voxels in the -th region. Here the correlations refer to to correlations between CSR coefficients within twins using eq. 8.
We expect that not all voxels in brain regions of interest possess high heritability and contribute equally to the classification accuracy. Thus, we filter out insignificant voxels with and average correlations only across significant voxels with . As the result of this heritability index-based voxel filtration, we suppress AAL regions that do not have significant voxels before classification. Note the heritability index filtration is only used in the training set in the subsequent classification.
2.5 Twin Classification with Artificial Neural Network
Given paired twin images represented as a vector of region-level correlations between CSR coefficients of the original twin fMRI, we use ANN to identify if they belong to a pair of MZ or DZ twins. Since we classify the relationship between twin fMRI instead of the original images themselves, this is not a traditional binary classification problem often done in brain imaging.
We use a two-layer feed-forward ANN, with 200 sigmoid hidden and 1 softmax output neurons to classify the vectors of average correlations. We train the network to minimize cross-entropy loss function, using the scaled conjugate gradient backpropagation method (SCGD). Under the assumption of possible dependencies between AAL regions and that AAL regions contribute differently to the classification accuracy, we introduce -regularization term to the loss function to increase the sparsity of the model
where is the class label for the -th twin pair (’1’ for MZ twin pair, ’0’ for DZ twin pair), is the predicted class label for the -th twin pair, is the number of twin pairs in the training set, is a vector of weights corresponding to a -th AAL region, is the number of regions of interest after heritability index filtration. In practice, solving a binary classification problem with artificial neural nets is equivalent to solving a regression problem with subsequent application of thresholding rule:
where is numerical output of the neural network for the -th pair of twins, and is the discrimination threshold.
Given a binary classification problem, we set a pair of MZ twin fMRI belonging to the positive class () and a pair of DZ twin fMRI belonging to the negative class (). We use the holdout method to split dataset randomly into training (70% of the data), validation (15%) and test (15%) subsets. The validation dataset is used to avoid overtraining of the model, and to fine tune hyperparameters of the neural network, e.g., the number of hidden neurons and discrimination threshold . Splitting data into training, validation and test subsets in the proportion of 70:15:15 is considered as the gold standard (Datasplit). When splitting dataset into subsets, we do not preserve class proportions in each subset. Thus, the performance of the model is highly dependent on the ratio of classes in subsets. To decrease the influence of random split on the results, we train 1000 independently initialized models and average across them. We employ classification accuracy, false-positive rate (FPR) and false-negative rate (FNR) to measure the performance of the proposed framework and report them in Table 3.
2.6 Boosting Classification Accuracy with Hill Climbing
Performing voxel-level filtration with respect to the heritability index followed by the -regularization modelling we obtain the sparse classification model. We further implement Hill Climbing variable selection procedure (considering each AAL parcellation as a variable) to infer on the contribution of different regions to the classification accuracy of the model. Hill Climbing is an iterative optimization technique that attempts to find a better solution by incrementally changing a single element of the solution. If the change produces a better solution, an incremental change is made to the new solution, repeating until no further improvements can be found (HillClimbing). With the application to variable selection, we implement Hill Climbing procedure as follows.
We start with the empty variable space and a pool of candidate variables. At each iteration, we test candidate variables by picking one variable at a time, adding it to the model and estimating the performance of the model. When all candidate variables are tested, the variable that provides the best performance, with respect to classification accuracy, FPR and FNR, is removed from the pool of candidates and added to the variable space of the model. Variables are ranked with respect to classification accuracy. If two or more variables have the same accuracy, the one with the smallest sum of FPR and FNR is selected. We continue the process iteratively until all variables are added to the model. We rank AAL regions with respect to the classification accuracy they provide to the model based on iterations, at which corresponding variables are added to the variable space. We assume that the higher the region is ranked by Hill Climbing procedure the more the region is affected by the genetic effect, in other words, the higher its heritability.
3 Simulation Study
|Study||Accuracy (%)||FPR (%)||FNR (%)|
We validate the proposed twin classification framework using a dataset generated with known ground truth. To perform the simulation, we first generate ground truth data that represents the underlying resting-state functional MRI signals from distinct brain regions of interest using degree 5 CSR. The CSR coefficients are
We generate twin data in -th region using the following statistical model
where is the vector of CSR coefficients of the first twin in the -th pair, is the vector of CSR coefficients of the second twin in the -th pair, is a twin-level noise distributed as and is an individual-level noise distributed as . After the twin data is generated, we process it according to the proposed framework: compute correlations between twins, train ANN to classify zygosity of twin pairs, and estimate the extent by which regions of interest are affected by the genetic effects using Hill Climbing variable selection procedure.
We perform three different simulation studies using the same ground truth coefficients and individual-level variance but different twin-level variances, depending on the zygosity of twin pairs. In all three simulations, we generate 50 MZ pairs and 50 DZ pairs. To obtain stable results we repeated the simulation 1000 times and the average results are reported. The simulation results are processed by our classification framework and averaged results are summarized in Table 2 and Figure 2.
Study 1: No twin difference
In the first study, we test if the method is detecting any false positive when there is no twin difference. By letting the DZ-variability equals the MZ-variability
we are generating twin data without any DZ- and MZ-twin difference. The classification accuracy is indicating we are not falsely classifying the groups.
Study 2: Twin difference
In the second study, we force all 5 regions to be equally highly heritable. We achieve this by simulating MZ twin data with smaller twin-level variance compared to DZ twin data:
The classification accuracy is indicating we are classifying zygosity.
Study 3: Differential twin difference
In the third study, we force 5 regions to have gradually decreasing heritability:
where is a sequence of gradually decreasing numbers that control heritability of each -th region of interest. As seen on Figure 2, ROI rankings provided by the Hill Climbing variable selection procedure correspond to the induced gradually decreasing heritability of regions (in contrast to cases of equally small (Study 1) or equally high (Study 2) heritability of regions). The classification accuracy in Study 3 is . When utilizing Hill Climbing procedure, the classification accuracy increased to indicating we are gaining advantage of employing variable selection in case of unequal distribution of heritability across regions.
We applied the proposed framework to the resting state fMRI of 208 twin pairs from the HCP dataset. We used the first session with oblique axial acquisition in a left-to-right direction as well as the second session with the right-to-left direction. We show that our results and findings are consistently observed in two datasets.
|Phase encoding||Method||Accuracy (%)||FPR (%)||FNR (%)|
|LR||Base model w/o CSR||54.15(9.24)||74.67(18.74)||28.36(17.09)|
|RL||Base model w/o CSR||52.98(10.55)||75.39(18.52)||29.73(16.87)|
|LR||Sparse model + Hill Climbing||94.19(3.53)||9.54(6.87)||3.69(3.47)|
|RL||Sparse model + Hill Climbing||93.33(3.59)||10.25(7.04)||4.39(3.53)|
According to the proposed framework, we performed dimensionality reduction of the original twin fMRI data using the 120-degree cosine series representation at each voxel. First, we mean-centered the fMRI signal at each voxel, computed 120-dimensional vectors of CSR coefficients and normalized them to the unit norm. In the next step, we obtained 116 feature vector representation of twin data, computing correlations at the voxel level between twins and averaging them in each brain parcellation. The resulting vector representation of twin data was used to train an artificial neural network to classify zygosity of twins. We trained a base model according to the description of the training procedure presented in Section 2.5. Its performance on the test subsets of the two datasets is reported in Table 3. For the dataset obtained during left-to-right scanning, we achieved classification accuracy with false-positive rate and false-negative rate. For the dataset obtained during right-to-left scanning, we achieved classification accuracy with false-positive rate and false-negative rate. Table 3 displays the results of twin classification employing region-level correlation between original twin fMRI as variables of the input data to the ANN.
We assume that brain regions of interest do not contribute equally to the differentiation of twin pairs with respect to the zygosity type. Thus, we aimed to reduce the dimensionality of the variable space from 116 brain ROIs. We achieved this by applying heritability index-based filtration at the voxel level followed by sparse modelling of the variable space by introducing an -regularization term to the loss function of the classification model. Both sparse modelling and heritability-based voxel filtration were performed only on the training data subsets and their results were used afterwards on the validation and test subsets. We trained a new sparse model and reported its accuracy on two datasets in Table 3. Processing functional MRI obtained with LR phase encoding with the proposed framework we achieved classification accuracy, with FPR and FNR; and using fMRI from RL phase encoding in the same framework we achieved overall classification accuracy, with FPR and FNR.
Figure 4 illustrates the results of the Hill Climbing variable selection procedure (here, variables represent AAL parcellations) accumulated across 1000 independently initialized models for left-to-right and right-to-left oblique axial data acquisition, respectively. Here, we show how often the Hill Climbing procedure selects a certain variable at a given iteration. Based on the results obtained during the variable selection procedure, we infer on the extent of how much these ROI are contributing to the zygosity classification. We quantitatively estimated the contribution of -th AAL parcellation using the following expression:
where is the number of times when the brain region has been added to the variable space at the -th iteration of the Hill Climbing procedure, and is the dimensionality of the variable space. Regions that have been selected by the Hill Climbing procedure at least once across all models are marked with bold font, and the most important -th percentile regions are marked in red in Figure 4.
Employing Hill Climbing variable selection procedure not only allowed us to estimate the importance of AAL parcellations with respect to the classification accuracy, it also provided us with a tool to find the optimal variable space with the highest possible classification accuracy. Figure 3 illustrates the box-plot (Boxplot) of the dimensionality of variable space after heritability index-based voxel filtration and sparse modelling and after Hill Climbing variable selection for both left-to-right and right-to-left oblique axial data acquisition. The median dimensionality of variable space after Hill Climbing variable selection procedure is 4 variables for both datasets. The average performance of the classification model for left-to-right acquired data is classification accuracy with FPR and FNR, and for right-to-left acquired data is classification accuracy with FPR and FNR.
|Left-to-right Phase Encoding|
|Region||Label||Frequency||Accuracy (%)||FPR (%)||FPR (%)|
|Left gyrus rectus||GRG||83/1000||86.36 (4.06)||23.04 (10.95)||8.68 (5.16)|
|Left middle frontal gyrus, lateral part||F2G||70/1000||85.03 (4.93)||22.43 (9.73)||10.66 (5.59)|
|Left superior frontal gyrus, dorsolateral||F1G||52/1000||81.63 (4.97)||20.74 (8.58)||16.92 (7.22)|
|Right supramarginal gyrus||GSMD||50/1000||85.83 (6.76)||22.39 (10.6)||8.87 (6.32)|
|Right area triangularis||F3TD||45/1000||86.48 (3.82)||18.98 (9.36)||9.98 (6.39)|
|Left middle temporal gyrus||T2G||40/1000||86.67 (5.44)||19.48 (11.79)||10.22 (4.97)|
|Right superior temporal pole||T1AD||39/1000||86.09 (6.56)||20.43 (13.14)||10.61 (7.68)|
|Left inferior occipital||O3G||37/1000||87.69 (5.56)||19.5 (13.02)||8.46 (5.01)|
|Left supramarginal gyrus||GSMG||36/1000||88.41 (4.76)||20.98 (12.27)||7.22 (4.36)|
|Left transverse temporal gyri||HESCHLG||36/1000||85.56 (7.1)||22.17 (13.67)||9.85 (5.3)|
|Right inferior occipital||O3D||34/1000||88.32 (4.78)||18.14 (11.88)||8.02 (5.13)|
|Left angular gyrus||GAG||34/1000||86.94 (4.9)||20.03 (8.23)||7.71 (6.7)|
|Right orbital part of inferior frontal gyrus||F3OD||30/1000||85.42 (5.24)||20.48 (11.06)||11.06 (6.19)|
|Left calcarine sulcus||V1G||28/1000||87.32 (5.11)||22.11 (11.69)||5.88 (4.75)|
|Right transverse temporal gyri||HESCHLD||27/1000||84.77 (7.22)||21.91 (12.99)||11.24 (7.92)|
|Right middle frontal gyrus, lateral part||F2D||25/1000||86.64 (4.3)||20.84 (11.6)||9.46 (4.23)|
|Right-to-left Phase Encoding|
|Region||Label||Frequency||Accuracy (%)||FPR (%)||FPR (%)|
|Right superior temporal pole||T1AD||79/1000||85.4 (6.94)||20.72 (11.91)||11.24 (7.31)|
|Right anterior cingulate gyrus||CIAD||58/1000||84.82 (5.44)||26.05 (9.41)||9.2 (6.8)|
|Left superior temporal pole||T1AG||54/1000||85.57 (6.14)||20.55 (8.69)||11.07 (9.1)|
|Right supramarginal gyrus||GSMD||52/1000||86.72 (3.72)||21.35 (7.53)||8.84 (6.29)|
|Left middle cingulate||CINMG||49/1000||84.09 (4.56)||21.04 (12.23)||13.09 (6.25)|
|Left middle temporal gyrus||T2G||43/1000||85.27 (5.46)||24.81 (8.76)||9.19 (7.64)|
|Left supramarginal gyrus||GSMG||41/1000||87.03 (4.3)||20.63 (8.7)||8.76 (4.35)|
|Left middle frontal gyrus, lateral part||F2G||40/1000||79.4 (7.24)||27.19 (10.21)||16.98 (7.91)|
|Left postcentral gyrus||PAG||38/1000||83.96 (9.2)||24.2 (12.74)||11.55 (9.51)|
|Left transverse temporal gyri||HESCHLG||38/1000||87.77 (6.08)||20.61 (10.57)||7.62 (5.85)|
|Right postcentral gyrus||PAD||34/1000||86.93 (3.31)||19.28 (8.9)||9.66 (3.95)|
|Left angular gyrus||GAG||31/1000||88.06 (3.4)||16.9 (7.77)||9.21 (4.78)|
|Right superior temporal gyrus||T1D||25/1000||87.17 (4.39)||22.34 (10.05)||7.61 (3.27)|
|Right fusiform gyrus||FUSID||23/1000||84.54 (5.14)||24.06 (10.75)||10.73 (7.87)|
|Right thalamus||THAD||23/1000||84.96 (5.46)||27.34 (12.66)||8.27 (5.4)|
|Right amygdala||AMYGDD||22/1000||86.02 (5.59)||16.74 (7.81)||12.46 (7.21)|
|Left superior temporal gyrus||T1G||22/1000||88.63 (4.74)||14.62 (7.25)||9.58 (6.62)|
For the left-to-right phase encoding, the following areas are identified as the most important AAL regions with respect to the overall contribution to the classification accuracy: Left middle frontal gyrus, lateral part, Left superior frontal gyrus, dorsolateral, Left gyrus rectus, Left middle temporal gyrus, Right supramarginal gyrus, Right superior temporal pole, Left supramarginal gyrus, Right inferior occipital, Left angular gyrus, Left inferior occipital, Right area triangularis, Left transverse temporal gyri, Left calcarine sulcus, Right calcarine sulcus, Right transverse temporal gyri, Right superior frontal gyrus, dorsolateral, Right middle temporal gyrus.
For the right-to-left phase encoding, the following areas are identified as the most important AAL regions with respect to the overall contribution to the classification accuracy: Left middle frontal gyrus, lateral part, Right superior temporal pole, Left superior temporal pole, Right anterior cingulate gyrus, Left middle temporal gyrus, Left middle cingulate, Right supramarginal gyrus, Left supramarginal gyrus, Left postcentral gyrus, Left transverse temporal gyri, Right postcentral gyrus, Right superior temporal gyrus, Right gyrus rectus, Left angular gyrus, Right amygdala, Right middle temporal gyrus, Right thalamus.
Based on the overlap of two encoding schemes, we infer that the most heritable brain regions using consistent findings in two datasets are: Left middle frontal gyrus, lateral part, Left supramarginal gyrus, Right supramarginal gyrus, Left angular gyrus Left transverse temporal gyri, Right superior temporal pole, Left middle temporal gyrus, Right middle temporal gyrus.
5 Conclusion and Discussion
In this paper, we address the problem of classifying the zygosity of a pair of twin fMRI. This is a more complex problem than the usual classification problem of labeling each image into distinct classes. Here, we are interested in learning if the relationship between pairs of images is associated with the zygosity of twins. We solved the stated problem within the proposed pipeline through sequential preprocessing and classification of the data.
There were two practical advantages for using CSR as a new feature representation of twin fMRI. First of all, correlating between CSR coefficients allows to correlating signal in the frequency domain. Furthermore, representing original fMRI signals as a linear combination of 120 cosine basis functions serves not only as a dimensionality reduction technique but also as a means of denoising high frequency noise. To emphasize the importance of the new feature representation for the performance of the proposed framework, we classified twin fMRI without cosine series representation and obtained and for left-to-right and right-to-left phase encoding, respectively (Table 3, first two rows). Correlation between twin fMRI in the frequency domain employing coherence of twin CSR allowed us to gain the classification accuracy of for LR phase encoding and for RL phase encoding, the increase of and , respectively.
We applied our framework to two datasets from the Human Connectome Project, acquired in two sequential scanning runs. We were able to achieve an average classification accuracy of in determining the zygosity of twins for the left-to-right acquired data, and classification accuracy for the right-to-left acquired data.
Regions with significant genetic influence
In order to infer on the conformity of our findings in the resting-state functional MRI with regard to the most genetically affected brain regions, we examined a cohort of twin studies. Cannon3228 found that genetic influences were isolated primarily to polar, inferior, and dorsolateral prefrontal brain areas, and also in the frontal region. Despite the present discrepancy in the identified brain regions, there is a certain conformity of our findings with results reported in a number of twin studies of both task-based and resting-state functional MRI (BLOKLAND200870; Blokland10882; Koten1737; TwinStudy2010; PARK20121132; Gao11288; SINCLAIR2015243; Yang2016).
MATTHEWS2007223 found the first functional imaging evidence that dorsal anterior cingulate cortex activation during interference processing is significantly influenced by genes. Olli2008 demonstrated high heritability of medial and dorsolateral prefrontal cortex in a study on genetic influence of the verbal learning and verbal memory in twins. BLOKLAND200870; Blokland10882 found that the inferior, middle, and superior frontal gyri, left supplementary motor area, precentral and postcentral gyri, middle cingulate cortex, superior medial gyrus, angular gyrus, superior parietal lobule, including precuneus, and superior occipital gyri are genetically affected in twins. Koten1737 investigated genetic influences and observed significant genetic influences on brain activation in visual cortex, temporo-parietal and frontal areas, and anterior cingulate cortex. PARK20121132 found neural activity in the left visual cortex and left motor cortex were significantly heritable. SINCLAIR2015243 found that 47 out of 116 AAL regions are significantly heritable which overlap with most of our regions. Based on the examined twin literature, it is apparent that there is a certain consistency in the reports about the heritability of brain function in different regions of interest, however some disparities are also present that suggest that genetic influences may vary with task paradigm and brain region.
In this paper we report that the most genetically affected brain regions, as measured by their contribution to the classification accuracy of the zygosity type of twins, are mainly located in the temporo-parietal and frontal brain regions with the dominance of heritable regions in the left hemisphere. We identify these regions performing independent model training on two datasets of the Human Connectome Project. Though there is an evident partial accordance between regions with the highest genetic influence of different phase encoding datasets, there are some regions that were identified as the most genetically affected only for one dataset. This phenomenon requires a thorough investigation and is a subject for the future research.
AAL parcellation template
In the proposed framework, we compute pairwise correlation between twin subjects and then average them within predefined AAL parcellation. However, from the review of several structural and functional voxel-based twin studies, it is apparent that genetic effects may carry across anatomical boundaries (Blokland10882; joshi2011contribution; VANSOELEN20123871). Therefore, voxel-based approaches may have preference in imaging genetic studies over ROI approaches that average measurements across brain parcels. On the other hand, the drawback of voxel-wise approaches in general is an extensive computational load resulting from the computation of pairwise correlations at voxel level. In the proposed framework, we overcome this drawback by utilizing 120-degree cosine series representation that not only performs denoising of the original fMRI signals and increase classification accuracy, but also provides a means to drastically decrease computational time by compact signal representation. We take the advantage of both methods in increasing the classification accuracy. Although we demonstrated relatively high classification accuracy with a provided solution, there is an obvious need for better parcellation method that balances the trade-off between pure voxel-wise and anatomical template-based approaches, yet we leave the search of such method as a future study.
Adaptation for Deep Learning
The twin classification framework proposed in this paper is aimed to provide an evidence of the possibility to identify zygosity of a twin pair given only resting-state functional MRI. Though our framework shows high classification accuracy, the exact practical implementation described in this paper, including cosine series representation for dimensionality reduction, heritability index-based voxel filtration, -regularization and Hill Climbing variable selection for sparse modelling, and two-layer artificial neural network for classification, is not necessary the optimal approach and is a subject to further improvements.
The aspect of considering pairs of input data employed in the proposed framework is akin to Siamese neural networks (Siamese1; Siamese2). As a result, our framework can be easily adapted to be implemented as Siamese neural networks with deep architecture. This can be achieved by using Siamese neural network to learn the representation of twin data, instead of computing it explicitly in a model-based approach through HI filtration of pairwise correlations between twins at the voxel level sparse modelling that is aimed to preserve only significant regions of interest. Utilizing deep neural networks, e.g., convolutional neural networks or deep belief networks, as building blocks of Siamese neural network may provide exceptionally high classification accuracy by means of extraction of high-level abstract features from high-level layers, though at the cost of decreased interpretability of these discriminative features. To improve interpretability of output of deep Siamese neural network, for instance to infer significance of the genetic control on certain brain regions of interest, a more sophisticated design of the twin classification framework is required and one can consider using discriminative feature descriptors in combination with deep neural networks (Cimpoi2016; Song2017). This is left as a future study.
To infer the extent of genetic control on AAL parcels, we use variable selection applying the Hill Climbing algorithm – a greedy search algorithm that considers regions of interest as independent variables and tests one variable at a time. Recent studies of both resting-state and task-related fMRI on the functional connectivity of brain regions has revealed that activity in some regions may have strong coherence (TwinStudy2015; Yang2016). To get a more accurate inference on the optimal variable space, one may consider applying the concept of ”connectivity of regions” and performing a group variable selection, i.e., combine variables in groups and test each group as a single instance (HillClimbing). Additionally, this will reduce computational load of variable selection procedure, whose computational complexity for brute-force approaches that test all possible combinations of variables is (Yang2016). Other type of variable selection methods are left as a future study.
This work was supported by NIH Brain Initiative grant R01 EB022856 and Clinical and Translational Science Award (CTSA) program through the National Center for Advancing Translational Sciences (NCATS) grant UL1TR002373. We would like to thank Guorong Wu of University of North Carolina, Chapel Hill and Hernando Ombao of King Abdullah University of Science and Technology for valuable discussions and supports.