Functional network estimation using multigraph learning with application to brain maturation study.

Although most dramatic structural changes occur in the perinatal period, a growing body of evidences demonstrates that adolescence and early adulthood are also important for substantial neurodevelopment. We were thus motivated to explore brain development during puberty by evaluating functional connectivity network (FCN) differences between childhood and young adulthood using multi-paradigm task-based functional magnetic resonance imaging (fMRI) measurements. Different from conventional multigraph based FCN construction methods where the graph network was built independently for each modality/paradigm, we proposed a multigraph learning model in this work. It promises a better fitting to FCN construction by jointly estimating brain network from multi-paradigm fMRI time series, which may share common graph structures. To investigate the hub regions of the brain, we further conducted graph Fourier transform (GFT) to divide the fMRI BOLD time series of a node within the brain network into a range of frequencies. Then we identified the hub regions characterizing brain maturity through eigen-analysis of the low frequency components, which were believed to represent the organized structures shared by a large population. The proposed method was evaluated using both synthetic and real data, which demonstrated its effectiveness in extracting informative brain connectivity patterns. We detected 14 hub regions from the child group and 12 hub regions from the young adult group. We show the significance of these findings with a discussion of their functions and activation patterns as a function of age. In summary, our proposed method can extract brain connectivity network more accurately by considering the latent common structures between different fMRI paradigms, which are significant for both understanding brain development and recognizing population groups of different ages.

INTRODUCTION

Modern neuroimaging techniques provide an opportunity to study the human brain quantitatively. Especially, functional magnetic resonance imaging (fMRI) enables noninvasive measurements of both brain structure and its functional activities with a temporal resolution of seconds and with a spatial resolution of millimeters. It measures the low‐frequency fluctuations caused by hemodynamic response in blood‐oxygen‐level dependent (BOLD) signals and maps to the neural activities in the brain or spinal cord (Huettel, Song, & McCarthy, 2004), which facilitates the detection of correlations among distinct brain regions (Power et al., 2011). In particular, task‐based fMRI (tfMRI) analyses help to identify and characterize functionally distinct nodes in the human brain, providing the ability to interrogate the neural basis of mental functions and representations (Barch et al., 2013).

The human brain has attracted a significant number of studies due to its broad functional repertoire including action enabling, perception, and cognition. Despite a fixed anatomical structure, it has been a widely accepted assumption that the brain network architecture organizes local interactions to cope with diverse environmental demands (Park & Friston, 2013). Considerable research has been conducted to reveal the development of the brain from their functional activation differences throughout childhood and adolescence (Fair et al., 2009; Jolles, van Buchem, Crone, & Rombouts, 2010). In particular, task fMRI has been widely studied to explore task specific functional behavioral difference between child and adult brains (Passarotti et al., 2003), where different activation patterns were identified and reported for different age groups (Holland et al., 2001).

Recently, the graph signal processing (GSP) methodologies have become promising to tackle neuroimaging problems where observations are viewed as signals residing on the nodes of a graph (Hu et al., 2013; Sandryhaila & Moura, 2013; Shuman, Narang, Frossard, Ortega, & Vandergheynst, 2013). Within this framework, the brain is modeled as a graph or network consisting of a set of brain regions and pairwise relationships between these regions including connectivity or similarity measures. Spectral analysis was implemented to reveal the intrinsic structures of these graphs via the eigen‐analysis of their associated adjacency matrix, the graph Laplacian and their variants (Hu et al., 2013; Sandryhaila & Moura, 2013). A remarkable number of researches have been conducted to investigate the frequency domain behavior of the graph signals residing on these matrices (Chung, 1997), which facilitates the development of GSP (Shuman et al., 2013). Among them, the graph Fourier transform (GFT) has demonstrated its effectiveness in signal denoising (Shuman et al., 2013), the study of mental illness (Hu et al., 2013), and brain network analysis (Huang et al., 2016). Variants of GFT such as alternative graph Fourier transform (AGFT) (Sandryhaila & Moura, 2013) have been successfully applied to signal compression and label propagation. The rationality and advantages of these GSP tools over conventional Fourier transform have been illustrated by Shuman et al. (Shuman et al., 2013). These methods rely on a preselected kernel to construct the graph, presuming that the prior inner relationships among data samples are known. However, such information is unavailable because of the complexity of the brain with distinct task‐specific regions activated and deactivated among individuals. To overcome this issue, a data‐driven graph learning framework (Dong, Thanou, Frossard, & Vandergheynst, 2016) was proposed to better extract the graph structures. The authors further applied graph learning with graph Fourier transform and found applications in studying brain development (Wang et al., 2020).

