Smoothed group information-guided independent component analysis for brain functional network analysis

Abstract

The present invention discloses a smoothed group information-guided independent component analysis method for brain functional network analysis, which belongs to the technical field of independent component analysis of brain images. The present invention includes performing independent component analysis to obtain components at the group level as reference signals; calculating voxel features to construct a graph regularization term; using voxel features and reference signals as guidance, using multi-objective functions to perform iterative solutions to estimate the independent components of individual subjects; and calculating the time series corresponding to each component in the individual subject based on the extracted components. The present invention overcomes the limitation of the group information-guided independent component analysis method currently widely used in the field of brain functional network extraction that does not optimize the smoothness of the extracted components, and obtains a more accurate brain functional network. Voxel features are introduced as a guide in the process of constructing the objective function, which enhances the spatial smoothness and functional correlation of the results, and can help the new method learn a network that is more in line with the actual working mechanism of the brain.

Description

Technical Field

[0001] The present invention belongs to the technical field of brain image independent component analysis, and in particular relates to a smoothed group information guided independent component analysis method for brain function network analysis. Background Art

[0002] Among matrix decomposition-based methods, independent component analysis (ICA) is a widely used method. The ICA method allows the extraction of signals of interest without any prior information about the task, that is, it allows us to perform operations on the data without any prior information. However, because data from different individuals have different time series, the components generated by ICA are randomly arranged. If we use ICA to process functional magnetic resonance imaging (fMRI) data from different subjects, there is no guarantee that the results will correspond one to one. Therefore, it is difficult to apply the ICA method to the processing of group-level fMRI data.

[0003] Researchers have proposed numerous approaches to address this problem, including group information-guided independent component analysis (GIG-ICA). Compared to other methods, GIG-ICA not only enables comparability across subjects but also explicitly optimizes the independence of extracted components and circumvents the difficulty of selecting threshold parameters. This method has been widely used in brain functional network analysis, but GIG-ICA does not optimize the smoothness of the extracted components. Due to the high noise content and high subject specificity of real-world data, the stable and accurate extraction of brain functional networks is extremely difficult, and there is still a significant gap between the extracted components and the actual components. A deeper understanding of brain function and its mechanisms, as well as accurate medical diagnosis of brain diseases, requires the extraction of more accurate and realistic brain functional networks. Summary of the Invention

[0004] In view of the limitation of the GIG-ICA method that does not optimize the smoothness of the extracted components, the present invention provides a more accurate and practical brain functional network, a new smooth group information guided independent component analysis method (Groupinformation guided smooth ICA, GIG-sICA) for brain functional network analysis.

[0005] In order to achieve the above object, the present invention adopts the following technical solutions:

[0006] The smoothed group information-guided independent component analysis method is used for brain functional network analysis, which includes the following steps:

[0007] Step S1, preprocessing functional magnetic resonance imaging (fMRI) data and representing the four-dimensional data as a two-dimensional matrix;

[0008] Step S2: concatenate data in the time direction and perform principal component analysis (PCA) dimensionality reduction on the data at the subject level and group level respectively;

[0009] Step S3, perform ICA on the dimensionality-reduced data to obtain independent components at the group level;

[0010] Step S4, calculating 16 voxel features to construct a graph regularization term;

[0011] Step S5, constructing a multi-objective function based on voxel features and composition as a reference, and normalizing the objective function;

[0012] Step S6, iteratively optimizing the data of each subject using a multi-objective function to obtain the independent components of the individual subject;

[0013] Step S7, calculating the time series corresponding to each independent component of the individual subject, that is, the corresponding activation pattern of the brain functional network.

[0014] Furthermore, the preprocessing specifically includes: removing brain function imaging data of the first few time points, and performing time layer correction, head motion correction, spatial normalization, and spatial smoothing on the brain function imaging data of the remaining time points.

[0015] Furthermore, the data representation refers to converting the four-dimensional fMRI data into a two-dimensional matrix; that is, pulling the three-dimensional image corresponding to each time point extracted from the original fMRI data into a row, and then concatenating these row vectors in the time direction to obtain a T*S matrix, where T is the number of time points and S is the number of voxels;

[0016] The concatenation in the time direction is based on the assumption that all subjects in the data have a common spatial structure, and the fMRI data of multiple subjects are connected along the time dimension.

[0017] Furthermore, the step S2 is to perform principal component analysis dimensionality reduction on the data at the subject level and the group level in series according to the time dimension. Specifically, the data matrix of each subject is subjected to PCA dimensionality reduction to n1 dimensions, and then these reduced matrices are connected in the time dimension to perform group-level PCA dimensionality reduction to n2 dimensions; wherein n2 must be the same as the final estimated number of components, and n1 can be greater than or equal to n2.

