Numerical Simulation of Higher-Order Nonlinearity of Human Brain Functional Connectivity Using Hypergraph p-Laplacian

Unravelling how the human brain structure gives rise to function is a central question in neuroscience and remains partially answered. Recent studies show that the graph Laplacian of the human brain’s structural connectivity (SC) plays a dominant role in shaping the pattern of resting-state functional connectivity (FC). The modeling of FC using the graph Laplacian of the brain’s SC is limited, owing to the sparseness of the Laplacian matrix. It is unable to model the negative functional correlations. We extended the graph Laplacian to the hypergraph p-Laplacian in order to describe better the nonlinear and high-order relations between SC and FC. First we estimated those possible links showing negative correlations between the brain areas shared across subjects by statistical analysis. Then we presented a hypergraph p-Laplacian model by embedding the two matrices referring to the sign of the correlations between the brain areas relying on the brain structural connectome. We tested the model on two experimental connectome datasets and evaluated the predicted FC by estimating its Pearson correlation with the empirical FC matrices. The results showed that the proposed diffusion model based on hypergraph p-Laplacian can predict functional correlations more accurately than the models using graph Laplacian as well as hypergraph Laplacian.


Introduction
Humans still know little regarding how the human brain works out cognitive tasks, due to its complex structure [1,2], in which hundreds of thousands of neurons, sets of neural populations and multiple brain regions are interconnected and interact together to yield diverse functions in the brain, both in the absence of external stimuli and performing cognitive tasks. Yet, the underlying mechanism behind this still remains incompletely understood. In the last decade, a growing body of research has been performed to explore the relationship between SC and FC with the help of two different kinds of invasive neuroimaging methods, diffusion magnetic resonance imaging (dMRI) measuring the fiber tracts of the white matter between brain regions [3,4] and functional magnetic resonance imaging (fMRI) recording the blood oxygenation level-dependent (BOLD) signals that characterize ongoing neural activities [5,6].
In the last decade, a great number of studies have been conducted to predict the human brain's SC using resting-state FC, such as network models on the basis of graph measures, neural mass models (NMMs) that describe the neural activity with partial differential equations, and learning-based regression models explicitly mapping FC from SC. The network models usually regard the brain as a graph in which the brain regions are treated as nodes and the connections between regions are defined as edges. A number of network metrics, such as node degrees [7,8], shortest path [9], as well as search information and path transitivity [10], etc., have been used to establish the relation between SC and resting-state FC. These models can only partially capture the nature of FC. Neural mass models [11][12][13] are used to mimic the interactions between neurons in the brain, without considering the signal exchanges between network modules, in which tens of physiological parameters need to be specified prior to obtaining the optimal solution. The regression learning models [14][15][16] simply relate SC with FC using weighted regression, in which the weights are trained using part of the empirical data. Nevertheless, the learned weight coefficients are only effective for the trained dataset while failing to extend to other datasets.
Some studies show that the graph Laplacian of the human brain's SC has been found to be able to highly capture the relationship between SC and FC, such as the graph diffusion (GD) model [17,18] proposed by Abdelnour, et al. and the hypergraph diffusion (HGD) model [19] proposed by us, in which the FC is described as an exponential decay function [20] of the Laplacian of SC. However, the GD model can yield merely the positive correlations in FC since the graph Laplacian reflects merely the direct connectivity between brain regions. The HGD model, defined by the hypergraph representation of SC, can capture more connection information than SC, thus showing better performance than GD model. However, the modeling of FC using the hypergraph Laplacian of the brain's SC is also limited to capturing the high-order and nonlinear relation between SC and FC, since the hypergraph Laplacian matrix is still sparse.
To overcome the limitations of GD and HGD models, we extend the graph Laplacian to hypergraph p-Laplacian, which better preserves the local structure [21] and has been successfully applied in pattern recognition and image processing [22,23]. The proposed hypergraph p-Laplacian diffusion model, termed as the HpGD model, is capable of better capturing the high-order relation between SC and FC, demonstrating better performance than GD and HGD models on simulating FC, including the modeling of the negative correlations.
The remainder of the paper is structured as follows. Section 2 gives the notations of the hypergraph network of the human brain adopted in this paper. In Section 3, the HGD model is first introduced in brief then the proposed HpGD model is presented in detail. Section 4 demonstrates comprehensive experiments on two empirical datasets and is following by some discussions in Section 5.

