Hill Climbing Optimized Twin Classification Using Resting-State Functional MRI

Hill Climbing Optimized Twin Classification Using Resting-State Functional MRI

Andrey Gritsenko gritsenko@wisc.edu Martin A. Lindquist Gregory R. Kirk Moo K. Chung mkchung@wisc.edu Department of Biostatistics and Medical Informatics, University of Wisconsin, Madison, WI, USA Department of Biostatistics, Johns Hopkins University, Baltimore, MD, USA Waisman Laboratory for Brain Imaging and Behavior, University of Wisconsin, Madison, WI, USA

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.

Twins, heritability, resting-state fMRI, sparse model, Hill Climbing, Human Connectome Project
journal: NeuroImage

1 Introduction

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.

Figure 1: Pipeline of the proposed framework. The original fMRI signals are displayed as dashed lines. The 120-degree cosine series representations are displayed as the solid lines. Correlation between CSRs are computed at each voxel. We incorporate the filtration of voxels at the group level using the modified heritability index in the training set. We use AAL brain atlas to parcellate the voxel-level correlations into 116 region-level correlations. The resulting vector of 116 correlations feed into a two-layer feedforward neural network to determine zygosity of twin pair automatically.

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

2.1 Dataset

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
  MZ 131 56/75
  DZ 77 30/47
MZ: monozygotic twins, DZ: dizygotic twins. Age is displayed in years as MeanSD.
Table 1: Demographic data of twin groups

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 (%)
1 w/o HC 48.63(12.37) 50.64(24.24) 49.98(24.04)
 2 w/o HC 79.79(10.66) 24.65(14.66) 13.51(14.43)
 3 w/o HC 81.88(10.43) 23.02(16.84) 12.98(14.40)
 3 with HC 87.65(8.11) 16.93(13.98) 7.24(9.78)
FPR: false-positive rate, FNR: false-negative rate and overall classification accuracy are provided as MeanSD. For Study 3 we present classification performance both with (’3 with HC’) and without (’3 w/o HC’) employing Hill Climbing procedure. For Study 1,2 we present results only once since no statistical difference has been observed when utilizing Hill Climbing.
Table 2: Simulation study results
Figure 2: Results of Hill Climbing variable selection for three simulation studies. For each region of interest , the value at the intersection of -th column and -th row represents how often this regions has been selected at the -th iteration of the Hill Climbing variable selection procedure in the corresponding study.

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.

4 Results

HCP Data

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 Base model 79.93(7.59) 36.98(17.29) 9.99(6.79)
  RL Base model 80.54(7.37) 37.11(17.01) 9.95(6.86)
  LR Sparse model 89.73(3.14) 16.03(8.82) 6.94(4.12)
  RL Sparse model 89.31(2.99) 17.56(8.65) 6.92(4.25)
  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)
Phase encoding represents the way the data was acquired. LR: left-to-right oblique axial acquisition, RL: right-to-left oblique axial acquisition, simulation: simulated data as described in Section 3. FPR: false-positive rate, FNR: false-negative rate. Overall classification accuracy, false-positive rate and false-negative rate are provided for the test subset in percents as MeanSD. The first row corresponds to the results obtained when region-level correlations between twin fMRI has been employed as variables of the input data fed to the classification model. Neither HI-filtration at the voxel level, nor sparse modeling has been used in this experiment.
Table 3: Performance of the proposed classification pipeline at different stages

Base Model

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.

Sparse Model

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.

Hill Climbing

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.

Figure 3: Dimensionality of variable spaces of the model before (first and third boxes) and after (second and forth boxes) Hill Climbing optimization. Oblique axial acquisition is encoded LR for left-to-right acquired data, and RL for right-to-left acquired data. Values are gathered across 1000 independently initialized models.
Figure 4: Results of Hill Climbing variable selection. Results are accumulated over 1000 independently initiated models for left-to-right (A) and right-to-left (B) phase encoding. For each region of interest , the value at the intersection of -th column and -th row represents how often this regions has been selected at the -th iteration of the Hill Climbing variable selection procedure. Regions in bold font appear at least once in the variable space of the sparse model, i.e. after heritability index filtration and -regularization. Both HI filtration and regularization serve to increase the sparsity of the model. As a result, the dimensionality of the variable space is much less than the original 116-dimensional variable space, and in most cases it consists of 10 AAL regions. Red colored regions represent the variable space of the most important AAL regions (above -th percentile). Importance of a region is estimated using eq. 18. Therefore, we can equally state that these regions are the most genetically affected regions.
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)
FPR: false-positive rate, FNR: false-negative rate. Overall classification accuracy, false-positive rate and false-negative rate are provided for the test subset in percents as MeanSD. They present performance of the model at the first iteration of the Hill Climbing procedure, i.e. when the variable space contains only one variable. Bold font highlights the best performance. Region names are given in the first column according to the AAL brain atlas specification (http://neuro.compute.dtu.dk/services/brededatabase/index_roi_tzouriomazoyer.html), while the second column displays regions’ labels as they appear on Figure 5. Regions are sorted according to the frequency (across 1000 independently initialized models) that they were selected at the first iteration of Hill Climbing procedure in the descending order. Cumulatively, presented regions are selected at the first iteration of Hill Climbing procedure for 666 out of 1000 models for left-to-right phase encoding (top) and 672 out of 1000 models for right-to-left phase encoding (bottom).
Table 4: Most frequent AAL regions selected at the first iteration of the Hill Climbing procedure
Figure 5: Most important AAL regions in the LR (A) and RL (B) phase encodings (reference to Figure 4). Regions of interest are ranked according to their importance computed using eq. 18 and only regions wit h high importance , -th percentile, are colored and correspond to regions marked with red font on Figure 4.
Figure 6: Most important AAL regions in the LR and RL phase encoding. We combined AAL regions in three groups: 1) regions that were detected as the most genetically influenced based on the analysis of LR phase encoding only (yellow color), 2) regions that were detected as the most genetically influenced based on the analysis of RL phase encoding only (blue color), and 3) regions that were detected as the most genetically influenced both based on the analysis of LR and RL phase encoding (green color).
Figure 7: AAL regions most frequently selected at the first iteration of Hill Climbing in the LR (A) and RL (B) phase encodings. Only AAL regions that are selected more frequently than -th percentile are colored. Refer to Table 4 for a detailed information on the contribution of each parcellation to the performance of classification model.

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.

Variable selection

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.


Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description