Group-oriented horizontal resting state fMRI human function brain region division method and electronic equipment

Abstract

The invention belongs to the field of neuroscience, and particularly relates to a group-oriented horizontal resting state fMRI human function brain region division method and electronic equipment. According to the method, an auto-encoder is adopted as a feature extractor, and compared with a traditional feature extraction algorithm, the method has high generalization ability and feature extraction ability, and the robustness of abnormal points in a feature space is high. Different from a current machine learning brain function partitioning algorithm, the method realizes adaptation of data sets of different scales by changing model parameters based on a deep learning technology. According to the method, the advantages of a traditional brain function partitioning algorithm are fully integrated, that is, global similarity and local gradient are considered at the same time, but the defects of the traditional brain function partitioning algorithm are avoided, that is, excessive induction bias is not introduced, so that the generated brain function partitioning has less hidden deviation.

Description

Technical Field

[0001] The present invention belongs to the field of neuroscience, and specifically relates to a method for dividing functional brain regions of human beings in a group-level resting-state fMRI and an electronic device. Background Art

[0002] The most common methods currently used are to use global similarity algorithms, local gradient algorithms, and probabilistic graphical model algorithms that combine global and local to divide human functional brain regions. The brain regions generated by the above three methods have high regional homogeneity, but there are still some disadvantages: ① The global similarity algorithm and the local gradient algorithm are sensitive to initialization; ② The global similarity algorithm tends to generate regions of the same size, thus affecting repeatability and physiological interpretability; ③ The probabilistic graphical model algorithm introduces too much inductive bias, which may induce implicit bias in brain region division. Summary of the invention

[0003] In view of the above problems, the present invention proposes a technical solution:

[0004] A method for group-level resting-state fMRI human functional brain area division, comprising the following steps:

[0005] S1: The human brain is divided into n voxels, where n is a natural number. The rs-fMRI signal corresponding to each voxel is input into the autoencoder model for training to obtain the trained autoencoder model. The autoencoder model includes an encoder and a decoder. The encoder in the trained autoencoder model is used as a feature extractor. The feature extractor is used to extract the feature vector of each voxel, where the feature vector of the i-th voxel is represented as h i ; When training the autoencoder model, the formula for the reconstruction loss is:

[0006]

[0007] Among them, L res Represents the reconstruction loss when training the autoencoder model;

[0008] r i Represents the input data corresponding to the i-th voxel in the autoencoder model;

[0009] Represents the output data corresponding to the i-th voxel in the autoencoder model;

[0010] S2: Add position coding to the feature vector of each voxel to obtain a feature vector containing position coding, where the feature vector containing position coding of the i-th voxel is represented as h' i ;

[0011] The clustering model is used to cluster the feature vectors containing position codes of all voxels. The number of clusters is determined by exhaustive method within the preset range to obtain the clustering results under each number of clusters. The clustering results under each number of clusters are tested using regional homogeneity test. The number of clusters corresponding to the clustering result with the best clustering effect is denoted as M. In the clustering result with the best clustering effect, the jth cluster center is denoted as u. j , where the value of j is 1, 2, 3, ..., M;

[0012] S3: Calculate the cluster distribution Q and auxiliary distribution P, where the cluster distribution Q is q ij The auxiliary distribution P is p ij The set of, the formula is:

[0013]

[0014] Among them, q ij represents the probability that the feature vector containing the position code of the i-th voxel belongs to the j-th cluster;

[0015] t represents the degrees of freedom of Student t distribution;

[0016]

[0017] f j =m j +v j +c

[0018]

[0019] Among them, p ij The auxiliary probability that the feature vector containing the position code of the i-th voxel belongs to the j-th cluster;

[0020] f j Represents the comprehensive information of soft cluster frequency information and voxel quantity information in the jth cluster;

[0021] m j represents the soft clustering frequency of the jth cluster, that is, the sum of the probabilities of all voxels being assigned to the jth cluster;

[0022] v j Represents the number of voxels in the jth cluster;

[0023] c represents control m j and v j The hyperparameters of the relative relationship between ;