[0018] Furthermore, the independent component analysis method uses the FastICA algorithm or the Infomax algorithm.

[0019] Furthermore, the 16 voxel features calculated in step S4 can be divided into three groups, specifically including: the first group of features is the activation state of the voxel in gray matter, white matter, and cerebrospinal fluid; the mean of the activation state of the voxels in the neighborhood in gray matter, white matter, and cerebrospinal fluid; the mean of the similarity of the activation state of the voxel and the voxels in the neighborhood in gray matter, white matter, and cerebrospinal fluid; the second group of features is whether the voxel is at the edge of the template; the ratio of the voxels in the neighborhood at the edge of the template; the mean of the similarity of the voxel and the voxels in the neighborhood at the edge of the template; the third group of features is the degree of the voxel; the mean of the degree of the voxels in the neighborhood; the mean of the similarity of the degree of the voxel and the voxels in the neighborhood; the size of the connected area of ​​the voxel;

[0020] The calculation of the activation status of voxels in gray matter, white matter, and cerebrospinal fluid includes: judging whether the voxel is activated based on whether it is in the corresponding templates of the three, with activation being 1 and vice versa;

[0021] The calculation of the mean activation state of the voxels in the neighborhood in gray matter, white matter, and cerebrospinal fluid includes: the sum of the activation states of all voxels in the neighborhood divided by the number of voxels in the neighborhood;

[0022] The calculation of the mean similarity of the activation states of the voxels and their neighboring voxels in the gray matter, white matter, and cerebrospinal fluid includes: if the activation states of two voxels in the gray matter, white matter, or cerebrospinal fluid are both 0 or both 1, they are similar and their similarity value is 1; otherwise, the similarity value is 0, and finally the mean of these similarities is calculated;

[0023] The calculation of whether a voxel is on the edge of the template includes: taking the binary image template as input, performing a planar convolution on the image using the Sobel operator (two 3*3 matrix templates), obtaining vertical and horizontal difference approximations, and selecting the largest one as the grayscale size of the point. If the value is greater than the given threshold, the point is considered an edge point and the value is 1, otherwise the value is 0.

[0024] The calculation of the ratio of the voxels in the neighborhood to the edge of the template includes: the sum of the characteristic values ​​of whether the voxels in the neighborhood are at the edge of the brain is divided by the number of voxels in the neighborhood;

[0025] The calculation of the mean similarity of the voxel and the voxels in the neighborhood at the template edge includes: if the feature values ​​of the two voxels at the template edge are both 0 or both 1, they are similar, and the similarity value is 1; otherwise, the similarity value is 0, and finally the mean of these similarities is calculated;

[0026] The calculation of the voxel degree includes: taking the reciprocal of the absolute value of the difference between the voxel and the voxels in the neighborhood, multiplying it by the reciprocal of the distance between the two voxels, and then summing them;

[0027] Among them, the Z value, also called the standard score, is the process of dividing the difference between a number and the mean by the standard deviation. It is the deviation from the mean with the standard deviation as the unit, that is, the distance of a raw score from the mean is measured with the standard deviation as the ruler.

[0028] The calculation of the mean value of the voxel degree in the neighborhood includes: dividing the sum of the degrees of the voxels in the neighborhood by the number of voxels in the neighborhood;

[0029] The calculation of the mean value of the voxel degree in the neighborhood includes: dividing the sum of the degrees of the voxels in the neighborhood by the number of voxels in the neighborhood;

[0030] The calculation of the mean similarity between the voxel and the voxel degrees in the neighborhood includes: taking the difference between the voxel and the degrees of other voxels in the neighborhood, normalizing them by an exponential function with base e, and finally taking the inverse and calculating the mean;

[0031] The calculation of the size of the connected region includes: using the region growing algorithm to aggregate a voxel into a larger region; specifically, defining a variable and a maximum Z value distance max_dist. When the absolute value of the difference between the Z value of the voxel to be added and the mean Z value of the voxels in the segmented region is less than max_dist, the voxel is added to the segmented region. Otherwise, the algorithm stops; the size of the connected region obtained in the end is used as the eigenvalue.

[0032] Furthermore, in step S5, the multi-objective function is constructed as a reference by using voxel features and composition. The multi-objective optimization function is expressed as:

[0033]

