Show
cerebral cortex, convolutional neural network, gyri, sulci, wavelet entropy One of the most prominent feature of human brain lies in the high-level convolution of cerebral cortex. During brain development, the cerebral cortex undergoes substantial folding, leading to intricate arrangements of convex gyri and concave sulci (Rakic 1988; Zilles et al. 1988; Armstrong et al. 1995; White et al. 2010). Over the past decades, neuroscientific communities have proposed various hypotheses concerning gyrification, such as differential laminar growth (Richman et al. 1975; Cartwright 2002; Tallinen et al. 2016; Razavi et al. 2017), genetic regulation (Beck et al. 1995; Rakic 2004; Barkovich et al. 2012; Stahl et al. 2013; Sun and Hevner 2014), and axonal tension or pulling (Van Essen 1997; Hilgetag and Barbas 2005; Nie et al. 2012; Zhang et al. 2017; Ge et al. 2018). Although controversy regarding gyrification mechanisms continues, there is a consensus that cortical convolution allows for a dramatic increase of cerebral size and compact wiring that enhances efficient information processing (White et al. 2010; Zilles et al. 2013). Recently, emerging studies have reported that gyri and sulci exhibit distinct patterns from the perspectives of anatomy (Fischl and Dale 2000; Smart et al. 2002), morphology (Magnotta et al. 1999), structural connectivity (Nie et al. 2012; Deng et al. 2014), and function (Hilgetag and Barbas 2005; Jiang et al. 2015). For example, microcosmic studies at the cellular level revealed that the subventricular zone (a secondary proliferative zone containing neural progenitor cells) during embryonic life is thicker in areas underlying gyrus formation and thinner in areas underlying sulcus formation (Smart et al. 2002; Kriegstein et al. 2006). A similar finding in the neocortex of both developing infants (Li et al. 2015) and healthy adults (Fischl and Dale 2000; Hilgetag and Barbas 2005) is that cortical thickness is significantly thicker at the crown of gyri and normally thinner in the depth of sulci. In addition, under the mechanical impact of cortical folding, dendritic shapes in deep gyral and sulcal layers are different (Hilgetag and Barbas 2005). Somata and arbors of pyramidal neurons residing in the deep cortical layers of gyri are clearly stretched radially, while somata in the deep sulcal layers are compressed tangentially and their apical dendrites are wavy (Hilgetag and Barbas 2005). Mechanically-induced variations of cortical morphology result in different modes of neuronal functioning in gyri and sulci (Hilgetag and Barbas 2005). For instance, the elongation of the gyral dendrites results in a greater attenuation of excitatory postsynaptic potentials and a potential failure to produce action potentials (Hilgetag and Barbas 2005). What’s more, gyri and sulci have different morphologic traits during aging process. Gyri become more sharply and steeply curved while the sulci become more flattened and less curved (Magnotta et al. 1999). Such morphological feature become abnormal in adolescents with schizophrenia spectrum disorders. Compared with healthy controls, adolescents with schizophrenia disorders have significantly more flattened curvature in the sulci and more steeped curvature in the gyri (White et al. 2003). More recently, investigations into the neuronal pathways demonstrated that axonal streamline fibers in the human brains tend to course horizontally around the sulcal surface but orient tangentially to the gyral surface (Takahashi et al. 2012; Budde and Annese 2013; Chen et al. 2013; Zhang et al. 2014). Similar observations suggest that the fiber terminations dominantly concentrate on gyri rather than sulci, which is evolutionarily preserved in the cortical architectures of primate brains (Nie et al. 2012; Chen et al. 2013; Li et al. 2015). Joint representation of multimodal neuroimaging data, including diffusion tensor imaging (DTI; Basser and Jones 2002; Beaulieu 2002), and fMRI (Biswal et al. 1995; Friston 2009), found that both structural and functional connectivity patterns are strong among gyri, weak among sulci, and moderate between gyri and sulci (Deng et al. 2014). These evidences support a functional model of cerebral cortex that gyri are functional connection centers that exchange information among remote structurally connected gyri and neighboring sulci, while sulci communicate directly with their neighboring gyri and indirectly with other cortical regions through gyri (Deng et al. 2014). Given the close relationship between structural substrates and neural activation patterns (Hagmann et al. 2008), sulco-gyral differences are also elucidated from the functional perspectives, such as functional network heterogeneity (Jiang, Li et al. 2018), functional signal reconstruction residuals (Jiang, Zhao et al. 2018), and graph-theoretic characteristics (Liu et al. 2017). No matter under temporal stationary or dynamics, task-based heterogeneous functional regions, i.e., the cortical regions that are activated during multiple tasks conditions, locate significantly more on gyral regions than on sulcal regions (Jiang et al. 2015; Jiang, Li et al. 2018). In addition, graph-theoretic analyses show that gyri have higher global efficiency of functional interaction than sulci in task-evoked paradigms (Liu et al. 2017). Taken together, a series of studies across multiple disciplines have consistently brought insights into fundamental differences between gyri and sulci. However, previous studies mainly unveiled the functional differences between gyri and sulci via interregional functional interactions (Deng et al. 2014; Liu et al. 2017) or functional network overlaps (Jiang et al. 2015; Jiang, Li et al. 2018). Relatively little progress has been made in characterizing the frequency-specific patterns of cortical neural activities. Neurosciences of brain rhythms suggest that high-frequency brain activity reflects local domains of cortical processing, while low-frequency oscillation rhythms synchronize across distributed brain regions and can bind specific assemblies together (von Stein and Sarnthein 2000; Buzsáki and Draguhn 2004; Canolty and Knight 2010; Siegel et al. 2012; Buzsáki et al. 2013). Combining the findings from neuroscience of brain rhythms (Buzsáki et al. 2013) with brain connectome (Deng et al. 2014), it is natural to propose hypotheses that there is frequency-specific difference between gyral and sulcal fMRI signals. To validate the proposed hypothesis, we adopted Human Connectome Project (HCP) grayordinate resting-state and task-based fMRI data, which have decent temporal and spatial resolution (Glasser et al. 2013; Van Essen et al. 2013). Notably, the decent temporal resolution of HCP grayordinate fMRI data (TR = 0.72) provides an prerequisite for our work, since high-speed fMRI can track oscillatory neural activity directly in higher frequencies than expected during human cognition (Lewis et al. 2016). Specifically, we characterized the frequency-specific feature of cortical fMRI signals using convolutional neural network (CNN) and wavelet entropy. We adopted CNN, one of the most powerful deep learning models, to model task-based fMRI signals and further to identify discriminative features of band pass filters in fMRI signals in data-driven fashion (Golik et al. 2015; Zhou et al. 2016; Huang et al. 2018). Notably, these convolutional filters in learned CNN models tend to activate class-specific cortical regions (i.e., gyri or sulci) and are nonlinear distributed in frequency band (Golik et al. 2015), which can be used to characterize the diversity and frequency of fMRI signals. As illustrated in Figure 1(a), we first labeled cortical vertices as gyral or sulcal type according to their average convexity. Then we designed a CNN-based classifier to differentiate cortical fMRI signals, where the fMRI signals and labels of cortical vertices are used as the input and output of the CNN, respectively. After that, we closely investigated the internal representations of learned CNN models in frequency domain. We also correlated such functional role characteristics (e.g., CNN predictions, Figure 1(d)) with known anatomic and structural connectivity attributes (Fischl 2012; Nie et al. 2012), such as fiber density (Fig. 1(b)) and cortical thickness (Fig. 1(c)) to examine possible underpinnings. Wavelet analysis is another alternative approach to analyzing temporal signals in frequency domain (Mallat 1989). Notably, wavelet entropy can quantify the degree of signal complexity of temporal signals (Coifman and Wickerhauser 1992; Sang et al. 2011). Specifically, we measured the wavelet entropy of cortical fMRI signals in different frequency scales and inspected if there are differences of wavelet representations between gyri and sulci. In general, above experiments allow us to assess the frequency-specific differences between gyri and sulci quantitatively, and provide novel insights into functional architecture of the cerebral cortex. An overview of structural, anatomical and functional differences between gyri and sulci. (a) Illustration of gyral/sulcal parcellation of cortical surfaces. According to average convexity, the cortical surfaces are labeled as gyri (red), sulci (blue) and in-between areas (green), respectively. (b) Spatial distributions of fiber density. (c) Spatial distributions of cortical thickness. (d) The CNN predictions. Areas in red and blue are predicted as gyri and sulci, respectively. The gray regions are discarded in CNN classification. Materials and MethodsData Acquisition and PreprocessingIn this work, we validated the proposed hypothesis using grayordinate resting-state and task-based fMRI data from the HCP 900 subjects (S900) data release (Van Essen et al. 2013). HCP provides publicly available multimodal neuroimaging data and behavioral measures of cognitive processes, which significantly facilitates exploitation of functional brain architecture (Uğurbil et al. 2013; Van Essen et al. 2013). The HCP resting-state fMRI data were acquired with eyes open with relaxed fixation on a projected bright cross-hair on a dark background. Each resting-state fMRI data consisted of 2 sessions of approximately 15 minutes each (1200 frames). In the HCP task-based fMRI data, 7 task paradigms were designed to activate a variety of cortical and subcortical networks within realistic time constraints, including emotion processing (176 frames), gambling (253 frames), language, motor (284 frames), relational processing (232 frames), social cognition (274 frames), and working memory (405 frames). See Barch et al. (2013) for more details about task designs. In the HCP S900 release, 897 healthy adult participants were recruited in data acquisition. These participants ranges from 22 to 35 years old. The acquisition parameters of fMRI scans are as follows: 90 × 104 matrix, 72 slices, TR = 0.72 s, TE = 33.1 ms, FOV = 220 mm, flip angle = 52°, BW = 2290 Hz/Px, in-plane FOV = 208 × 18 mm and 2.0 mm isotropic voxels. The grayordinate fMRI data were preprocessed according to HCP MR Data preprocessing pipelines version 3 (Glasser et al. 2013; Robinson et al. 2014), including spatial artifact/distortion removal, cortical surface generation, within-subject multimodal surface matching registration, and cross-subject alignment to standard space. In addition to minimal preprocessing pipeline, we normalized each vertices’ fMRI time-series to zero mean and unit variance, as we did in previous CNN autoencoder model (Huang et al. 2018). Parcellation of Cortical SurfacesThe HCP grayordinate fMRI data take advantage of geometry information of convoluted cortical sheet and reliably map the fMRI time series in volume space to cortical space, which allows for combined cortical surface and subcortical volume analyses (Glasser et al. 2013). Specifically, each grayordinate (i.e., spatial coordinates in gray matter) has associated geometric information and corresponding fMRI time series. A grayordinate can be annotated as gyri or sulci according to its geometric information while corresponding fMRI signals can serve as the input of following CNN classification tasks. Before fMRI signal classification, it is essential to establish an accurate ground truth database of gyral and sulcal annotations. To minimize the risk of wrong annotations, we focused on neighborhood around the crown of primary gyri and the fundi of primary sulci. Specifically, we segmented cortical surfaces according to vertices’ average convexity, since average convexity (i.e., “sulc” map in FreeSurfer) is the signed distance that a vertex moves during the inflation process and mainly depicts the primary folding patterns (Destrieux et al. 2010; Fischl 2012). As illustrated in Figure 1(a), cortical vertices were divided into gyri (red), sulci (blue) and in-between (green) of 3 types. The 30% vertices with the lowest convexity value were considered sulcal, and the 30% vertices with the highest convexity value were considered gyral. The remaining 40% vertices were regarded as transitional or intermixed areas between gyri and sulci that either do not belong to gyral or sulcal regions or we are not sure about. These in-between vertices were excluded in our subsequent classification task, which guarantees sufficient geodesic distances between gyral and sulcal regions. Although one half to two-thirds of the cortical surface are buried inside the sulci and in the lateral fissure of the brain (Zilles et al. 1997; Van Essen 2005), our labeling method would yield balanced dataset for gyri and sulci. Another geometric property often used to delineate gyri and sulci is mean curvature (Li et al. 2010), which lays special emphasis on the degree of local cortical folding. To figure out that whether folding or burying of sulci derive the sulco-gyral differences, we adopted a similar 3-bin parcellation procedure using mean curvature. In general, curvature-based findings were similar with convexity-based findings. Given the limited space, we mainly report these convexity-based results in this paper. To characterize class-specific CNN features for gyri and sulci, we encoded the labeled grayordinates by one-hot method. That is, sulcal and gyral examples were encoded by [1 0] and [0 1], respectively. Configuration of Convolutional Neural NetworkArchitecture of Convolutional Neural NetworkWe proposed a simple but effective CNN architecture (Fig. 2) to differentiate the labeled gyral/sulcal fMRI signals. It is noteworthy that we replaced conventional fully connected layer with global averaging pooling (GAP) layer, since introducing GAP has been an well-established approach for object detection and investigating discriminative features in the field of computer vision (Zhou et al. 2016). Introducing GAP layer can also dramatically reduce the number of trainable parameters, which greatly mitigates potential risk of overfitting (Lin et al. 2014). The non-linearity in the convolution layer is fulfilled by the rectified linear unit activation function (Nair and Hinton 2010). The dropout layer acts as a regularizer to prevent the proposed model from overfitting (Srivastava et al. 2014). Through softmax layer, the feature map of each filter channel connect to the 2 target units that are representative of the gyrus and sulcus labels. The proposed CNN model is implemented based on Keras (http://keras.io) with a backend of Theano (http://deeplearning.net/software/theano/). The architecture of the proposed CNN model. This model consists sequentially of 1 convolution layer, 1 max pooling layer, 1 GAP layer, 1 dropout layer, and 1 softmax layer. The rows in the convolution layer represent different convolution channels. Notably, there are correspondences between filters and weights, as illustrated by the line styles. In general, the proposed model has 2 compelling advantages. First, the convolutional network can extract features of neural activities in the temporal axis (Wang et al. 2017; Huang et al. 2018), since the input of the network is original fMRI time-series. These temporal features are nonlinearly distributed in frequency scale (Golik et al. 2015). Second, the proposed model contains only 1 convolution block for easy interpretability. As the average values of the feature map of each filter directly connects to softmax layer (Lin et al. 2014), there is one-to-one correspondence between convolution filters and softmax weights, which enables easy exploitation of class-specific features of gyral/sulcal fMRI signals. In brief, the proposed model provides us a promising approach to capture the time-frequency signatures and class-specific features of cortical fMRI signals (Lin et al. 2014; Zhou et al. 2016). Training FashionsGiven the inter-subject variability of fMRI data (Dubois and Adolphs 2016), we trained the proposed model at 2 levels, e.g., individual subject level and group subjects level. Since there are significant correlations between neighboring fMRI signals from the same scan, it is hard to split training or testing set at individual subject level. In addition, it is not appropriate to use other individuals’ fMRI signals as testing data due to the large inter-subject variability. Thus, at the individual level, we just trained CNN classifier without testing procedure for pattern discovery. In details, for a subject under a specific task, we trained a unique CNN model using all available individual fMRI signals from that subject under the corresponding task. In this way, each subject has 7 models corresponding to the 7 different tasks. These individual-level models were used to test the consistency of our findings across different tasks. Also, we conducted additional experiments to validate the CNN models at individual level. Specifically, we adopted resting-state fMRI data in 2 different sessions as an additional validation test. We trained the CNN model on the fMRI data of 1 session and tested the model on the fMRI data of the other session. At group level, we randomly divided all subjects into 2 disjoint subsets. In training procedure, groups of subjects were excluded from the training dataset, and the trained model was tested on these subjects. Notably, different training fashions could validate our hypotheses comprehensively. Generally, individual-level experimental results can validate consistency of our findings while group-level experiments can detect the commonly shared patterns from multi-subject fMRI signals. Hyperparameter SettingsThe labels of cortical vertices annotated in Section 2.2 are used as the output of the CNN, while their corresponding fMRI signals are used as input of CNN. One crucial issue concerning CNN models is hyperparameter settings. Based on our prior experiences (Huang et al. 2018), we empirically set the filter number, filter length, pooling stride, and batch size to 64, 32, 2, and 128, respectively. The dropout layer discards 20% hidden nodes randomly. The loss function is the categorical cross-entropy. We randomly initialize the trainable weights with unit normal distribution. The learning epoch is set as 100 at individual level and 10 at group level, which ensures convergence of training procedure. We use stochastic gradient descent (SGD) optimizer to train the proposed CNN model. Although SGD is intrinsically stochastic, the results are consistent in different runs of training CNN on a same dataset (Supplementary Fig. S1). To demonstrate the empirical rigor and validity of the experimental results, we checked the consistency of the experimental findings using different hyperparameter settings of filter number (64, 128) and filter length (32, 64) on 1 subset of WM fMRI data at group training level. Generally, the experimental results under different hyperparameter settings are consistent with previous findings (Supplementary Fig. S2). Perhaps there are better hyperparameters for our tasks, our current parameter setting works well enough to provide a qualified demonstration of sulco-gyral differences in frequency domain. Exploration of Learned CNN ModelsGiven easy interpretability of the proposed CNN model, we interpret the internal representations of learned CNN models from 3 aspects, including visualization of the CNN predictions, class-specific division of convolution filters and time-frequency analysis of the filters. First, we mapped the CNN predictions back to the cortical surface and examined their spatial distribution patterns within a neuroanatomical context. These spatial patterns will informatively exhibit the geometrical distributions of the classification accuracy across the entire cortical surface. As mentioned in Section 2.3, each filter in convolution layer correspond to the weights in softmax layer, each of which indicates the potential target class of these filters. Specifically, the filters can be divided into gyral, sulcal or uncertain types according to the sign of softmax weights, as formulated in Equation (1). K={gyral,ifwG>0,wS<0sulcal,ifwG〈0,wS〉0uncertain,ifwG∗wS>0 (1)In Equation (1), K denote the type of filter, while wG and wS denotes the weight connecting to gyral unit and sulcal unit, respectively. Literally, given a filter, if its gyral weight wG is positive while its sulcal weight wS is negative, it belongs to gyral type. While wG is negative and ws is positive, this filter is considered as sulcal type. If a given filter’s weight wG and ws has some sign, this filter is considered as uncertain type. After division of convolution filters, we can quantify the histogram of these class-specific filters. Finally, we performed time-frequency analysis of learned filters and examined the frequency-specific characterizations. Specifically, we aggregated all filters from all learned models of a specific task, calculated their power spectrum using fast Fourier transformation and examined power spectrum differences of gyral and sulcal filters. There are 3 primary reasons for time-frequency analysis. First, human brains are likely to be a multi-frequency oscillation system (Deco et al. 2017). Second, fast fMRI acquisition techniques has significantly increased fMRI bandwidths, which make it possible to scan human brains with higher sampling rate (Lewis et al. 2016). The time repetition of HCP fMRI data equals to 0.72 s with sampling frequency of 1.39 Hz, which overlap with multiple frequency bands according to earlier electrophysiological studies (Penttonen and Buzsáki 2003; Buzsáki et al. 2013), including slow-5 (0.01–0.027 Hz), slow-4 (0.027–0.073 Hz), slow-3 (0.073–0.198 Hz), slow-2 (0.198–0.5 Hz), and slow-1 (0.5–0.69 Hz). Third, CNN-based models are very powerful in learning class-specific features which are nonlinearly distributed in frequency scale (Zhou et al. 2016; Huang et al. 2018). Underpinnings of Frequency-Specific Differences Between Gyri and SulciThere are known anatomical and structural differences between gyri and sulci, such as cortical thickness (Fischl and Dale 2000; Li et al. 2015) and fiber density (Nie et al. 2012; Chen et al. 2013). We made initial attempts to examine their relationships with CNN predictions, as already illustrated in Figure 1. Cortical thickness information is available in HCP grayordinate data, whereas fiber density information is unavailable. Fiber density of a vertex is defined as the number of fibers penetrating its neighboring areas (Nie et al. 2012). Specifically, we extracted fibers from DTI data using fiber tractography method (Fillard et al. 2003), reconstructed DTI surfaces from white/gray matter segmentation of DTI images (Liu et al. 2007), and warped DTI surfaces to HCP grayordinate space using spherical registration method (Yeo et al. 2010). After that, we registered the fiber density information in DTI space to grayordinate space using a nearest neighbor method. Finally, we quantitatively measure the similarity between CNN predictions and the cortical thickness and fiber density patterns using mutual information (defined in Equation (2)). MI(X;Y)=∑x∈X∑y∈Yp(x,y)lnp(x,y)p(x)p(y) (2)where p(x,y) is the joint probability function of X and Y, and p(x) and p(y) are the marginal probability distribution functions of X and Y. In this paper, X denotes fiber density or cortical thickness, and Y represents CNN predictions. We carried out Monte Carlo simulations (1 000 simulations per subject) to access the statistical power of calculated mutual information.Wavelet Entropy MeasurementAmong wavelet-based entropy, Shannon entropy is the most typical representative, which quantifies the degree of signal complexity of signals (Sang et al. 2011). Formally, the Shannon entropy of a given vector t=[t1,t2,…,tn] is defined in Equation (3) (Coifman and Wickerhauser 1992). Generally, higher entropy reflects more random and complicated systems, and vice versa (Sang et al. 2011). In this work, we measured Shannon entropy of gyral/sulcal fMRI signals in the different frequency scales. We first performed 1D wavelet decompositions of original fMRI time series using discrete wavelet transformation based on Daubechies wavelet (Mallat 1989; Daubechies 1992), and then computed the Shannon entropy for the low-pass and high-pass components, respectively. Finally, we inspected whether there are significant differences between entropy distributions on gyri and sulci. Furthermore, to independently verified the CNN-based models, we employed wavelet entropy to discern cortical fMRI signals between gyri and sulci using random forest classification (RFC; Breiman 2001). Considerations of Physiological Noise and Signal-to-Noise Ratio DifferencesBlood-oxygen-level dependent (BOLD) hemodynamic response (HRD) is a low-pass filter in canonical HRD model (Glover 1999). High-frequency (>0.2 Hz) neural oscillatory activity may not be reflected at the oscillation frequency in fMRI data. Therefore, we call into question whether our observations are attributable to physiological noise or possible time course signal-to-noise ratio (tSNR) differences. On the one hand, the physiological noise, i.e., respiration and especially cardiac fluctuations related signals that exist in the fMRI time-series of the HCP data, will contribute to the high-frequency regime of the fMRI signals (Dagli et al. 1999). It is important to remove these noises before ascribing neural phenomena of sulco-gyral differences to high-frequency signal. To address this issue, we compared the group-level results with those obtained after appropriate physiological noise reduction with retrospective correction (RETROICOR, Glover et al. (2000)), on a subset of data from motor and working memory task. In addition, the higher frequency signal in sulci could also arise from cardiac fluctuations related fMRI signal in large sulci like Sylvian fissure (Dagli et al. 1999). We compared the results from the CNN models which were trained on the fMRI database without the signals in Sylvian fissure. On the other hand, different cortical areas might exhibit time course tSNR differences (Krüger and Glover 2001). Therefore, we measured the tSNR of grayordinates’ fMRI signals to inspect possible tSNR differences between gyri and sulci using the same procedure in our previous paper (Jiang et al. 2015). Specifically, we randomly chose 50 subjects from HCP S900 release and calculated the tSNR values of fMRI time series as Equation (4) (Krüger and Glover 2001; Murphy et al. 2007). We examined if there is statistical tSNR difference between gyri and sulci via unpaired 2-sample t-test (P < 0.01). tSNR=mean(timeseries)std(timeseries) (4)ResultsCNN Classification PerformanceWe summarize the classification performance of the proposed CNN model in Table 1. On average, the classification accuracies range from 74.0% to 76.2% in individual training fashion, and from 62.5% to 64.8% in group training fashion. Notably, the recall of sulcal examples are larger than that of gyral ones at both individual and group levels. The accuracy gap between individual level and group level might ascribe to fMRI noise confounds (Murphy et al. 2013) and inter-subject variability of hemodynamic responses (Dubois and Adolphs 2016). In addition, the CNN model achieved 73.6% accuracy on the training session and 67.2% accuracy on the testing session at individual subject level using HCP resting-state fMRI data. Generally, the accuracies in both training conditions are statistically beyond the chance level of 50%, indicating that there exist true discriminable temporal patterns between gyri and sulci. Table 1 The classification performance of the proposed CNN model at individual and group levels.
Spatial Distributions of CNN PredictionsFigure 3 displays the cortical representation of CNN predictions from individual training fashion (Fig. 3(a)) and group training fashion (Fig. 3(b)). Visual inspection of spatial distributions of CNN predictions suggests that there are only a small percentage of misclassified vertices over entire cortical surface. The majority of vertices are accurately classified, especially those on cingulate sulcus, the central sulcus, precentral sulcus, superior temporal sulcus, superior occipital gyrus, superior frontal gyrus, and superior temporal gyrus. In general, the reasonably good classification performance in both training fashions suggest the different functional roles that gyri and sulci may play. The lateral and medial cortical representation of CNN predictions from individual and group learning in working memory task. Red/blue indicates the vertices being predicted as gyri and sulci, respectively. Gray regions were neglected in the training procedure. (a) The spatial pattern of CNN predictions from one randomly chosen subject in individual training fashion. (b) The spatial pattern of average CNN predictions from all subjects in group training fashion. Frequency-Specific Differences Between Gyri and SulciIn this section, we present some qualitative examples and quantitative measurements, offering an intuitive understanding of the frequency-specific differences between gyri and sulci. Intuitive Visualizations of Gyral and Sulcal FiltersWe present some examples of gyral and sulcal components, including softmax weights (Fig. 4(a)) and time-series representations of convolution filters (Fig. 4(b)) and 2 representative gyral and sulcal components of GAP activations (Fig. 5). In Figure 4(a), filter #11 belongs to sulcal type, since its sulcal weight is positive (0.56) while its gyral weight is negative (−0.64). On the contrary, filter #24 belongs to gyral type, since its gyral weight is positive (1.01) and its sulcal weight is negative (−0.79). Quantitatively, there are 34 sulcal filters and only 15 gyral filters. In terms of time-series shapes, typical sulcal filter contains more high-frequency oscillations, while there are relatively less oscillations for gyral filters. The gyral filter (Fig. 5(a)) activates a majority of gyral areas in GAP layer, especially those on the dorsal occipital gyrus and bilateral supramarginal gyrus. The power spectrum of this gyral filter mainly concentrates in the low-frequency band, albeit with a high-frequency peak around 0.6 Hz. In contrast, the sulcal filter (Fig. 5(b)) activates a majority of sulcal areas, especially those on the central sulcus and precentral sulcus. The spectrum of this sulcal filter are heavy-tailed in the high-frequency band. Intuitive examples of gyral and sulcal components from one randomly selected subject in the working memory task derived from individual training fashion. (a) CNN weights in the softmax layer. (b) Time-series representations of corresponding CNN filters. The weights with black points are connecting to gyral units, while the weights with white triangles are connecting to sulcal units. Top panel (R1–R5 (#1–#2)), middle panel (R5 (#3–#8)–R7 (#1)) and bottom panel (R7 (#2–#8)–R8) represent sulcal, uncertain, and gyral filters, respectively. Two representative components from one randomly selected CNN model of in the working memory task. (a) Gyral components. (b) Sulcal components. Each panel includes GAP activations on cortical surface (left), time-series representation of corresponding filter (top right) and power spectrum of the filter (bottom right). The cortical areas in red and blue indicate activated and deactivated states, respectively. Quantitative Differences Between Gyral and Sulcal FiltersIn this section, we will demonstrate that above intuitive sulco-gyral differences are statistically significant. The bar chart in Figure 6 presents the numbers of sulcal and gyral filters in the learned CNN models from 7 HCP tasks. No matter at individual (Fig. 6(a)) or group (Fig. 6(b)) level, the sulcal filters are consistently more than gyral filters. Importantly, this finding is consistent across different tasks. Compared with individual-level results, the number of gyral and sulcal filters in group-level results changes across different tasks. Notably, the less similar number of filters in group-level results could result from inter-subject variances. This phenomenon remains to be explored in future. Tested by right-tailed paired-samples t-test, the quantity of sulcal filters is significantly large than the quantity of gyral filters (P < 10−6). Since these filters are inferred from original fMRI signals in a data-driven fashion, the statistically significant evidences implicitly reflect the intrinsic compositions differences between gyral and sulcal fMRI signals to some degree. In other words, sulcal fMRI signals might have more diversified temporal patterns than gyral fMRI signals. Quantitative differences between the numbers of sulcal and gyral filters in 7 HCP tasks. (a) Experimental results of individual training fashion. (b) Experimental results of group training fashion. The error bar indicates one standard deviation. Time-Frequency Analysis of Convolution FiltersWe present average frequency spectrums of gyral and sulcal filters in Figure 7. In general, the power spectrums of gyral filters concentrate in low-frequency band, and decrease quickly against increasing frequency. By contrast, sulcal filters exhibit broadband power spectrum, which is more heavy-tailed in high-frequency band and decrease slowly against increasing frequency. Notably, in the group-level classification, the “1/f” characteristic of the gyral power spectral density (PSD) plot is more pronounced while the sulcal PSD seems to exhibit a bimodal behavior with 2 high power PSD regimes. We performed one-tailed t-test on the total power spectrums of gyral and sulcal filters in low-frequency (0.01 Hz–0.08 Hz) and high-frequency (0.5 Hz–0.65 Hz) band to check whether the magnitude of gyral and sulcal filters are significantly different. The statistical results of individual training and group-level training are similar. For simplicity, we only summarize the group-level results in Table 2. In the low-frequency band, the magnitude of gyral filters are significantly larger than that of sulcal filters. In high-frequency band, the opposite is the case. The average frequency spectrums of convolutional filters in resting-state fMRI data are consistent with the spectrums in task-based fMRI data (Supplementary Fig. S3). Importantly, the time-frequency characteristics of gyral/sulcal filters demonstrated promising consistency across 7 task-based and resting sate fMRI data. The average power spectrum of gyral and sulcal filters over the HCP 7 tasks. (a) Results from individual training fashion. (b) Results from group training fashion. Black points and black triangles represent the spectrums of gyral and sulcal filters, respectively. Table 2 One-tailed t-test on the power spectrums of gyral and sulcal filters derived form group-level training fashion in 7 tasks (E: emotion; G: gambling; L: language; M: motor; R: relational; S: social; W: working memory).
Possible Underpinnings of Sulco-Gyral DifferencesFigure 8 displays the spatial distributions of CNN predictions, fiber density and cortical thickness from 3 randomly selected subjects in motor task (Fig. 8(a)) and the mutual information between them (Fig. 8(b)). It is clear to see that these spatial patterns well correspond to each other. The average mutual information for fiber density and cortical thickness is larger than 0.56 and 0.22 across all HCP tasks, respectively. Generally, fiber density have larger correlation with CNN predictionos than cortical thickness. Although the concordance between cortical thickness and CNN predictions is relatively lower, they are statistically significant (P < 10−3), tested by Monte Carlo simulations. Spatial patterns of fiber density and CNN predictions from 3 randomly chosen subjects in motor task and their mutual information. (a) Spatial patterns of CNN predictions, fiber density and cortical thickness. For CNN predictions, vertices in red and blue represent that they predicted as gyri and sulci, respectively. The vertices in green are uncertain samples and excluded before CNN classification. (b) Mutual information (1) between CNN predictions and cortical thickness (in orange) and (2) between CNN predictions and fiber density (in blue) (E: emotion; G: gambling; L: language; M: motor; R: relational; S: social; W: working memory). The dashed line indicates the random level. Wavelet Entropies of Cortical fMRI SignalsIn this section, we present the spatial patterns and detailed distribution histograms of wavelet entropy in Figure 9. The t-test results on the wavelet entropy of gyral and sulcal fMRI signals in low- and high-frequency bands are listed in Table 3. It is evident that gyral signals exhibit higher wavelet entropy values than sulcal in low-frequency band and vice versa, especially for those areas in the occipital lobe, frontal lobe and temporal lobe. Notably, this finding is also consistent in other HCP tasks (Supplementary Fig. S4). The above findings indicate that neural signals on gyri are more complex in low-frequency band while neural signals on sulci are more complex in the high-frequency band. Based on the wavelet entropy values in different frequency scales, the RFC classifiers achieved an average accuracy of 66.2% for the HCP 7 tasks. Although the accuracy of RFC is slightly worse than that of CNN, it is statistically beyond the chance level of 50%. This classification approach independently verified our CNN-based results. Cortical representation of wavelet entropies for low-frequency components (a) and high-frequency components (b) in HCP emotion task. Corresponding histograms of wavelet entropies for low-frequency components (c) and high-frequency components (d). Table 3 One-tailed t-test on the wavelet entropy of gyral and sulcal fMRI signals in low- and high-frequency band.
Physiological Noise and Signal-to-Noise Ratio DifferencesWe present the experimental results from CNN models with physiological noise reduction and without lateral sulcus in Supplementary Figs 5 and 6, respectively. These results with physiological noise reduction or without lateral sulcus are generally consistent with previous results in Figures 6 and 7, suggesting the identified sulco-gyral differences do not ascribe to physiological noise. The cortical tSNR differences are shown in Figure 10. We can see that average tSNR value of gyri is close to that of sulci, though gyri have larger standard deviation than sulci. Using unpaired 2-sample t-test, the results showed that the average tSNR values of gyri and sulci are statistically equal in all 7 tasks, which is similar to the findings using HCP Q1 data sets in our previous paper (Jiang et al. 2015). In conclusion, the identified sulco-gyral differences do not result from physiological noise or tSNR differences. The tSNR differences between gyri and sulci. (a) Cortical representation of tSNR values from one subject in emotion task. (b) Average tSNR value across 50 randomly chosen subjects in all 7 tasks (E: emotion; G: gambling; L: language; M: motor; R: relational; S: social; W: working memory). Discussion and ConclusionImplications of Functional Architecture of Cerebral CortexIn this paper, we performed a series of experiments to explore the frequency-specific differences between gyri and sulci based on HCP task-based fMRI data. Task-based fMRI data capture functional brain activities under specific task paradigms and can be used to identify brain regions and networks that are functionally involved in task conditions (Friston et al. 1994; Friston 2009). Task-based fMRI data is also a great source for us to discover and study the differences between gyri/sulci, however, powerful learning algorithms are urgently needed to better understand the differences. In this work, we adopted powerful learning model, such as CNN and wavelet entropy, to achieve more adequate understanding of the sulco-gyral differences. To the best of our knowledge, our work is one of the earliest studies to explore the potential frequency-specific pattern differences of cortical neural activities. The experimental results have shown that the proposed CNN approach is effective for classifying gyral/sulcal signals and inferring discriminative class-specific features. In addition, we find that sulcal fMRI signals are more diverse and more high-frequency than gyral signals. The wavelet entropy analysis further validated above findings. Notably, the detected sulco-gyral differences is not due to physiological noise or SNR differences in HCP fMRI data, and might truly reveal novel functional architecture of cortical gyri and sulci. In general, our obersvations are coherent with the inference of previous findings in different domains (Buzsáki et al. 2013; Deng et al. 2014). On the one hand, neuroscience of brain rhythms suggest that high-frequency brain activity reflects local domains of cortical processing, while low-frequency brain activity synchronize across distributed brain regions (Buzsáki et al. 2013). On the other hand, gyri are structural and functional hubs for global information exchange while sulci are functional units for local information processing (Nie et al. 2012; Deng et al. 2014). Combining with our findings, we could make inference that gyri and sulci play different functional roles in the working mechanisms of human cerebral cortex. Specifically, gyri with lower rhythms might play major roles in reflecting long distance communication and performing global functional integration, while sulci with higher rhythms tend to be in charge of local information processing and are critical for segregation of brain areas. In conclusion, these results revealed novel functional architecture of cortical gyri and sulci, supporting a new hypothesis that the cerebral cortex is bisectionally segregated into 2 fundamentally different units of gyri and sulci. Possible Determinants of the Frequency-Specific Differences Between Gyri and SulciWe correlated the functional role differences of gyri and sulci with the known attributes (e.g., fiber density and cortical thickness). Given the close relationship between the sulco-gyral differences and these known attributes, the distinctive functional activity patterns on gyri and sulci might be deeply rooted in the substrates of axonal fiber wiring patterns and cellular mechanisms. First, the cell type ratios and cortical thickness of cerebral cortex in different regions are different. It was found in Hilgetag and Barbas (2005) that gyri have more neuron numbers than sulci, especially in the deep layers (V + VI). Because different layers consist of different types of cells (for example, it was found in Lam and Sherman (2010) that abundant pyramidal neuron axons in Layer V run to subcortical structures (such as the basal ganglia), while small spindle-like pyramidal and multiform neurons in Layer VI sends efferent axons to the thalamus), the differences between functions and rhythms of gyri and sulci could be affected by the components of different cell types. Also, the cell morphology could also affect the function of cell. In Hilgetag and Barbas (2005), dendrites of pyramidal cells in gyri were found to be elongated, resulting a greater attenuation of excitatory postsynaptic potentials. Therefore, gyro-sulcal functional activity differences fundamentally result from the difference at cellular level, while how to explicitly relate the cell signaling to more global signal coupling and oscillation patterns is still to be better elucidated. Second, the subcortical axonal wiring patterns are different between gyri and sulci. Axons appear to be connected to gyri with greater density (Deng et al. 2014). This hypothesis is proposed in our previous work based on DTI dataset (Nie et al. 2012), and gains a growing number of supports from a variety of field, such as dissection studies (Xu et al. 2010) and histological studies (Budde and Annese 2013). This hypothesis could induce a speculation that gyri have stronger structural and functional connections to other cortices and subcortical regions. These observations and induced speculations implicitly suggest that gyri could be hubs to segregate and integrate information from multiple regions in contrast to sulci. Given the larger correlations between CNN predictions and fiber density, among the 2 potential determinants, deployment of long-range axonal fibers could be a key factor in the preserved time management in brains (Buzsáki et al. 2013). The underlying determinants is still to be better elucidated. Identifying signaling pathways of the described frequency fluctuations of gyri and sulci deserves further study to gain explicit insight into the cell signaling to more global signal coupling and oscillation patterns. Considerations and LimitationsThe works in this paper have some major limitations and can be enhanced in following directions. First, we mainly explored the frequency differences at macro gyri/sulci level, while the cytoarchitectonics of human brains are much more sophisticated (Amunts et al. 2007). For example, primary, secondary and tertiary gyri develop at different pre and postnatal periods and might have different mechanisms of morphogenesis (Sun and Hevner 2014; Budday et al. 2015). Whether our findings hold for these fine-grained brain structures remain unanswered. Second, BOLD mechanism tracks hemodynamic response in an indirect and sluggish manner. Hence, high-frequency neural oscillatory activity may not be reflected at the oscillation frequency in fMRI data (Chen et al. 2017). Our results need to be further verified using multimodal imaging data, such as intracranial electroencephalography (iEEG, Engel et al. (2005)) and magnetoencephalogram. Last but not the least, while impressive, our designed CNN model in this work are relatively simple. We made an initial attempt to train deeper CNN models and achieved higher accuracy performance (as shown Supplementary Figure S7). As we can see that as the layer goes deeper, the classification performance on training data grows from 66.1% to 69.7% while accuracy on testing data increase from 65.0% to 66.1%. In general, the deeper the convolutional networks, the higher the classification performance. An interesting extension to the present work would be applying more powerful representation learning models, such as much deeper CNN model (He et al. 2016), recurrent neural network (Mikolov et al. 2010), and generative adversarial networks (Goodfellow et al. 2014), which could reveal hierarchical feature differences between gyri and sulci. FundingX. Jiang was supported by the China National Science Foundation (NSFC) 61 703 073, and the China Special Fund for Basic Scientific Research of Central Colleges ZYGX2017KYQD165. T. Liu was supported by the National Institutes of Health (DA-033393, AG-042599) and by the National Science Foundation (NSF CAREER Award IIS-1149260, CBET-1302089, BCS-1439051, and DBI-1564736). NotesWe gratefully acknowledge anonymous reviewers for their constructive comments, which have helped the authors to greatly improve the quality of this article. Conflict of Interest: None declared. References, , , , , , , , , , et al. . Function in the human connectome: task-fMRI and individual differences in behavior . . :–., , , , . . A developmental and genetic classification for malformations of cortical development: update 2012 . . :–., , , , . . Igf1 gene disruption results in reduced brain size, CNS hypomyelination, and loss of hippocampal granule and striatal parvalbumin-containing neurons . . :–., , , , , , , , , . . Coevolution of gyral folding and structural connection patterns in primate brains . . :–.. Ten lectures on wavelets. Society for Industrial and Applied Mathematics. , , , , , . . Single or multiple frequency generators in on-going brain activity: a mechanistic whole-brain model of empirical MEG data . . :–., , . Quantitative Analysis of White Matter Fiber Properties along Geodesic Paths. Medical Image Computing and Computer-Assisted Intervention. 16–23. , , , , , . . Statistical parametric maps in functional imaging: a general linear approach . . :–., , , , , , , , , , et al. . The minimal preprocessing pipelines for the Human Connectome Project . . :–., , . Convolutional Neural Networks for Acoustic Modeling of Raw Time Signal in LVCSR. 16th Annual Conference of the International Speech Communication Association. 26–30. , , . Generative adversarial nets. International Conference on Neural Information Processing Systems. 2672–2680. , , . Deep Residual Learning for Image Recognition. 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). 770–778. , , , , , , . . Sparse representation of HCP grayordinate data reveals novel functional architecture of cerebral cortex . . :–., , , , , , , , , . . Temporal dynamics assessment of spatial overlap pattern of functional brain networks reveals novel functional architecture of cerebral cortex . . :–., , , , , . . Spatiotemporal patterns of cortical fiber density in developing infants, and their relationship with cortical thickness . . :–., , . Network in network. International Conference on Learning Representations. 1–8. , , , , , , , . . Elucidating functional differences between cortical gyri and sulci via sparse representation HCP grayordinate fMRI data . . :–., , , , , , , . . Quantitative in vivo measurement of gyrification in the human brain: changes associated with aging . . :–., , . Recurrent neural network based language model. Conference of the International Speech Communication Association. 1045–1048. , . Rectified linear units improve restricted boltzmann machines. International Conference on International Conference on Machine Learning. 807–814. , , , , , , , , , , et al. . Axonal fiber terminations concentrate on gyri . . :–., , , , , , , , , . . Radial structure scaffolds convolution patterns of developing cerebral cortex . . :–., , , , , , , , . . MSM: a new flexible framework for multimodal surface matching . . :–., , , , . . Wavelet-based analysis on the complexity of hydrologic series data under multi-temporal scales . . :–.. . Unique morphological features of the proliferative zones and postmitotic compartments of the neural epithelium giving rise to striate and extrastriate cortex in the monkey . . :–., , , , , , , , , , et al. . Trnp1 regulates expansion and folding of the mammalian cerebral cortex by control of radial glial fate . . :–., , , , , , , , , , et al. . Pushing spatial and temporal resolution for functional and diffusion MRI in the human connectome project . . :–., . . Different frequencies for different scales of cortical integration: from local gamma to long range alpha/theta synchronization . . :–., , . Time series classification from scratch with deep neural networks: A strong baseline. International Joint Conference on Neural Networks. 1578–1585. , , . Learning Deep Features for Discriminative Localization. 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). 2921–2929. , , , , , , , , , , et al. . Quantitative analysis of sulci in the human cerebral cortex: development, regional heterogeneity, gender difference, asymmetry, intersubject variability and cortical architecture . . :–. |