It has been demonstrated that integrating multi‐modality or multi‐omics data can yield additional informative knowledge for the improvement of predictive accuracy (Kim et al., 2014; Sundermann, Herr, Schwindt, & Pfleiderer, 2014; Tsuda, Shin, & Schölkopf, 2005). In this work, we extend our model in (Wang et al., 2020) to incorporate multiple graphs, which can integrate multiple pieces of information to make diagnoses for patients. Specifically, we propose a multigraph learning framework to jointly learn the graph structure from multiple graphs. Afterward, we conduct GFT and design filters to detect essential brain regions in a specific frequency range to study the mechanisms of brain maturation. Finally, we apply the proposed model to neuroimaging data collected from the Philadelphia Neurodevelopmental Cohort (PNC) (Satterthwaite et al., 2014). It is a large‐scale collaborative study between the Brain Behavior Laboratory at the University of Pennsylvania and the Center for Applied Genomics at the Children's Hospital of Philadelphia, including healthy developing volunteers of age 8 to 22 years. The database contains 652 subjects with task fMRI observations available for both emotion identification and working memory tasks.

The contributions of this work can be summarized as follows.

• We propose a multigraph Laplacian learning framework to jointly estimate the graph structures from neuroimaging studies involving multiple paradigms. This is different from the previous graph integration methods (Kim et al., 2014; Tsuda et al., 2005), graph concatenated methods (Hu et al., 2013), and other multi task learning models (Wee, Yap, Zhang, Wang, & Shen, 2014), which combined the graphs estimated from single modality data linearly, calculated the graph from concatenated data, or estimated the identical graph among all the subjects, respectively. Instead, our model estimates multiple graphs with respect to different paradigms simultaneously by considering the overall brain network structures. The complementary information from multiple paradigms enables a more reliable construction of the brain network, resulting in improved detection of hub regions.

• We evaluate the proposed learning framework first on synthetic data, and then conduct extensive experiments on the real PNC data. The results demonstrate the superiority of the proposed method over traditional GSP approaches. A list of hub regions were identified revealing the different activation patterns in child and young adult group. A detailed discussion of the detected regions is given regarding their functions and relationships with brain maturation.

• We summarize our findings to demonstrate the consistency and discrepancy with previous studies and explore the potential biological patterns of the brain maturation process. Therefore, our work benefits the neuroimaging field by providing additional biomarkers related to brain development during adolescence.

The rest of the paper is organized as follows: in the Methods section, we provide a detailed description of signals defined on the graph and review the developments of applying GFT on these graph signals. We further propose a multigraph learning model to estimate the graph Laplacian matrix in the multi‐paradigm setting, followed by graph Fourier analysis and filter design. In the Results section, we validate the effectiveness of the proposed model on both synthetic and real data. Finally, we summarize the work in the Conclusion section.

METHODS

Signals defined on a graph

We use a graph to describe the pairwise relationships between a set of objects. For example, given a set of nodes , a corresponding graph can be denoted as , where refers to a adjacency matrix of the relationships between each node‐pair. Brain functional connectivity aims to study the connection and activation pattern between different brain regions. As a result, we treat each brain region as a node and consequently the adjacency matrix denotes the connections between different brain regions. A graph signal is defined on the vertex set of the graph , where the ‐th element of represents the strength of the neural activity of the corresponding brain region.