[0034] The optimization of this multi-objective function includes the independence between different components of the test, the similarity between the components and the reference signal, and the smoothness of the components; is the estimated independent component, represents the random vector after whitening, v is a Gaussian variable with zero mean and unit variance, R i is a reference signal with zero mean and unit variance; G(.) represents any non-quadratic function, E(.) is the mathematical expectation, Tr(.) is the trace of the matrix, L represents the Laplace matrix, Represents Y i and R i The Pearson correlation coefficient is i and R i All have zero mean and unit variance, so this term can be expressed as Y i and R i expectations.

[0035] Among them, the similarity measure F(Y i ) includes calculating the mathematical expectation of the component and the reference product component, and expressing it in samples as the Pearson correlation coefficient between the component and the reference signal, so that the results can be comparable between different subjects; the independence measure J(Y i ) calculates the negative entropy of the component. The larger the negative entropy value, the stronger the independence of the result. The smoothness measure S(Y i ) is the trace of the product of the estimated component and the Laplacian matrix. By introducing the features of the voxels to construct the Laplacian matrix, a nearest neighbor graph is obtained to simulate its manifold structure, which ensures that spatially adjacent voxels or those with similar functions and features are located in the same functional network, thus enhancing spatial smoothness and correlation.

[0036] Among them, the specific method of constructing the Laplacian matrix L is: first calculate the voxel features to obtain a feature matrix of size S*F, where S is the number of voxels and F is the number of calculated voxel features; then standardize the feature matrix by column; for the first two groups of features, we construct the adjacency matrices Q1 and Q2 by multiplying the features within the group, and for the third group of features, we construct the adjacency matrix Q3 by calculating its Pearson correlation coefficient. The adjacency matrix Q used by the three groups of features at the same time is Q=Q1+Q2+Q3; finally, calculate the degree matrix D, and the values ​​of the diagonal elements of the matrix D are the column sums of the corresponding adjacency matrix Q; the Q matrix is ​​a symmetric matrix, and the D matrix is ​​a diagonal matrix; the Laplacian matrix is ​​expressed as: L=DQ.

[0037] Furthermore, the specific method for normalizing the objective function in step S5 includes normalizing using an inverse tangent function or a sigmoid function, and normalizing using an exponential function with e as the base; normalizing multiple objective functions is to avoid the optimization process being dominated by a larger objective function.

[0038] The normalized multi-objective optimization function is expressed as:

[0039]

[0040] The multi-objective function includes the independence constraints K(Y i ), the similarity constraint F(Y i ), the smoothness constraint of the component T(Y i );in, is the estimated independent component, represents the random vector after whitening, v is a Gaussian variable with zero mean and unit variance, R iis a reference signal with zero mean and unit variance; G(.) represents any non-quadratic function, E(.) is the mathematical expectation, Tr(.) is the trace of the matrix, L represents the Laplace matrix, Y i ' indicates Y i is the transpose of ; c and t are automatically determined parameters so that the three objective functions can have the same magnitude.

[0041]