[0024] N j represents the number of voxels corresponding to the feature vector containing position encoding assigned to the jth cluster;

[0025] S4: Use L as the total loss function to iteratively train the autoencoder model and clustering model, where the formula of L is:

[0026] L=L res +αL clu +βL dis

[0027] Among them, L represents the total loss function;

[0028] L clu represents the intra-class distance loss term; α represents the coefficient of the intra-class distance loss term;

[0029] L dis represents the inter-class loss term; β represents the inter-class loss term coefficient;

[0030]

[0031] Among them, p' ij represents the auxiliary probability that the feature vector containing the position code of the i-th voxel updated using the momentum method belongs to the j-th cluster;

[0032] When the first iteration, p' ij =p ij ;

[0033] In the second and subsequent iterations, in Indicates the previous iteration p' ij The value of

[0034]

[0035] S5: traverse all connected sub-regions in each cluster, and when there is more than one connected sub-region in a cluster, mark these connected sub-regions;

[0036] Find the connected subregion with the smallest area among the marked connected subregions as the first connected subregion, and calculate the connected subregion with the highest similarity among the connected subregions adjacent to the first connected subregion as the second connected subregion;

[0037] Merging the first connected subregion into the second connected subregion to form a merged connected subregion, and the cluster to which the merged connected subregion belongs is the cluster corresponding to the second connected subregion; and unmarking the merged connected subregion;

[0038] Repeat this step until each cluster has only one connected subregion, thus completing the brain region division;

[0039] The definition of a connected subregion is as follows: the surface of the cerebral cortex is modeled using a triangular mesh, with each voxel corresponding to a vertex of a triangle; within the same cluster, if a group of voxels are connected to each other through the edges of the triangle to form a closed area without breakpoints, then this part of the area is called a connected subregion;

[0040] The definition of connected subregions being adjacent is: if there is at least one voxel in each of the two connected subregions connected to each other through the edge of a triangle, then the two connected subregions are adjacent;

[0041] The method for comparing similarities between two connected sub-regions is:

[0042] The rs-fMRI signals of all voxels in the two connected sub-regions were averaged to obtain the rs-fMRI mean signal corresponding to each connected sub-region; the Pearson correlation coefficient of the fMRI mean signals of the two connected sub-regions was used to compare the similarities of the two connected sub-regions.

[0043] Preferably, h' i The calculation formula is:

[0044] h' i =h i ×cosx i ×cosy i ×cosz i

[0045] Among them, x i ,y i 、z i They are the three dimensional values ​​of the three-dimensional coordinates of the i-th voxel.

[0046] Preferably, the clustering model is a K-means++ algorithm.

[0047] Preferably, the method for testing regional homogeneity is a null model method.

[0048] Preferably, the value of t is 1.

[0049] Preferably, the value of c is 10.

[0050] Preferably, the values ​​of α and β are both 0.01.

[0051] The present invention also provides an electronic device,

[0052] include:

[0053] Memory, used to store computer programs;

[0054] A processor, for implementing the method steps described in any one of claims 1 to 7 when executing a program stored in a memory.

[0055] Beneficial effects:

[0056] (1) This method uses an autoencoder (AE) as a feature extractor. This algorithm is a deep learning method. Compared with traditional feature extraction algorithms, it has stronger generalization and feature extraction capabilities, and is more robust to outliers in the feature space.

[0057] (2) Different from the current machine learning brain function partitioning algorithm, this method is based on deep learning technology to adapt to data sets of different sizes by changing the model parameters.

[0058] (3) This method fully incorporates the advantages of the traditional brain function partitioning algorithm, that is, it considers both global similarity and local gradients, but avoids its disadvantages, that is, it does not introduce too much inductive bias, so the resulting brain function partitioning has less implicit bias. Specifically, improving global similarity is achieved by reducing the KL divergence of the distribution Q and P; reducing local gradients is achieved by increasing the difference in the cluster center vectors of each region; reducing inductive bias is achieved by optimizing the design of the autoencoder and this algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 is a flow chart of the method of the present invention;

[0060] Figure 2 It is the algorithm flow chart of the present invention; DETAILED DESCRIPTION