We further introduce the graph Laplacian defined below:where is the degree matrix, which is diagonal with the ‐th diagonal entry being the sum of the entries in the ‐th row of the weighted adjacency matrix . The graph Laplacian in (1) is real and symmetric, and thereby has a complete set of orthonormal eigenvectors {. By eigendecomposition, we have , where Λ is the diagonal eigenvalue matrix Λ, and is the corresponding eigenvector matrix.

To further illustrate the behavior of the graph signal defined on the network, we define the smoothness (Laplacian regularization) of a graph signal with respect to the Laplacian as (Belkin & Niyogi, 2003; Chan, Osher, & Shen, 2001; Shuman et al., 2013)where represents the element of at ‐th node. It serves as an indicator reflecting the smoothness properties of a signal on the graph.

Related works

In this subsection, we review the recent developments in GSP especially graph Fourier transforms.

GFT and AGFT

Suppose we have subjects in total. For each subject, we have graph signals or [ available, where and denote the number of nodes and graph signals. is the ‐th column of denoting the ‐th graph signal residing on the vertices, and is the ‐th row of denoting the observations of the ‐th vertex across the graph signals. Let denote the association or similarity matrix for a subject (Hu et al., 2013; Shuman et al., 2013).where indicates an edge between node and node , and is the hyper‐parameter that determines the width of Gaussian distribution. The GFT of with respect to is defined as (Hu et al., 2013; Huang et al., 2016; Shuman et al., 2013)The inverse GFT of with respect to is defined asNote that is orthonormal, that is, , thus and form a GFT pair.

The eigenvalues of the graph Laplacian carry the similar meaning as “frequency” to that of conventional Fourier transform. In fact, the eigenvalue reflects how much the corresponding eigenvector oscillates over the graph if we calculate the Laplacian regularization of the eigenvector :In (Sandryhaila & Moura, 2013, 2014), an alternative GFT was proposed to expand graph signals with respect to the eigenfunctions of the corresponding adjacency matrix. Suppose the adjacency matrix is available and its Jordan decomposition , where is the Jordan normal form and is the matrix of generalized eigenvectors. The alternative GFT and its inverse of a given graph signal can be defined with respect to as in (4), (5).

Graph Laplacian learning based Fourier transform (GLFT)

In order to avoid the deviation of pre‐selected kernels in constructing the network, the authors in (Dong et al., 2016) proposed a data‐driven learning framework to better estimate the topology of the graph.

Suppose we have graph signals observed on a graph with nodes. The smooth signal and the graph Laplacian were estimated by minimizing the Laplacian regularization term as well the loss (Dong et al., 2016), where Frobenius norm penalty of and constrains were applied to ensure the learned graph Laplacian is nontrival and valid. Detailed description can be found in File S1.

Notably, the proposed learning algorithm was further combined with GFT in (Wang et al., 2020) to split the graph signals into different frequency components. Given a graph signal and its learned graph Laplacian, the GLFT and the inverse of can be defined with respect to as in (4), (5).

Multigraph analysis

Multigraph Laplacian learning

Graph signals from different views or paradigms may provide complementary information for better characterizing the network structure. In this section, we extend the learning framework into multiple graphs setting. To begin, we introduce the notions of the multigraph signals. For a subject with graph signals defined on modalities, we use to represent the set of graph signals from the ‐th( modality, where is the number of graph signals from the ‐th modality. For each modality, by taking the smoothness into consideration, the graph Laplacian can be easily estimated (File S2).

In multigraph scenario, we still adhere to the assumption that the graph signals have a small total variation on the estimated graph for each modality. This assumption reflects a better fit between the learned graph structure and the graph residing on it. Moreover, since the multiple signals come from the same subject, common structures reflecting interrelationships between modalities should exist. Based on these assumptions, it is natural to generalize the graph Laplacian learning framework into a multigraph setting using the following objective function.Compared with the single variate learning framework (File S2), we considered the overall graph structures by regulating their difference among modalities. The positive parameters are the key factors reflecting the latent structure shared by different views or modalities. If is set to 0, which means that the graph structures are considered independently with each other, the solutions are equivalent to learning the graph structure from each modality separately as in (File S2). If is set to infinity, which forces the graph structures to be the same among modalities, the solutions are equivalent to learning the graph with the concatenated method. Ideally, should be set as the degree of similarity measurement between modalities with certain prior information. However, in our case, the inner‐relation between modalities is rather sophisticated, so is tuned through cross‐validation in this work. is a positive regulator that controls the distribution of the learned graphs and we tune it through cross‐validation as well. Additional constraints are enforced to prevent the null solutions and to ensure that the learned Laplacian matrices are valid. The optimization problem can be effectively solved with CVX toolbox (Grant, Boyd, & Ye, 2015).

Graph Fourier transform

After obtaining graph Laplacian for each modality, we are able to conduct GFT of in each modality with respect to aswhere is the eigenvector matrix of . In addition, the inverse GFT of with respect to is defined asAs is orthonormal, the original signal can be recovered through the inverse GFT of . Thus, and form a GFT pair.

Remarkably, the GFT can be reduced to conventional Fourier transform for specific graphs such as a cycle graph (Huang et al., 2016), where the eigenvectors in are expressed as . On the other hand, as defined in (6), the eigenvalue reflects the smoothness of the corresponding eigenvector over the graph. Thus, the GFT naturally divides the original graph signal into different “frequency” components: small eigenvalues oscillate slowly, appearing to be smooth, whereas the large eigenvalues account for high frequencies, and the associated eigenfunctions oscillate rapidly.

Graph filtering

As discussed above, the transformed signals can be manipulated into a range of frequencies to facilitate extracting information and exploring their behaviors. Following the standard procedure from signal processing, we design a low‐pass filter, a high‐pass filter, and a band‐pass filter based on eigenvalues to split the frequency components (see supplementary section 3).

Based on the designed filters, the graph signals were naturally divided into low, intermediate, and high frequency components as . Accordingly, the original graph signals can be recovered by with being the reconstructed signals from the respective frequencies.

Brain network theory has demonstrated the co‐existence of disorganized behavior and ordered regularity in the brain functional network (Sporns, 2010). The organized part reveals the functional or structural information of the brain network in a large population, such as the default mode network (DMN) (Buckner, Andrews‐Hanna, & Schacter, 2008), sensory motor network (Gallese & Lakoff, 2005), while the disorganized part reflects the individual differences among subjects. In our multigraph learning framework, the ordered part corresponds to low frequency components that are smooth across the graph representing the main contribution of the brain functional operations, while the disordered portion accounts for the high frequency component with a higher chance of capturing noise or individual differences. The low variability has proven to be essential in the analysis of neurological disease and behavior and the learning process (Garrett, Kovacevic, McIntosh, & Grady, 2012; Huang et al., 2016). As the main motivation of this work is to analyze the brain functional connectivity patterns at different pubertal stages, we put our emphasis on the learning of the low frequency subspace in the following section.

RESULTS

We evaluated the performance of our proposed multigraph Fourier transform framework using both synthetic and real data analysis in this section. The comparison between other graph Fourier transform related methods such as GFT, alternative GFT, and the graph learning based Fourier transform was demonstrated (Sandryhaila & Moura, 2013; Shuman et al., 2013; Wang et al., 2020). Meanwhile, hub regions related to brain maturation were identified from task‐related fMRI, followed by a significant test of the detected regions. Finally, we calculated the consistency between subjects for each group and conducted a classification test to further validate the performance of our proposed method. In short, our proposed multigraph learning method enabled a more precise characterization of the brain FCN with multi‐paradigm fMRIs, and thus facilitates the study of brain development. For clarification, we illustrated the major steps involved in real data analysis in Figure 1.

FIGURE 1

A flow diagram of multigraph signal analysis of one subject for brain network construction and hub regions identification. There are four major steps involved. In the preprocessing step: task fMRI measurements from multi‐paradigm were preprocessed to become graph signals defined on the parcellated brain regions according to the Power atlas (Power et al., 2011). In the multigraph learning step: we applied the proposed multigraph Laplacian learning algorithm to construct the brain connectivity networks simultaneously for each paradigm. In the GFT‐filters step: we conducted graph Fourier transform and designed proper filters to transform the original signals into the frequency domain and put our emphasis on the low frequency components. Finally, in the mapping back step: we analyzed the energy distribution of the eigenfunctions and mapped these back to their corresponding brain regions, where regions with significant contribution were identified

Synthetic data analysis

Artificial graphs construction

To evaluate the performance of our proposed model, we conducted an experiment on synthetic multigraph data. Specifically, we synthesized artificial graphs with 20 vertices following the Erdos‐Renyi (ER) model (Erdos & Rényi, 1960) and the Barabasi‐Albert (BA) model (Barabási & Albert, 1999). The ER model first generates the whole node set, and then each edge randomly with a .2 probability. The BA model, by contrast, adds nodes and edges one by one, according to the degree distribution of the existing nodes, based on the assumption that there is a higher chance for a new vertex connecting to an existing node with higher degree. The ground truth graph Laplacian was obtained from (1), where the edges were binarized. For each artificial graph, a total of 100 signals were yielded following a zero‐mean multivariate Gaussian distribution as followsWe set and for two paradigms, respectively.

Result and comparison

Given the synthesized multi‐paradigm signals on each graph, we first applied the proposed learning framework to construct the graph Laplacian matrix. For comparison, we also calculated the graph Laplacian based on single‐paradigm data, which demonstrated superiority over GFT and AGFT (Wang et al., 2020). For ease of comparison, the learned graph Laplacian matrices were binarized at a threshold 10. Precision, recall, and F‐measure were adopted as the criteria to evaluate the performance of these methods statistically. We displayed the averaged results and their SD from 20 repeated experiments in Table 1. The results demonstrated that the constructed graph networks using the multigraph learning framework outperformed those using the single‐paradigm based method. In other words, the overall structural information from multiple paradigms offers advantages in better constructing the graph than in single paradigm. The comparatively poorer performance of the ER model was most likely caused by the higher randomness against smoothness condition.

TABLE 1

Simulation results

BA model
Graphs	F‐measure	Precision	Recall
2D‐L1	0.9870 ± 0.0215	0.9474 ± 0.0292	0.9921 ± 0.0193
2D‐L2	0.9936 ± 0.0114	0.9875 ± 0.0222	1.0000 ± 0.0000
L1	0.8723 ± 0.0954	0.8539 ± 0.0920	0.8921 ± 0.1017
L2	0.8545 ± 0.0815	0.8366 ± 0.0801	0.8737 ± 0.0861
ER model
Graphs	F‐measure	Precision	Recall
2D‐L1	0.7892 ± 0.0337	0.7848 ± 0.0335	0.7937 ± 0.0344
2D‐L2	0.8033 ± 0.0240	0.7984 ± 0.0248	0.8083 ± 0.0240
L1	0.7043 ± 0.0288	0.6993 ± 0.0300	0.7094 ± 0.0283
L2	0.7380 ± 0.0272	0.7324 ± 0.0273	0.7438 ± 0.0280

Note: 2D‐L1 and 2D‐L2 represent the graphs jointly estimated from two paradigm signals. L1 and L2 are the graphs recovered from graph learning algorithm for single paradigm.

For clarity, we also visualized the learned graph Laplacian matrices with respect to both BA and ER graphs in Figure 2. The experiments are repeated 20 times and the visualization is based on the results from the last time. Overall, our proposed method gives a better reconstruction of the ground truth graph compared with the single view based method.

FIGURE 2

The ground truth graphs versus the reconstructed graphs from our proposed multigraph learning framework and the graph learning method for single view data. Note that the edges in the last four sub figures represent the falsely predicted edges. The upper part is the simulation from BA model and the lower part is the ER model. 2D‐L1 and 2D‐L2 represent the Laplacian matrices estimated from the proposed multi‐paradigm learning framework. GL‐L1 and GL‐L2 represent the Laplacian matrices learned from single view graph signals

Real data analysis

We further applied the proposed model to fMRI observations collected from the Philadelphia Neurodevelopmental Cohort (PNC) (Satterthwaite et al., 2014) with a focus on the brain maturation process. MRI examinations were conducted on a single 3 T Siemens TIM Trio whole‐body scanner. The images were collected using a single‐shot, interleaved multi‐slice, gradient‐echo, echo planar imaging sequence.

Data acquisition and preprocessing

In this study, we consider two types of task MRI data from PNC (Satterthwaite et al., 2016) to reveal the brain network difference between child and young adult given the same cognitive task. The two tasks involve emotion identification and working memory examination, and a detailed description of the experiments can be found in the paper by Satterthwaite et al. (2014). For each subject, both task based fMRI time series were available and parcellated based on the same atlas.

The emotion identification task employs a fast event‐related design with a jittered inter‐stimulus interval (ISI). Participants were required to label the emotions displayed by the trained actors (50% female) with neutral, happy, sad, angry, or fearful expressions. There were 60 faces in total in color photographs that were rated and selected by professional directors. Each face was displayed for 5.5 s followed by a variable ISI from 0.5 to 18.5 s, during which a complex crosshair (that matched the faces' perceptual qualities) was displayed. The task duration was 10.5 min (Satterthwaite et al., 2014).

The working memory task involved a fractal version of the standard n‐back task, which has proved to be a reliable probe of the executive system and was beneficial to avoid lexical processing abilities that evolve during development (Brown et al., 2004; Schlaggar et al., 2002). The task involved presentation of complex geometric figures (fractals) for 500 ms, followed by a fixed ISI of 2,500 ms. There are three levels of working memory load in total: 0‐back, 1‐back, and 2‐back. In 0‐back setting, participants responded to a specified target fractal; in 1‐back setting, participants responded if the current fractal was identical to the previous one; in 2‐back setting, participants responded if the current fractal was identical to the previous two trails. The task duration was 11.6 min (Satterthwaite et al., 2014).

We followed a standard pipeline to preprocess the data using SPM12 (Ashburner et al., 2014), consisting of motion correction, spatial normalization to standard MNI space (adult template) and spatial smoothing with a 3 mm FWHM Gaussian kernel. We applied regression to remove the effect of motion (6 parameters), followed by a band‐pass filter in 0.01 Hz to 0.1 Hz frequency range. Next, we further mapped the imaging data into 264 regions of interest based on the Power atlas (Power et al., 2011) with a sphere radius 5 mm. Finally, we extracted the ROI‐level fMRI time series by averaging the time sequences of all voxels in the same ROI to get a 264 ×  matrix for each subject with and denoting the number of time points for emotion identification and working memory tasks, respectively, due to the total duration of the MRI scan.

Since we focus on the brain development analysis as a function of age, we subdivide the subjects into child and young adult groups. Specifically, subjects whose age was below 12 years were considered to belong to the child group, while subjects over 18 years belong to the young adult group, as shown in Table 2. Chi‐Square statistic test (McHugh, 2013) was performed and there was no significant difference in distribution of males and females between the child and young adult groups ().

TABLE 2

Characteristics of the subjects in this study

Group	Age (Mean ± SD)	Male/female
Child	10.2926 ± 0.9358	61/74
Young adult	19.2606 ± 0.9952	57/92

Parameter tuning/ effect of parameters

We illustrated the determination of the important parameters involved in the proposed framework, including that controls the similarity between paradigms, that regulates the off diagonal entrees of the graph Laplacian, and that determines the cutoff of the low frequency components. The parameters were determined separately and 10‐fold cross validation was applied to evaluate the performance of the parameters. Specifically, was initially set in the range of [10, was initially set in the range of [10, and cutoff was set in [20, 60]. The cross validation criteria was the classification accuracy of the estimated graph Laplacian towards their corresponding age groups.

Based on the objective function (7), a large (over 10) leads to a uniform distribution of the learned graph Laplacian (the diagonal close to one and the off‐diagonal tends to be the same) and a small (below 1, compared with the scale of the Laplacian regularization term controlling the smoothness) leads to an extremely sparse solution where a clear majority of off‐diagonal elements equal 0. The parameter reflects the inner structural relationships between paradigms by controlling the shared latent information. We visualized the relationship of the parameters () versus cross‐validation accuracy to evaluate the effect of the parameters in Figure 3. The cutoff splitting the frequency was also selected via cross‐validation, where 10‐fold cross validation was performed on the train set with in [20, 60]. We finally adopted , and as the low frequency boundary and parameter for the optimization algorithm. In fact, the classification performance was robust across a range of both parameters.

FIGURE 3

The visualization of the sets of parameters versus their corresponding classification accuracy in classifying different age groups. The is set to 10 and the x‐axis is the . For each pair of the parameters, the classification task was repeated 10 times and the averaged cross‐validated accuracy and standard derivation were displayed. We concluded that the classification performance is solid and robust (over 95% accuracy and fluctuation is within 3 × 10) with a proper range of the parameters. A similar result is observed when is at a proper range, which demonstrated the robustness of our model

Detection and visualization of hub regions

As illustrated in Figure 1, we applied our proposed framework to detect age related biomarkers using imaging data from two task‐based fMRI paradigms. Specifically, fMRI observations at each time point were treated as graph signals defined on the parcellated brain regions. Multigraph Laplacian learning algorithm was applied on both working memory and emotion identification tasks to estimate the graph structure for each subject. We conducted graph Fourier transform and designed filters afterwards to divide the graph signals into low, intermediate and high frequency components. Based on the explanations at the end of method section, we put our emphasis on the low frequency portion. In order to identify hub regions for each age group, we calculated the projection energy on eigenbasis in the low frequency range, where the spikes of the energy concentration were selected and further mapped back to their corresponding cortical regions. Specifically, the inner products of the transformed signals in the low frequency domain on each eigenbasis were computed and two standard deviations higher than mean value (Bassett et al., 2008) was adopted as the threshold to select hub regions. The energy distribution of the transformed signals provided additional evidence on the superiority of our proposed method (Figure 4). The relatively higher energy concentration over competing method in the low frequency ranges indicated that our proposed method could facilitate the extraction of the information in the organized brain regions, which in turn improved the subsequent hub identification and classification. We investigated the group differences of the detected hub regions using standard two‐sample t test to estimate their significance. We summarize the results in Table 3, where previous reported regions are marked with *. Detailed illustrations of these regions and their literature references are given in Discussion Section. For ease of visualization, we display the locations of the identified hub regions related to brain maturity in Figure 5 with the help of BrainNet Viewer (Xia et al., 2013).

FIGURE 4

The comparison of projection energy between our proposed method and the RBF (Radial basis function) kernel based method involved in the traditional graph Fourier transform (Shuman et al., 2013). The projection energy was calculated based on the averaged projection power across all the subjects and normalized to have unified area under the curve. The eigenbases are sorted in an increasing manner with corresponding eigenvalues except for the first one, which has constant value in its eigenfunction. Our proposed method enforced the extraction of the features in the low frequency components compared with the traditional methods

TABLE 3

Identified hub regions in different age groups

	Child		Young adult
Common regions	Location	p‐value	Location	p‐value
Left cingulate gyrus*	2.6606e−4	Left cingulate gyrus*	2.6606e−4
Right cingulate gyrus*	1.8612e−10	Right cingulate gyrus*	1.8612e−10
Right postcentral gyrus*	.0439	Right postcentral gyrus*	.0439
Right medial frontal gyrus *	2.3261e−6	Right medial frontal gyrus *	2.3261e−6
Left posterior cingulate*	2.1909e−4	Left posterior cingulate*	2.1909e−4
Right inferior parietal lobule*	.0019	Right inferior parietal lobule*	.0019
Left Claustrum	4.4006e−4	Left Claustrum	4.4006e−4
Left superior frontal gyrus*	1.3694e−5	Left superior frontal gyrus*	1.3694e−5
Right inferior frontal gyrus*	.0035	Right inferior frontal gyrus*	.0035
Left middle frontal gyrus*	7.7311e−4	Left middle frontal gyrus*	7.7311e−4
	Left middle temporal gyrus*	2.3678e−4	Left Postcentral gyrus	.0796
Different regions	Right cuneus*	.1352	Right paracentral lobule*	8.9501e−5
	Left insula	.0072
	Right inferior temporal gyrus*	4.73e−10

Note: Previously reported regions are marked with *.

FIGURE 5

The visualization (sagittal, axial, and coronal views) of the detected hub regions from task fMRI analysis. The regions in the upper part of the figure are identified from the child group while the regions in the lower part of the figure are detected from the young adult group. The visualization of the brain network is by BrainNet Viewer (Xia, Wang, & He, 2013)

Subject consistency

To validate our previous statement that organized brain functional and structural information should be consistent within a large population, the subject consistency can be adopted as an indicator to compare the subjects' behavior in the same group (Li et al., 2019). It is well known that Dice similarity coefficient (DSC) is a commonly used metric to quantify the subject consistency:where and denote the true positive, false positive, and false negative, respectively. A high DSC score indicates a high consistency of the functional connectivity networks among subjects from the same population. We calculated the DSC values based on our proposed method and the competing methods for child and young adult group, respectively, shown in Figure 6. The DSC scores of our proposed method were higher than the DSC sores of all the comparison methods, demonstrating an increase of the subject consistency in the learned brain connectivity networks. Furthermore, the DSC values of the young adult group were always higher than those of the child group, suggesting an increase of consistency with increased age. In other words, the brain functional network continues to change and stabilize through adolescence, which is in line with the previous studies and common assumptions of brain development (Arain et al., 2013; Fair et al., 2009; Jolles et al., 2010).

FIGURE 6

Dice similarity coefficient (DSC) values for the proposed methods and the competing methods

Classification

To further illustrate our findings, we conducted classification experiments to separate different age groups with the learned brain network. Specifically, we considered the connectivity matrix jointly estimated for each subject as the input to a linear support vector machine (SVM) classifier (Schölkopf & Smola, 2002) to divide subjects into child or young adult groups. We also performed classification tasks with the constructed networks using graph Laplacian learning based Fourier transform, GFT and alternative GFT for comparison (Sandryhaila & Moura, 2013; Shuman et al., 2013; Wang et al., 2020).Sensitivity, specificity and accuracy were adopted as the criteria to evaluate the performance. We repeated the prediction experiments 10 times with 10‐fold cross‐validation to avoid occasionality and overfitting and displayed the result in Table 4.

TABLE 4

Classification Results for fMRI measurements

emoid‐fMRI
Method	ACC	SEN	SPE
MGFT	0.9606 ± 0.0022	0.9385 ± 0.0036	0.9805 ± 0.0021
GLFT	0.9401 ± 0.0080	0.9448 ± 0.0106	0.9349 ± 0.0085
GFT	0.9187 ± 0.0047	0.9230 ± 0.0078	0.9140 ± 0.0087
AGFT	0.9415 ± 0.0033	0.9491 ± 0.0064	0.9332 ± 0.0073
nback‐fMRI
Method	ACC	SEN	SPE
MGFT	0.9592 ± 0.0034	0.9378 ± 0.0038	0.9785 ± 0.0053
GLFT	0.9050 ± 0.0106	0.9331 ± 0.0071	0.8743 ± 0.0187
GFT	0.9066 ± 0.0079	0.9069 ± 0.0105	0.9067 ± 0.0104
AGFT	0.9237 ± 0.0057	0.9309 ± 0.0049	0.9164 ± 0.0107

Note: ACC, SEN, and SPE are the abbreviations of accuracy, sensitivity, and specificity, respectively. We display the result for mean value ± SD. MGFT, GLFT, and AGFT is the short for multi GFT analysis, graph Laplacian learning based Fourier transform and alternative GFT, respectively.

In summary, the proposed multigraph learning method outperformed the other GFT GSP methods. Our approach was able to extract useful shared latent information from multi‐paradigms, resulting in improved classification performance, which increased the reliability of the detected regions associated with brain maturation during puberty.

DISCUSSION

The most discriminative regions

In this section, we provided a detailed investigation of the most discriminative regions selected by our method to explore their relationships with brain maturation. The cingulate cortex was detected for both groups, which has been reported to be essential for integrating inputs from various sources (including cognitive and emotional networks) and processing emotional information. Previous work has provided evidence that cingulate cortex undergoes a long period of development from infancy to late childhood (Bush, Luu, & Posner, 2000).

Most cortical regions detected in both age groups were located in frontal lobes. The frontal lobes are among the last cortex to mature anatomically and functionally, and may develop until mid‐twenties or even later (Rubia et al., 2000). Specially, the middle frontal gyrus identified on the left hemispheric prefrontal region has a stronger generic activation in the adults group compared with the adolescents. A significant increase in power with age has been reported for both left middle inferior frontal gyrus and right inferior frontal gyrus (Rubia et al., 2000). The inferior frontal cortex (both in left and right cerebrum) has also been reported to be activated during the emotional tasks for both adolescents and young adults (Iidaka et al., 2002). The superior frontal gyrus (SFG) is detected for both groups, which is thought to contribute to a variety of cognitive and motor functions. The lateral part of SFG plays an especially important role in the working memory network (Boisgueheneuc et al., 2006; Li et al., 2013). The middle frontal gyrus has been demonstrated to be related to memory retrieval and is observed in both groups, which is consistent with the studies illustrating the co‐activation of inferior and middle frontal gyrus for both young and older adults (Gunning‐Dixon et al., 2003; Rajah, Languay, & Grady, 2011). The right paracentral lobule is observed only in the young adult group, implying a stronger activation compared with the younger age, which provides additional evidence of the study by Langan (Langan & Seidler, 2011).

In the parietal lobe, increased functional specialization related to memory had been reported for inferior parietal cortex and this maturation happens during puberty (Gogtay et al., 2004; Rivera, Reiss, Eckert, & Menon, 2005). The identified postcentral gyrus is believed to be the first region to mature (Gogtay et al., 2004).

In the temporal lobe, the inferior and middle temporal gyrus were observed in the child group. According to a longitudinal study conducted by Gogtay (Gogtay et al., 2004), the inferior temporal lobe, similar to the maturation pattern with inferior frontal lobe, matures early and does not change much thereafter. The middle temporal gyrus is thought to be one of the heteromodal association sites involved in integration of memory, audiovisual association and object‐recognition functions (Gogtay et al., 2004). Moreover, the stronger extensive activation in dealing with the emotional task has been compared with the older group (Gunning‐Dixon et al., 2003), which is in accordance with our findings in the temporal lobe.

Major findings and contributions

We propose a novel framework to incorporate multi‐paradigm fMRI data to detect biomarkers related to brain maturation. From the most frequently identified brain regions, our work is consistent with the common knowledge about brain maturation during adolescence. Especially the abovementioned middle, inferior, and superior frontal gyrus, inferior parietal cortex, and inferior and middle temporal gyrus were identified to have similar activation pattern associations with age in previous work. We also found newly identified regions such as claustrum and insula, which have not previously been reported to be associated with development. These regions deserve further investigation on their roles in the emotion identification or working memory task.

Apart from the brain region detection, our proposed model also demonstrated its superior performance in other aspects. The higher concentration of energy distributed over the low frequency band of the learned graph Laplacian (Figure 4) demonstrated better extraction of well‐organized brain networks over the frequently used kernel based methods. This indicates that different paradigms provide complementary information for better characterization of brain networks based on an appropriately designed model. In addition, the robustness of our model was also evaluated and a high classification performance was achieved (Figure 3). Finally, both the high accuracy of classifying different age groups and subject consistency provided further evidence regarding the validity of our findings.

Limitations and future directions

One potential limitation is the sample size and scope of paradigms (284 subjects and two paradigms). Thus, a future direction is the testing of our proposed model on larger number of data sets from different modalities and paradigms, with further applications to mental disorder diagnosis and prognosis, and other clinical applications. Recent research also shows that the parcellation strategies can affect the constructed functional networks and the predictive power (de Reus & Van den Heuvel, 2013; Pervaiz, Vidaurre, Woolrich, & Smith, 2020). Therefore, another interesting direction would be to evaluate the stability of the proposed model with different parcellation atlases. Other concerns include the test–retest reliability of task fMRI (Elliott et al., 2020), and the reliability measured by intraclass correlation coefficient is poor compared with the resting state fMRI. Opposite result was observed in a recent study where predictive models built from task fMRI outperformed those from resting‐state fMRI in the classification of intelligence levels (Greene, Gao, Scheinost, & Constable, 2018). Because of this, the combination of multi‐paradigm fMRI studies including both resting state and tasks may increase the reliability. Recent work by Gao et al. (Gao, Greene, Constable, & Scheinost, 2019) demonstrated improved classification performance with the combination of resting state and task fMRI, and the reliability varies with different models and different data sets (the conclusion was based on the experiments on both PNC and Human Connectome Project [HCP] [van Essen et al., 2013] data set). We also conducted experiments based on resting state fMRI in our study and most of the identified regions were the same between two age groups and within the default mode network, which suggests that the DMN has similar patterns from childhood through puberty into adulthood. Thus resting‐state fMRI has limited contribution in identifying discriminative regions related to brain maturation and the classification performance of using resting‐state fMRI was worse compared with task fMRI. Therefore, we do not include the result of using resting‐state fMRI while believe the combination of multi‐paradigm fMRI is a more promising research direction. Finally, we recently propose a hypergraph learning based framework to capture higher order relationships in the brain network (Xiao et al., 2019). Accordingly, the proposed model can be expanded by the use of hypergraph Laplacian matrix, which can characterize more sophisticated structure of the brain.

CONCLUSION

In this paper, we proposed a multigraph Fourier analysis framework to study brain functional connectivity differences between children and young adults. Specifically, fMRI data were parcellated into regions of interest (ROI) according to a predefined functional atlas (Power et al., 2011). The two tasks of fMRI observations were treated as graph signals defined on the ROIs within the brain. The brain connectivity networks were jointly estimated based on the proposed model for both paradigms, followed by graph Fourier analysis with filtering to decompose the original signals into several frequency bands. Hub regions were detected in the low frequency domain based on the biological significance and energy concentration. In both synthetic and real data analysis, our proposed framework demonstrated the superiority over other GSP tools such as GFT, AGFT, and GLFT for the construction of brain connectivity network by considering the shared latent graph structures from both paradigms. From the analysis results of PNC data set, we identified the hub regions that were activated in child and young adult groups and explored their implication for brain maturation during puberty. We verified our findings with several existing studies and interpreted the significance of the detected regions.

CONFLICT OF INTERESTS

The authors have no conflict of interest to disclose.

Appendix S1: Supporting information

Click here for additional data file.