Hypergraph Network Notation of the Human Brain
The hypergraph representation of human brain is defined as G = (V, E, W), where V = {v i | i = 1, 2, · · · , n } is a set of n vertices referring to grey matter brain areas, E represents a set of m hyperedges, E = e j | j = 1, 2, · · · , m , in which each hyperedge is composed of several brain areas, v k (k = 1, 2, · · · , n), that connect to the corresponding brain area v j , and W = w j | j = 1, 2, · · · , m corresponds the weights of each hyperedge defined by summing all the edge weights in the hyperedge.
According to the above notations, we can define the incidence matrix of G by a |V| × |E| matrix H, in which H ij = h (v i , e j ) = 1 if v i is a vertex of hyperedge e j and 0 otherwise. The degree of a vertex v ∈ V in the hypergraph is defined as d v (v) = ∑ e∈E w (e) h (v, e), while the degree of a hyperedge e ∈ E in the hypergraph is defined by δ(e) = ∑ v∈V h (v, e). Let D v and D e denote the diagonal vertex degree matrix and hyperedge degree matrix, respectively, and W d indicates the diagonal matrix containing the weights of each hyperedge. Then the normalized hypergraph Laplacian can be expressed as for the case of simple graph when the hyperedge degree matrix D e reduces to 2I [24].

The Hypergraph Laplacian Diffusion Model
In what follows, we briefly describe the hypergraph diffusion model (HGD) in our recent work [19]. Firstly, consider an isolated cortical region v i . There are k(k = 1, 2, · · · , n) brain areas linking to the brain area v i via edges, with p(p = 1, 2, · · · , k) brain areas,v p k , with positive correlations and q(q = 1, 2, · · · , k) brain areas,v q k , with negative correlations to v i , respectively, where k = p + q. Then the first-order dynamics of the signal transmission between cortical regions v i and all the other brain regions (v j ) can be modified as where β is the decay rate. Assume that M p indicates the matrix whose entries are 1 if the correlations between brain areas are positive, otherwise 0, while M q indicates the matrix whose entries are 1 if the correlations between brain areas are negative, otherwise 0, Then the hypergraph network dynamics between all the brain regions can be depicted as where M p and M q can be estimated by analyzing the averaged experiment FC matrices, respectively. "•" denotes the Hadamard product. The solution to the above equation can be solved, i.e., is the normalized hypergraph network Laplacian. Suppose the initial condition x 0 = I, then the above solution can be used to simulate the functional interactions between brain areas.

The Hypergraph p-Laplacian Diffusion Model
Recently, the p-Laplacian has been found to be able to describe the nonlinear relationship in a graph and exhibit better performance than Laplacian in image classification [21]. To this end, we extend the Laplacian in the above model to p-Laplacian and the hypergraph p-Laplacian diffusion (HpGD) model can be formulated as: where L p denotes the hypergraph p-Laplacian, which satisfies, where g(·) is a real-valued function defined on the vertices of a weighted hypergraph G, and w h p i j ∈ W hp denotes the weight strength between vertex i and j in the hypergraph. The detailed calculation of the hypergraph p-Laplacian is presented in the following section.

The Computation of Hypergraph p-Laplacian
For any given p, it is complex and difficult to directly obtain the hypergraph p-Laplacian matrix L p from Equation (5). Here we propose an alternative approach to determine L p by eigen-decomposition.
Thomas et al. [21] proposed a method to define the eigenvectors and eigenvalues for the nonlinear operator L p . The real number λ p is named as an eigenvalue of the operator L p if the function f : V → R satisfies the relationship: where ϕ p is defined by ϕ p (x) = |x| p−1 sign(x), and the function f is called an eigenvector of L p corresponding to the eigenvalue λ p . When p = 2, the operator L p is equivalent to the standard hypergraph Laplacian L. According to Equation (5) Similarly, a functional F p : R V → R is defined as From the Theorem 3.1 in [21] we learned that the functional F p has a critical point at f ∈ R V if, and only if f is a p-eigenfunction of L p . The corresponding eigenvalue λ p is given as The above theorem can be applied to the computation of p-Laplacian. It is extremely time-consuming, which limits its applications. Here we propose an efficient approximation algorithm by converting the problem of complex solving eigenvectors into the following optimizing problem: min where Γ = ( f 1 , f 2 , · · · , f n ) denotes the optimal approximate eigenvector of L p . According to Equations (7) and (8), it is equivalent to: where Γ is the eigenvector matrix of the optimal approximate L p . Since the partial derivative with respect to f k v i can be written as, then the gradient change can be obtained from Equation (10), as follows, Meanwhile, the eigenvalues λ = (λ 1 , λ 2 , · · · , λ n ) can be obtained by Finally, the hypergraph p-Laplacian can be approximated by the eigenvalues and eigenvectors.
The estimation algorithm for hypergraph p-Laplacian is summarized in Algorithm 1, in which the step length α is chosen by α = β ∑ i k |Γ i k | ∑ i k |G i k | , where β is a parameter to control the step length.