[0061] In order to make the purpose, features and advantages of the present invention more obvious and easy to understand, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the embodiments described below are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0062] like Figure 1 to Figure 2 As shown, a method for group-level resting-state fMRI (rs-fMRI) human functional brain area division includes the following steps:

[0063] S1: The human brain is divided into n voxels, where n is a natural number. The rs-fMRI signal corresponding to each voxel is input into the autoencoder model (AE) for training to obtain the trained autoencoder model. The autoencoder model includes an encoder and a decoder. The encoder in the trained autoencoder model is used as a feature extractor. The feature extractor is used to extract the feature vector of each voxel, where the feature vector of the i-th voxel is represented as hi ; When training the autoencoder model, the formula for the reconstruction loss is:

[0064]

[0065] Among them, L res Represents the reconstruction loss when training the autoencoder model;

[0066] r i Represents the input data corresponding to the i-th voxel in the autoencoder model;

[0067] Represents the output data corresponding to the i-th voxel in the autoencoder model;

[0068] Specifically, the encoder and decoder in the autoencoder model are both multi-layer perceptrons (MLPs), and the output of each layer of the feature extractor is expressed as:

[0069]

[0070] in, Represents the activation function. The optional activation function is the Relu function.

[0071] Represents the output data of the feature extractor layer l, where l is a natural number greater than or equal to 1;

[0072] Represents the input data of the feature extractor;

[0073] represents the weight matrix of the lth layer in the feature extractor;

[0074] represents the bias term of the lth layer in the feature extractor;

[0075] S2: Add position coding to the feature vector of each voxel to obtain a feature vector containing position coding, where the feature vector containing position coding of the i-th voxel is represented as h' i ; In one possible implementation, the formula is:

[0076] h' i =h i ×cosx i ×cosy i ×cosz i

[0077] Among them, x i ,y i 、z i are the three dimensional values ​​of the three-dimensional coordinates of the i-th voxel;

[0078] The clustering model is used to cluster the feature vectors containing position codes of all voxels. The number of clusters is determined by exhaustive method within the preset range to obtain the clustering results under each number of clusters. The clustering results under each number of clusters are tested using regional homogeneity test. The number of clusters corresponding to the clustering result with the best clustering effect is denoted as M. In the clustering result with the best clustering effect, the jth cluster center is denoted as u. j , where the value of j is 1, 2, 3, …, M.

[0079] In one possible implementation, the clustering model is a K-means++ algorithm; the method for testing regional homogeneity is a null model method;

[0080] The null model method is a common algorithm for testing regional homogeneity. The model consists of randomly rotating partitions of the same size, shape and relative position as the original partitions. This involves randomly rotating the original partitions around the x, y and z axes on a spherical dilation of the fs_LR_32k cortical surface. This process is repeated 1000 times to generate an average homogeneity distribution. We use the z-score to represent the functional homogeneity of the partition. The z-score is calculated as: [(original partition homogeneity - mean of the randomly rotated partition homogeneity) / standard deviation of the randomly rotated partition homogeneity]. The larger the z-score, the better the quality and more consistent connectivity pattern within the partition.

[0081] S3: Calculate the cluster distribution Q and auxiliary distribution P, where the cluster distribution Q is q ij The auxiliary distribution P is p ij The set of, the formula is:

[0082]

[0083] Among them, q ij represents the probability that the feature vector containing the position code of the i-th voxel belongs to the j-th cluster;

[0084] t represents the degree of freedom of Student t distribution. In this embodiment, the value of t is 1. The present invention calculates q ij When using Student t distribution as the kernel;

[0085]

[0086] f j =m j +v j +c

[0087]

[0088] Among them, p ijThe auxiliary probability that the feature vector containing the position code of the i-th voxel belongs to the j-th cluster;

[0089] f j Represents the comprehensive information of soft cluster frequency information and voxel quantity information in the jth cluster;

[0090] m j represents the soft clustering frequency of the jth cluster, that is, the sum of the probabilities of all voxels being assigned to the jth cluster;