[0042] t=-Tr(Y i *L*Y i ') / ln(F(Y i ))

[0043] Furthermore, the step S6 performs iterative optimization on the data of each subject using a multi-objective function. Specifically, the multi-objective function optimization is achieved using a linear weighted sum method. The total objective function of the optimization problem is expressed as:

[0044]

[0045] sta+b+d=1

[0046] Where a, b, and d are weight parameters. Although choosing a weight value can only find one point in the Pareto optimal set, if the weights are strictly positive and add up to 1 in the solution method of the multi-objective optimization problem, then for general multi-objective optimization problems, it is effective to optimize a linear weighted sum of the result of the cost function and its solution in a less complex problem. By adjusting the weights a, b, and d used in the linear weighted sum objective function, different points in the Pareto set can be explored.

[0047] Among them, the method for iterative solution of multi-objective functions is the gradient descent method.

[0048] Furthermore, the step S7 of calculating the time series corresponding to each independent component of the individual subject is as follows: the time series corresponding to each component in the individual subject is calculated based on the extracted components and the observed data, and the formula is: TC = E(X*Y -1 );TC is the time series, X is the two-dimensional data matrix, Y is the estimated independent component, Y -1 represents the inverse of Y; E(.) is the mathematical expectation; the estimated component Y represents the brain functional network, and the time series TC represents the activation state of the brain functional network.

[0049] Compared with the prior art, the present invention has the following advantages:

[0050] The present invention improves the GIG-ICA method and uses three objective functions to optimize the independent components of individual subjects. Its features and innovations mainly lie in: adding a graph regularization term for constraint, achieving a smoothing function on the basis of explicitly optimizing component independence, and overcoming the limitation of the GIG-ICA method that does not optimize the smoothness of the extracted components; secondly, the voxel characteristics are introduced as a guide in the construction of the objective function, making voxels connected in space or voxels with similar functions and characteristics more likely to be classified into the same functional network, thereby enhancing the spatial smoothness and correlation of the results. In addition, most of the noise in white matter, gray matter, cerebrospinal fluid and brain edges is removed, providing support for subsequent further research and analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 This is a flow chart of a smoothed group information-guided independent component analysis method provided by the present invention for brain functional network analysis.

[0052] Figure 2 It is a schematic diagram of the results obtained by a smoothed group information guided independent component analysis method and a GIG-ICA method provided by the present invention on simulated data with block noise added.

[0053] Figure 3 This is a schematic diagram of the results of a brain functional network obtained based on real data using a smoothed group information guided independent component analysis method and a GIG-ICA method provided by the present invention. DETAILED DESCRIPTION

[0054] In order to make the purpose, technical solutions and advantages of the present invention more clear, the present invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described are only used to explain the present invention, but are not used to limit the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary persons in this field without making creative work will fall within the scope of protection of the present invention.

[0055] Example 1

[0056] Refer to the attached Figure 1 The present invention provides a smoothed group information guided independent component analysis method for brain function network analysis, comprising:

[0057] Step S1, preprocessing the fMRI data and representing the four-dimensional data as a two-dimensional matrix;

[0058] Step S2: concatenate data in the time direction and perform PCA dimensionality reduction on the data at the subject level and group level respectively;

[0059] Step S3, perform ICA on the dimensionality-reduced data to obtain independent components at the group level;

[0060] Step S4, calculating 16 voxel features to construct a graph regularization term;

[0061] Step S5, constructing a multi-objective function based on voxel features and composition as a reference, and normalizing the objective function;

[0062] Step S6, iteratively optimizing the data of each subject using a multi-objective function to obtain the independent components of the individual subject;

[0063] Step S7, calculating the time series corresponding to each independent component of the individual subject, that is, the corresponding activation pattern of the brain functional network.

[0064] Among them, the steps of preprocessing the brain imaging data of the subjects include: removing the brain function imaging data of the first few time points, and performing time layer correction, head movement correction, spatial normalization, and spatial smoothing on the brain function imaging data of the remaining time points.

[0065] Data representation refers to converting four-dimensional fMRI data into a two-dimensional matrix. That is, the three-dimensional image corresponding to each time point extracted from the original fMRI data is pulled into a row, and then these row vectors are concatenated in the time direction to obtain a T*S matrix, where T is the number of time points and S is the number of voxels.

[0066] Among them, the concatenation in the time direction is based on the assumption that all subjects in the data have a common spatial structure, and the fMRI data of multiple subjects are connected along the time dimension.

[0067] The steps of performing subject-level and group-level PCA dimensionality reduction in series along the time dimension include: performing PCA dimensionality reduction on the data matrix of each subject to n1 dimensions, then concatenating these reduced matrices in the time dimension, and performing group-level PCA dimensionality reduction to n2 dimensions; where n2 must be the same as the final estimated number of components, and n1 can be greater than or equal to n2.

[0068] Among them, the ICA method uses the FastICA algorithm or the Infomax algorithm.

[0069] Among them, the 16 voxel features calculated can be divided into three groups, including: the first group of features is the activation state of the voxel in gray matter, white matter, and cerebrospinal fluid; the mean of the activation state of the voxels in the neighborhood in gray matter, white matter, and cerebrospinal fluid; the mean of the similarity of the activation state of the voxel and the voxels in the neighborhood in gray matter, white matter, and cerebrospinal fluid; the second group of features is whether the voxel is at the edge of the template; the ratio of the voxels in the neighborhood at the edge of the template; the mean of the similarity of the voxel and the voxels in the neighborhood at the edge of the template; the third group of features is the degree of the voxel; the mean of the degree of the voxel in the neighborhood; the mean of the similarity of the degree of the voxel and the voxels in the neighborhood; the size of the connected area of ​​the voxel;

[0070] The calculation of the brain, gray matter, white matter, and cerebrospinal fluid templates used in feature calculations involved setting thresholds of 0.5, 0.2, 0.9, and 0.7 for the initial templates of the brain, gray matter, white matter, and cerebrospinal fluid, respectively. Values ​​below the corresponding thresholds were set to 0, and otherwise were set to 1. The new templates were then registered based on the experimental data.

[0071] The calculation of the activation status of voxels in gray matter, white matter, and cerebrospinal fluid includes: judging whether the voxel is activated based on whether it is in the corresponding templates of the three, with activation being 1 and vice versa;

[0072] The calculation of the mean activation state of the voxels in the neighborhood in gray matter, white matter, and cerebrospinal fluid includes: the sum of the activation states of all voxels in the neighborhood divided by the number of voxels in the neighborhood;

[0073] The calculation of the mean similarity of the activation states of the voxels and their neighboring voxels in the gray matter, white matter, and cerebrospinal fluid includes: if the activation states of two voxels in the gray matter, white matter, or cerebrospinal fluid are both 0 or both 1, they are similar and their similarity value is 1; otherwise, the similarity value is 0, and finally the mean of these similarities is calculated;

[0074] The calculation of whether a voxel is on the edge of the template includes: taking the binary image template as input, performing a planar convolution on the image using the Sobel operator (two 3*3 matrix templates), obtaining vertical and horizontal difference approximations, and selecting the largest one as the grayscale size of the point. If the value is greater than the given threshold, the point is considered an edge point and the value is 1, otherwise the value is 0.

[0075] The calculation of the ratio of the voxels in the neighborhood to the edge of the template includes: the sum of the characteristic values ​​of whether the voxels in the neighborhood are at the edge of the brain is divided by the number of voxels in the neighborhood;

[0076] The calculation of the mean similarity of the voxel and the voxels in the neighborhood at the template edge includes: if the feature values ​​of the two voxels at the template edge are both 0 or both 1, they are similar, and the similarity value is 1; otherwise, the similarity value is 0, and finally the mean of these similarities is calculated;

[0077] The calculation of the voxel degree includes: taking the reciprocal of the absolute value of the difference between the voxel and the voxels in the neighborhood, multiplying it by the reciprocal of the distance between the two voxels, and then summing them;

[0078] Among them, the Z value, also called the standard score, is the process of dividing the difference between a number and the mean by the standard deviation. It is the deviation from the mean with the standard deviation as the unit, that is, the distance of a raw score from the mean is measured with the standard deviation as the ruler.

[0079] The calculation of the mean value of the voxel degree in the neighborhood includes: dividing the sum of the degrees of the voxels in the neighborhood by the number of voxels in the neighborhood;

[0080] The calculation of the mean similarity between the voxel and the voxel degrees in the neighborhood includes: taking the difference between the voxel and the degrees of other voxels in the neighborhood, normalizing them by an exponential function with base e, and finally taking the inverse and calculating the mean;

[0081] The calculation of the size of the connected region includes: using the region growing algorithm to aggregate a voxel into a larger region; specifically, defining a variable and a maximum Z value distance max_dist. When the absolute value of the difference between the Z value of the voxel to be added and the mean Z value of the voxels in the segmented region is less than max_dist, the voxel is added to the segmented region. Otherwise, the algorithm stops; the size of the connected region obtained in the end is used as the eigenvalue.

[0082] Among them, the multi-objective function of the reference structure is divided into voxel features and composition:

[0083]

[0084] The optimization of this multi-objective function includes the independence between different components of the test, the similarity between the components and the reference signal, and the smoothness of the components; is the estimated independent component, represents the random vector after whitening, v is a Gaussian variable with zero mean and unit variance, R i is a reference signal with zero mean and unit variance; G(.) represents any non-quadratic function, E(.) is the mathematical expectation, Tr(.) is the trace of the matrix, L represents the Laplace matrix, Represents Y i and R i The Pearson correlation coefficient is i and R i All have zero mean and unit variance, so this term can be expressed as Y i and R i expectations.

[0085] Among them, the specific method of constructing the Laplacian matrix L is: first calculate the voxel features to obtain a feature matrix of size S*F, where S is the number of voxels and F is the number of calculated voxel features; then standardize the feature matrix by column; for the first two groups of features, we construct the adjacency matrices Q1 and Q2 by multiplying the features within the group, and for the third group of features, we construct the adjacency matrix Q3 by calculating its Pearson correlation coefficient. The adjacency matrix Q used by the three groups of features at the same time is Q=Q1+Q2+Q3; finally, calculate the degree matrix D, and the values ​​of the diagonal elements of the matrix D are the column sums of the corresponding adjacency matrix Q; the Q matrix is ​​a symmetric matrix, and the D matrix is ​​a diagonal matrix; the Laplacian matrix is ​​expressed as: L=DQ.

[0086] In order to prevent the optimization process from being dominated by the larger objective function, the multi-objective function is normalized. This includes normalization using the inverse tangent function or the sigmoid function, or using an exponential function with base e. The normalized multi-objective optimization function is:

[0087]

[0088] The multi-objective function includes the independence constraints K(Y i ), the similarity constraint F(Y i ), the smoothness constraint of the component T(Y i );in, is the estimated independent component, represents the random vector after whitening, v is a Gaussian variable with zero mean and unit variance, R i is a reference signal with zero mean and unit variance; G(.) represents any non-quadratic function, E(.) is the mathematical expectation, Tr(.) is the trace of the matrix, L represents the Laplace matrix, Y i ' indicates Y i is the transpose of ; c and t are automatically determined parameters so that the three objective functions can have the same magnitude.

[0089]

[0090] t=-Tr(Y i *L*Y i ') / ln(F(Y i ))

[0091] The data of each subject is iteratively optimized using a multi-objective function. Specifically, the linear weighted sum method is used to implement the optimization problem of the multi-objective function. The overall objective function of the optimization problem can be expressed as:

[0092]

[0093] sta+b+d=1

[0094] Here, a, b, and d are weight parameters. Different values ​​of a, b, and d reflect the user's preference for component independence, similarity between a component and a reference signal, and smoothness of a component. In our experiments, we set a = 0.3, b = 0.4, and d = 0.3. In optimization problem solving methods, weights are strictly positive and sum to 1. Therefore, for general multi-objective optimization problems, optimizing a linear weighted sum of the cost function and its solution in a less complex problem is effective. By adjusting the weights a, b, and d used in the linear weighted sum objective function, different points in the Pareto set can be explored.

[0095] Among them, the method for iteratively solving the multi-objective function in the present invention is the gradient descent method.

[0096] The specific calculation of the time series corresponding to each independent component of the individual subject is as follows: the time series corresponding to each component in the individual subject is calculated based on the extracted components and the observed data, that is, TC = E(X*Y -1 );TC is the time series, X is the two-dimensional data matrix, Y is the estimated independent component, Y -1 represents the inverse of Y; E(.) is the mathematical expectation; the estimated component Y represents the brain functional network, and the time series TC represents the activation state of the brain functional network.

[0097] Data Generation for This Example: Simulated data was generated using the SimTB toolbox. Five subjects were simulated, each with 48x48 voxels, 150 time points, and eight independent components. Several random blocks of noise were then added to each component. The real data for this example comes from fMRI data collected from 25 healthy subjects at the FBIRN site.

[0098] Figure 2 These are the individual-level components estimated in the simulated data by this method. (A) shows the simulated data after adding noise; (B) shows the components estimated for each subject using the GIG-sICA method and the GIG-ICA method. The results show that the GIG-sICA method effectively removes noise from the simulated data.

[0099] Figure 3 These are some of the test results of our method on real data. (A) shows the five components; (B) shows the five components extracted by the GIG-ICA method; and (C) shows the five components extracted by the GIG-sICA method. Compared with the GIG-ICA method, the method proposed in this paper successfully removes most of the noise in gray matter, white matter, and cerebrospinal fluid.

[0100] Any matters not described in detail in this specification are prior art known to those skilled in the art. Although the above description of the present invention is based on specific embodiments to facilitate understanding of the present invention by those skilled in the art, it should be understood that the present invention is not limited to the scope of the specific embodiments. As long as various modifications are within the spirit and scope of the present invention as defined and determined by the appended claims, such modifications will be obvious to those skilled in the art, and all inventions and creations utilizing the concepts of the present invention are protected.

Claims

1. Smoothed group information guided independent component analysis method for brain functional network analysis, characterized by: The following steps are involved: Step S1, preprocessing the functional magnetic resonance imaging (fMRI) data and representing the four-dimensional data as a two-dimensional matrix; Step S2: concatenate data in the time direction and perform principal component analysis dimensionality reduction on the data at the subject level and group level respectively; Step S3, performing independent component analysis on the data after dimensionality reduction to obtain independent components at the group level; Step S4, calculating 16 voxel features to construct a graph regularization term; Step S5, constructing a multi-objective function by taking the voxel features and the independent components at the group level as references, and normalizing the objective function; Step S6, iteratively optimizing the data of each subject using a multi-objective function to obtain the independent components of the individual subject; Step S7, calculating the time series corresponding to each independent component of the individual subject, that is, the corresponding activation pattern of the brain functional network; In step S5, a multi-objective function is constructed by taking voxel features and independent components at the group level as references. The multi-objective optimization function is expressed as: The optimization of this multi-objective function includes the independence between different components of the test, the similarity between the components and the reference signal, and the smoothness of the components; is the estimated independent component, represents the random vector after whitening, v is a Gaussian variable with zero mean and unit variance, R i is a reference signal with zero mean and unit variance; G(.) represents any non-quadratic function, E(.) is the mathematical expectation, Tr(.) is the trace of the matrix, L represents the Laplace matrix, Represents Y i and R i The Pearson correlation coefficient is i and R i All have zero mean and unit variance, so this term is represented by Y i and R i expectations; Among them, the similarity measure F(Y i ) includes calculating the mathematical expectation of the component and the reference product component, and expressing it in samples as the Pearson correlation coefficient between the component and the reference signal, so that the results are comparable between different subjects; the independence measure J(Y i ) calculates the negative entropy of the component. The larger the negative entropy value, the stronger the independence of the result. The smoothness measure S(Y i ) is the trace of the product of the estimated component and the Laplacian matrix. By introducing the features of the voxels to construct the Laplacian matrix, a nearest neighbor graph is obtained to simulate its manifold structure.

2. The smoothed group information guided independent component analysis method according to claim 1 is used for brain functional network analysis, characterized in that: The preprocessing specifically includes: removing brain function imaging data of the first few time points, and performing time layer correction, head motion correction, spatial standardization, and spatial smoothing on the brain function imaging data of the remaining time points.

3. The method of smoothed group information guided independent component analysis according to claim 1 is used for brain functional network analysis, characterized in that: The data representation is to convert the four-dimensional fMRI data into a two-dimensional matrix; that is, the three-dimensional image corresponding to each time point extracted from the original fMRI data is pulled into a row, and then these row vectors are concatenated in the time direction to obtain a T*S matrix, where T is the number of time points and S is the number of voxels; The concatenation in the time direction is based on the assumption that all subjects in the data have a common spatial structure, and the fMRI data of multiple subjects are connected along the time dimension.

4. The method of smoothed group information guided independent component analysis according to claim 1 is used for brain functional network analysis, characterized in that: The step S2 is to perform principal component analysis dimensionality reduction on the data at the subject level and the group level in series in the time direction. Specifically, the data matrix of each subject is subjected to principal component analysis to reduce the dimension to n1 dimension, and then these reduced-dimensional matrices are connected in the time dimension to perform principal component analysis at the group level to reduce the dimension to n2 dimension; wherein n2 must be the same as the final estimated number of components, and n1 can be greater than or equal to n2.

5. The method of smoothed group information guided independent component analysis according to claim 1 is used for brain functional network analysis, characterized in that: The independent component analysis method uses the FastICA algorithm or the Infomax algorithm.

6. The method of smoothed group information guided independent component analysis according to claim 1 for brain functional network analysis, characterized in that: The 16 voxel features calculated in step S4 are divided into three groups, specifically including: the first group of features is the activation state of the voxel in gray matter, white matter, and cerebrospinal fluid; the mean of the activation state of the voxels in the neighborhood in gray matter, white matter, and cerebrospinal fluid; the mean of the similarity of the activation state of the voxel and the voxels in the neighborhood in gray matter, white matter, and cerebrospinal fluid; the second group of features is whether the voxel is at the edge of the template; the ratio of the voxels in the neighborhood at the edge of the template; the mean of the similarity of the voxel and the voxels in the neighborhood at the edge of the template; the third group of features is the degree of the voxel; the mean of the degree of the voxel in the neighborhood; the mean of the similarity of the degree of the voxel and the voxels in the neighborhood; the size of the connected area of ​​the voxel; The calculation of the activation status of voxels in gray matter, white matter, and cerebrospinal fluid includes: judging whether the voxel is activated based on whether it is in the corresponding templates of the three, with activation being 1 and vice versa; The calculation of the mean activation state of the voxels in the neighborhood in gray matter, white matter, and cerebrospinal fluid includes: the sum of the activation states of all voxels in the neighborhood divided by the number of voxels in the neighborhood; The calculation of the mean similarity of the activation states of the voxels and their neighboring voxels in the gray matter, white matter, and cerebrospinal fluid includes: if the activation states of two voxels in the gray matter, white matter, or cerebrospinal fluid are both 0 or both 1, they are similar and their similarity value is 1; otherwise, the similarity value is 0, and finally the mean of these similarities is calculated; The calculation of whether a voxel is on the edge of the template includes: taking the binary image template as input, performing a planar convolution on the image using the Sobel operator, obtaining vertical and horizontal difference approximations, and selecting the largest one as the grayscale size of the point. If the value is greater than the given threshold, the point is considered to be an edge point and the value is 1, otherwise the value is 0. The calculation of the ratio of the voxels in the neighborhood to the edge of the template includes: the sum of the characteristic values ​​of whether the voxels in the neighborhood are at the edge of the brain is divided by the number of voxels in the neighborhood; The calculation of the mean similarity of the voxel and the voxels in the neighborhood at the template edge includes: if the feature values ​​of the two voxels at the template edge are both 0 or both 1, they are similar, and the similarity value is 1; otherwise, the similarity value is 0, and finally the mean of these similarities is calculated; The calculation of the voxel degree includes: taking the reciprocal of the absolute value of the difference between the voxel and the voxels in the neighborhood, multiplying it by the reciprocal of the distance between the two voxels, and then summing them; Among them, the Z value, also called the standard score, is the process of dividing the difference between a number and the mean by the standard deviation. It is the deviation from the mean with the standard deviation as the unit, that is, the distance of a raw score from the mean is measured with the standard deviation as the ruler. The calculation of the mean value of the voxel degree in the neighborhood includes: dividing the sum of the degrees of the voxels in the neighborhood by the number of voxels in the neighborhood; The calculation of the mean similarity between the voxel and the voxel degrees in the neighborhood includes: taking the difference between the voxel and the degrees of other voxels in the neighborhood, normalizing them by an exponential function with base e, and finally taking the inverse and calculating the mean; The calculation of the size of the connected region includes: using the region growing algorithm to aggregate a voxel into a larger region; specifically, defining a variable and a maximum Z value distance max_dist. When the absolute value of the difference between the Z value of the voxel to be added and the mean Z value of the voxels in the segmented region is less than max_dist, the voxel is added to the segmented region. Otherwise, the algorithm stops; the size of the connected region obtained in the end is used as the eigenvalue.

7. The method of smoothed group information guided independent component analysis according to claim 6 is used for brain functional network analysis, characterized in that: The specific method of constructing the Laplacian matrix L is: first calculate the voxel features to obtain a feature matrix of size S*F, where S is the number of voxels and F is the number of calculated voxel features; then standardize the feature matrix by column; for the first two groups of features, we construct the adjacency matrices Q1 and Q2 by multiplying the features within the group. For the third group of features, we construct the adjacency matrix Q3 by calculating its Pearson correlation coefficient. The adjacency matrix Q used by the three groups of features at the same time is Q=Q1+Q2+Q3; finally, calculate the degree matrix D, and the values ​​of the diagonal elements of the matrix D are the column sums of the corresponding adjacency matrix Q; the Q matrix is ​​a symmetric matrix, and the D matrix is ​​a diagonal matrix; the Laplacian matrix is ​​expressed as: L=DQ.

8. The method of smoothed group information guided independent component analysis according to claim 1 is used for brain functional network analysis, characterized in that: The specific method of normalizing the objective function in step S5 includes normalizing using an inverse tangent function or a sigmoid function, or normalizing using an exponential function with e as the base; The normalized multi-objective optimization function is expressed as: The multi-objective function includes the independence constraints K(Y i ), the similarity constraint F(Y i ), the smoothness constraint of the component T(Y i );in, is the estimated independent component, represents the random vector after whitening, v is a Gaussian variable with zero mean and unit variance, R i is a reference signal with zero mean and unit variance; G(.) represents any non-quadratic function, E(.) is the mathematical expectation, Tr(.) is the trace of the matrix, L represents the Laplace matrix, Y i ′ represents Y i The transpose of ; c and t are automatically determined parameters so that the three objective functions have the same magnitude; t=-Tr(Y i *L*Y i ′) / ln(F(Y i ))。 9. The method of smoothed group information guided independent component analysis according to claim 8 is used for brain functional network analysis, characterized in that: The step S6 performs iterative optimization on the data of each subject using a multi-objective function. Specifically, the multi-objective function optimization is achieved using a linear weighted sum method. The total objective function of the optimization problem is expressed as: sta+b+d=1 Where a, b, and d are weight parameters. Although choosing a weight value can only find one point in the Pareto optimal set, if the weights are strictly positive and add up to 1 in the solution method of the multi-objective optimization problem, then for general multi-objective optimization problems, it is effective to optimize a linear weighted sum of the result of the cost function and its solution in a less complex problem. By adjusting the weights a, b, and d used in the linear weighted sum objective function, different points of the Pareto set can be explored. Among them, the method for iterative solution of multi-objective functions is the gradient descent method.

10. The method of smoothed group information guided independent component analysis according to claim 1 is used for brain functional network analysis, characterized in that: The step S7 of calculating the time series corresponding to each independent component of the individual subject is as follows: the time series corresponding to each component of the individual subject is calculated based on the extracted components and the observed data. The formula is: TC = E(X*Y -1 );TC is the time series, X is the two-dimensional data matrix, Y is the estimated independent component, Y -1 represents the inverse of Y; E(.) is the mathematical expectation; the estimated component Y represents the brain functional network, and the time series TC represents the activation state of the brain functional network.

Images