Algorithm 1 The estimation of hypergraph p-Laplacian
Input: The structural connectivity network and three parameters k, p, β Output: Hypergraph p-Laplacian: L p Step 1: Initialize the hypergraph Laplacian matrix L and adjacent matrix W hp from the structural connectivity network.

Experiments
We demonstrate the performance of the extended HpGD model in contrast to the GD model and HGD model in predicting FC. Two commonly used experiment datasets, one with 90 regions of interest (ROIs) derived from eight subjects [17,18] and another with 246 ROIs derived from 50 subjects [25,26], are employed to test the models. The FC and SC of all subjects are rendered as two connectome matrices, with each element in the FC matrix representing the Pearson's correlation coefficients. Assume that X and Y are BOLD time series of two brain areas recorded using fMRI at rest, each time series has M points, then the Pearson correlation between the two brain areas is calculated by Note that all the p values of the correlations estimated in our approach are P << 10 −6 .

Correlations with the Experiment FC
For each subject, the prediction result is evaluated using the Pearson correlation between the predicted FC and the empirical FC. The higher of the correlation value is, the more similar the predicted FC is to the empirical FC. Figure 1 shows the time-varying Pearson correlations between the predicted FC and the empirical FC for each subject of the two datasets using the GD model, HGD model and the proposed HpGD model, respectively. We can observe that the maximum Pearson correlation of the eight subjects of the 90 ROI dataset changed from 0. It should be noted that the maximal correlations (critical points) appeared at βt crit = 2 for the 90 ROIs dataset, βt crit = 2.5 for the 246 ROIs dataset with the GD model, while at βt crit = 5 and βt crit = 5.5 with HGD model and HpGD model, respectively. Figure 2 demonstrates the scatter plots between the mean predicted FC and the mean empirical FC for the two datasets using the three models. We can observe that the proposed HpGD model showed the best effect in mapping the empirical FC because the relation between the predicted FC using HpGD model and the empirical FC was more linear.    Fig. 3 shows the mean predicted FC in contrast to the mean empirical FC rendered as matrix using the three models, respectively. Fig. 4 shows the mean predicted FC in  Figure 3 shows the mean predicted FC in contrast to the mean empirical FC rendered as a matrix using the three models, respectively. Figure 4 shows the mean predicted FC in contrast to the mean empirical FC rendered as a network [27] using the three models. It can be clearly observed that the FC predicted by the proposed HpGD model was much closer to the empirical FC than the other two models.

Stability and Robustness Analysis
To demonstrate that the predicted FC was not acquired by happenstance, for each dataset, we scrambled the mean SC 100 times randomly, then we estimated the FC for each permuted SC matrix with the three models. Figure 5 plots the histogram of the Pearson correlations between the averaged predicted FC and the averaged empirical FC. The histogram of the Pearson correlations without scrambling using the three models is also shown for the two datasets, respectively. The Pearson correlations were calculated between the mean empirical FC and the mean predicted FC (blue) as well as the mean predicted FC without permuting (orange). 14 correlations between the averaged predicted FC and the averaged empirical FC. The histogram of the Pearson correlations without scrambling using the three models are also shown for the two datasets, respectively. The Pearson correlations are calculated between the mean empirical FC and the mean predicted FC (blue) as well as the mean predicted FC without permuting (orange).

Conclusions
In this paper, we present HpGD as an extended graph diffusion model to map FC from SC relied on hypergraph p-Laplacian. HpGD allows the network to capture the non-linear and high-order relations between brain regions. Furthermore, two matrices

Conclusions
In this paper, we presented HpGD as an extended graph diffusion model to map FC from SC relying on hypergraph p-Laplacian. HpGD allowed the network to capture the nonlinear and high-order relations between the brain regions. Furthermore, the two matrices that indicate whether the correlations are positive or negative were embedded to entail the model simulating negative functional correlations. To validate the effectiveness of the HpGD model, we performed tests on two experiment connectome datasets with different resolutions. The experimental results from both the 90 ROI dataset and the 246 ROI dataset showed that a hypergraph p-Laplacian-based HpGD model could predict functional correlations more accurately than both the graph Laplacian-based GD model and the hypergraph Laplacian-based HGD model.
However, there are still some limitations. Firstly, the computation of p-Laplacian is time consuming and the three parameters k, p and β need to be optimized beforehand. Secondly, there is still a scarcity of biologically plausible interpretations to the mechanism behind the negative correlations. We will continue to address these issues in our future work.