[0091] v j Represents the number of voxels in the jth cluster;

[0092] c represents control m j and v j The hyperparameter of the relative effect relationship, in this embodiment, the value of c is 10;

[0093] N j represents the number of voxels corresponding to the feature vector containing position encoding assigned to the jth cluster;

[0094] S4: Use L as the total loss function to iteratively train the autoencoder model and clustering model, where the formula of L is:

[0095] L=L res +αL clu +βL dis

[0096] Among them, L represents the total loss function;

[0097] L clu represents the intra-class distance loss term; α represents the coefficient of the intra-class distance loss term;

[0098] L dis represents the inter-class loss term; β represents the inter-class loss term coefficient;

[0099] In this embodiment, the values ​​of α and β are both 0.01;

[0100]

[0101] Among them, p' ij represents the auxiliary probability that the feature vector containing the position code of the i-th voxel updated using the momentum method belongs to the j-th cluster;

[0102] When the first iteration, p' ij =p ij ;

[0103] In the second and subsequent iterations, in Indicates the previous iteration p'ij The value of

[0104] Considering that the drastic changes in the auxiliary distribution P are not conducive to model convergence, we use the momentum method to update the auxiliary distribution P in each iteration;

[0105]

[0106] S5: traverse all connected sub-regions in each cluster, and when there is more than one connected sub-region in a cluster, mark these connected sub-regions;

[0107] Find the connected subregion with the smallest area among the marked connected subregions as the first connected subregion, and calculate the connected subregion with the highest similarity among the connected subregions adjacent to the first connected subregion as the second connected subregion;

[0108] Merging the first connected subregion into the second connected subregion to form a merged connected subregion, and the cluster to which the merged connected subregion belongs is the cluster corresponding to the second connected subregion; and unmarking the merged connected subregion;

[0109] Repeat this step until each cluster has only one connected subregion, thus completing the brain region division;

[0110] The definition of a connected subregion is as follows: the surface of the cerebral cortex is modeled using a triangular mesh, with each voxel corresponding to a vertex of a triangle; within the same cluster, if a group of voxels are connected to each other through the edges of the triangle to form a closed area without breakpoints, then this part of the area is called a connected subregion;

[0111] The definition of connected subregions being adjacent is: if there is at least one voxel in each of the two connected subregions connected to each other through the edge of a triangle, then the two connected subregions are adjacent;

[0112] The method for comparing similarities between two connected sub-regions is:

[0113] The rs-fMRI signals of all voxels in the two connected sub-regions were averaged to obtain the rs-fMRI mean signal corresponding to each connected sub-region;

[0114] The Pearson correlation coefficient of the fMRI mean signals of two connected sub-regions was used to compare the similarity of the two connected sub-regions.

[0115] The pseudo code for this step is as follows:

[0116]

[0117]

[0118] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A method for group-level resting-state fMRI human functional brain area segmentation, characterized in that: The following steps are involved: S1: Divide the human brain into n voxels, where n is a natural number, and input the rs-fMRI signal corresponding to each voxel into the autoencoder model for training to obtain a trained autoencoder model; the autoencoder model includes an encoder and a decoder, and the encoder in the trained autoencoder model is used as a feature extractor; use Feature extractor, extracts the feature vector of each voxel, where the feature vector of the i-th voxel is represented as h i ; When training the autoencoder model, the formula for the reconstruction loss is: Among them, L res Represents the reconstruction loss when training the autoencoder model; r i Represents the input data corresponding to the i-th voxel in the autoencoder model; Represents the output data corresponding to the i-th voxel in the autoencoder model; S2: Add position coding to the feature vector of each voxel to obtain a feature vector containing position coding, where the feature vector containing position coding of the i-th voxel is represented as h' i ; The clustering model is used to cluster the feature vectors containing position codes of all voxels. The number of clusters is determined by exhaustive method within the preset range to obtain the clustering results under each number of clusters. The clustering results under each number of clusters are tested using regional homogeneity test. The number of clusters corresponding to the clustering result with the best clustering effect is denoted as M. In the clustering result with the best clustering effect, the jth cluster center is denoted as u. j , where the value of j is 1, 2, 3, ..., M; S3: Calculate the cluster distribution Q and auxiliary distribution P, where the cluster distribution Q is q ij The auxiliary distribution P is p ij The set of, the formula is: Among them, q ij represents the probability that the feature vector containing the position code of the i-th voxel belongs to the j-th cluster; t represents the degrees of freedom of Student t distribution; f j =m j +v j +c Among them, p ij The auxiliary probability that the feature vector containing the position code of the i-th voxel belongs to the j-th cluster; fj represents the comprehensive information of soft cluster frequency information and voxel quantity information in the jth cluster; m j represents the soft clustering frequency of the jth cluster, that is, the sum of the probabilities of all voxels being assigned to the jth cluster; v j Represents the number of voxels in the jth cluster; c represents control m j and v j The hyperparameters of the relative relationship between ; N j represents the number of voxels corresponding to the feature vector containing position encoding assigned to the jth cluster; S4: Use L as the total loss function to iteratively train the autoencoder model and clustering model, where the formula of L is: L=L res +αL clu +βL dis Among them, L represents the total loss function; L clu represents the intra-class distance loss term; α represents the coefficient of the intra-class distance loss term; L dis represents the inter-class loss term; β represents the inter-class loss term coefficient; Among them, p' ij represents the auxiliary probability that the feature vector containing the position code of the i-th voxel updated using the momentum method belongs to the j-th cluster; When the first iteration, p' ij =p ij ; In the second and subsequent iterations, in Indicates the previous iteration p' ij The value of S5: traverse all connected sub-regions in each cluster. When there is more than one connected sub-region in a cluster, mark these connected sub-regions. Find the connected subregion with the smallest area among the marked connected subregions as the first connected subregion, and calculate the connected subregion with the highest similarity among the connected subregions adjacent to the first connected subregion as the second connected subregion; Merging the first connected subregion into the second connected subregion to form a merged connected subregion, and the cluster to which the merged connected subregion belongs is the cluster corresponding to the second connected subregion; and unmarking the merged connected subregion; Repeat this step until each cluster has only one connected subregion, thus completing the brain region division; The definition of a connected subregion is as follows: the surface of the cerebral cortex is modeled using a triangular mesh, with each voxel corresponding to a vertex of a triangle; within the same cluster, if a group of voxels are connected to each other through the edges of the triangle to form a closed area without breakpoints, then this part of the area is called a connected subregion; The definition of connected subregions being adjacent is: if there is at least one voxel in each of the two connected subregions connected to each other through the edge of a triangle, then the two connected subregions are adjacent; The method for comparing similarities between two connected sub-regions is: The rs-fMRI signals of all voxels in the two connected sub-regions were averaged to obtain the rs-fMRI mean signal corresponding to each connected sub-region; the Pearson correlation coefficient of the fMRI mean signals of the two connected sub-regions was used to compare the similarities of the two connected sub-regions.

2. The method for group-level resting-state fMRI human functional brain region segmentation according to claim 1, characterized in that: h' i The calculation formula is: h’ i =h i ×cosx i ×cosy j ×cosz i Among them, x i ,y i 、z i They are the three dimensional values ​​of the three-dimensional coordinates of the i-th voxel.

3. The method for group-level resting-state fMRI human functional brain region segmentation according to claim 1, characterized in that: The clustering model is K-means++ algorithm.

4. The method for group-level resting-state fMRI human functional brain region segmentation according to claim 1, characterized in that: The method for regional homogeneity test is the null model method.

5. The method for group-level resting-state fMRI human functional brain region segmentation according to claim 1, characterized in that: The value of t is 1.

6. The method for group-level resting-state fMRI human functional brain region segmentation according to claim 1, characterized in that: The value of c is 10.

7. The method for group-level resting-state fMRI human functional brain region segmentation according to claim 1, characterized in that: The values ​​of α and β are both 0.

01.

8. An electronic device, characterized in that: include: Memory, used to store computer programs; A processor, for implementing the method steps described in any one of claims 1 to 7 when executing a program stored in a memory.

Images