CN103870710B - Tensor grouping method for multi-subject fMRI data analysis - Google Patents

Abstract

The invention discloses a grouping tensor method for multi-subject fMRI data analysis, which belongs to the field of fMRI data analysis. It is characterized in that, based on the principle of minimizing the difference between subjects, all the fMRI data of all subjects originally constructed as a large tensor are divided into multiple subgroups in units of subjects, so that the time course components of each subject in each subgroup There was the greatest cross-correlation between the spatial and brain-spatial activation regions. Compared with the original large tensor, the subgroup tensor constructed from the subgroup fMRI data is more compatible with the CP model of tensor decomposition, so it can decompose the time course components and brain space activation area components shared by multiple subjects with improved performance , if the correlation coefficient between the task-related time course component and the prior task stimulus process can be increased by about 0.1, the number of noise voxels in the task-related brain space activation area component can be reduced by about 23%, while the number of expected activation voxels remains basically unchanged. The invention has a reference function for solving the mismatch problem between other types of high-dimensional data and CP models.

Description

A grouped tensor method for multi-subject fMRI data analysis

technical field

The invention relates to an analysis method for multi-subject fMRI data, in particular to a tensor decomposition method for multi-subject fMRI data analysis.

Background technique

Functional magnetic resonance imaging (fMRI) is known as an effective window to observe the brain, because fMRI technology can capture the brain function of the subject when completing a specific task (task, such as vision, hearing, movement, etc.) data, and has the advantages of no damage and high spatial resolution. By using data-driven analysis methods such as blind source separation (BSS) and independent component analysis (ICA), it is possible to estimate multiple ( Usually dozens of) brain spatial activations (spatial activations) components and their time course (time courses) components provide detailed evidence for brain function analysis and clinical diagnosis.

fMRI data are high-dimensional data. Among them, single-subject fMRI data is 4-dimensional, such as 165×53×63×46, including 3-dimensional whole-brain data (53×63×46) and 1-dimensional scan times (165); multi-subject fMRI data is 5-dimensional , such as 165×53×63×46×16, that is, the number of 1-dimensional subjects (16) is added to the 4-dimensional single-subject fMRI data. Generally speaking, people will expand the 3D whole-brain data into one-dimensional voxels (voxels) data (53×63×46=153594), and the multi-subject fMRI data at this time also have as many as 3 dimensions (16×153594×165 ). Therefore, tensor decomposition methods suitable for processing high-dimensional (≥3-dimensional) data have great potential in the analysis of multi-subject fMRI data.

Currently, tensor-decomposed CP (canonical polyadic) models have been used for fMRI data analysis. The CP model assumes that the tensor data has a parallel factor structure. For multi-subject fMRI data expanded into 3 dimensions, it is assumed that multiple subjects share the same brain spatial activation area components and their time processes at different intensities. ingredients, as follows:

In the formula, Represents the tensor composed of all the multi-subject fMRI data, a _j and s _j represent the time course components and brain space activation area components shared between the subjects, J is the number of components, cj represents the intensity difference, ε represents the residual Difference. Among J time course components a ₁ , a ₂ ,…,a _J and J brain spatial activation area components s ₁ ,s ₂ ,…,s _J , there is a task-related ) components, which we denote as a _* and s _* . The following takes a _* and s _* as examples to illustrate the mismatch between multi-subject fMRI data and the CP model. Assuming a total of M subjects, the CP model assumes $a_{*}^{\left(1\right)}=a_{*}^2=...=a_{*}^{\left(m\right)},$ ${\mathrm{the\; s}}_{*}^{\left(1\right)}={\mathrm{the\; s}}_{*}^{\left(2\right)}=...={\mathrm{the\; s}}_{*}^{\left(m\right)}.$ However, existing studies have shown that even under the same task, the time-course components and brain-spatial activation components of different subjects will show large differences, resulting in $a_{*}^{\left(1\right)}\≠a_{*}^2\≠...\≠a_{*}^{\left(m\right)},$ ${\mathrm{the\; s}}_{*}^{\left(1\right)}\≠{\mathrm{the\; s}}_{*}^{\left(2\right)}\≠...\≠{\mathrm{the\; s}}_{*}^{\left(m\right)},$ In other words, there is a mismatch problem between the multi-subject fMRI data and the CP model, so the tensor decomposition method is difficult to exert its performance advantages in the current multi-subject fMRI data analysis.

Contents of the invention

The purpose of the present invention is to provide a solution to solve the mismatch between multi-subject fMRI data and CP model, minimize the difference between subjects, and improve the performance of tensor decomposition method in multi-subject fMRI data analysis.

The technical solution of the present invention is to divide a large tensor x constructed by M subjects' fMRI data into K subgroups in units of subjects Each subgroup contains N _k subjects, 2≤N _k <M; the sum of the mean values of cross-correlation coefficients between subjects of K subgroups It is the maximum value among the M grouping schemes, is the average value of the cross-correlation coefficients between the task-related time-course components and brain space activation area components of each subject in subgroup k, k=1,...,K; CP decomposition is performed on the tensors of each subgroup respectively:

The time course component a _j(k) and the brain space activation area component s _j(k) shared among the subjects of each subgroup, the intensity difference c _j(k) and the residual _(k) between the subjects were obtained. The basic block diagram is shown in Figure 1. Still taking task-related time course components and brain space activation area components as an example, since the time course components and brain space activation area components of each subject in subgroup k have the greatest correlation, then $a_{*\left(k\right)}^{\left(1\right)}\&ap;a_{*\left(k\right)}^{\left(2\right)}\&ap;...\&ap;a_{*\left(k\right)}^{\left(N_k\right)},$ ${\mathrm{the\; s}}_{*\left(k\right)}^{\left(1\right)}\&ap;{\mathrm{the\; s}}_{*\left(k\right)}^{\left(2\right)}\&ap;...\&ap;{\mathrm{the\; s}}_{*\left(k\right)}^{\left(N_k\right)},$ established. Visible, subgroup tensor data The degree of matching with the CP model is better than that of the original large tensor data In this way, the performance of the subgroup shared time course component and the brain space activation area component obtained by subgroup tensor decomposition will be improved. Through further analysis and processing, such as averaging, the performance can be significantly improved. For all The time-course component and the brain-spatial activation area component shared by the subjects.

The time course components and brain space activation area components of each subject's fMRI data need to be used when grouping, and the present invention uses the ICA method to separate and obtain them. As mentioned above, ICA can estimate task-specific brain spatial activation area components and their time-course components from fMRI data without any prior information. Specifically, the object of the cross-correlation coefficient calculation is the task-related time-course component of the fMRI data of M subjects and) Brain Spatial Activation Area Components ask for and) The covariance matrix R _a and R _s :

$R_a=\left[\begin{array}{cccc}1& r_{a_{*}^{\left(1\right)}{a}_{*}^{\left(2\right)}}& \&CenterDot;\&CenterDot;\&Center\; Dot;& r_{a_{*}^{\left(1\right)}{a}_{*}^{\left(m\right)}}\\ r_{a_{*}^{\left(2\right)}{a}_{*}^{\left(1\right)}}& 1& \&CenterDot;\&Center\; Dot;\&Center\; Dot;& r_{a_{*}^{\left(2\right)}{a}_{*}^{\left(m\right)}}\\ \&Center\; Dot;& \&Center\; Dot;& \&Center\; Dot;& \&Center\; Dot;\\ \&CenterDot;& \&CenterDot;& 1\&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ r_{a_{\text{*}}^{\left(m\right)}{a}_{*}^{\left(1\right)}}& r_{a_{*}^{\left(m\right)}{a}_{*}^{\left(2\right)}}& \&Center\; Dot;\&Center\; Dot;\&Center\; Dot;& 1\end{array}\right]$

Obtain the cross-correlation coefficient between different subjects p, q, p, q∈{1,2,...M} time course components and brain space activation area components In order to reduce the differences of time course components and brain space activation area components among subjects in the subgroup at the same time, a comprehensive matrix R=0.5|R _a |+0.5|R _s | was constructed from the covariance matrix R _a and R _s |, R The elements in are equal to Between and The mean value of the positive correlation coefficient between them.

The basic operation of grouping is to find the two most related subjects in terms of time course components and brain space activation area components, which are called two related subjects {p,q}. Set the elements of row p to zero (that is, find other subjects except p), find the number of rows corresponding to the maximum value in the elements of column p, record it as q, and then find the new subject q most related to subject p ;Repeat {p,q} to find (N _k -1) times, the first search makes the initial subject p=the smallest number among the ungrouped subjects, and then searches for p=q, then the construction is completed and contains Nk subjects a subgroup of . For example, N _k = 4, the initial subject of the first search is p, then through three {p,q} search: {p,q}→{q,u}→{u,v}, can construct { p, q, u, v} four subgroups of subjects. Let w be the smallest number among ungrouped subjects, then the initial subject of the next subgroup is p=w, and repeat the subgroup construction operation to construct a new subgroup, such as {w, x, y, z}. Repeat the construction of subgroups K times, and when K subgroups are constructed, a grouping scheme is constructed.

For K subgroups of a grouping scheme, the optimality of grouping cannot be guaranteed yet. Therefore, for the case of M subjects, the present invention constructs M kinds of grouping schemes. The specific method is to change the initial subjects for the first {p, q} search in the first subgroup to 1, 2, ..., M, and then calculate the average cross-correlation coefficient among K subgroups of each grouping scheme Sum is the average value of the cross-correlation coefficient between the task-related time course components and brain space activation area components of each subject in subgroup k, k=1,...,K, that is, to calculate the covariance matrix of K subgroups for each grouping scheme The sum of the average values of the triangular elements on R ₁ , R ₂ ,…,R _K α ₁ ,α ₂ ,…,α _m , and the grouping scheme with the largest sum value is the final grouping result.

The effect and benefits achieved by the present invention are that, compared with the existing method of using all multi-subject fMRI data to construct a tensor, the present invention can improve the matching degree of each subgroup fMRI data and the CP model, thereby significantly improving Performance of tensor decomposition methods. Taking the extraction of common task-related components of multi-subject fMRI data as an example, compared with the original large tensor method, the correlation coefficient between the task-related time process components obtained by the present invention and the prior task stimulus process is increased by about 0.1, However, the noise voxel number of the acquired task-related brain spatial activation area components decreased by about 23%, while the expected activation voxel number remained basically unchanged. In addition, this method has a reference role in solving the mismatch between other types of high-dimensional data (such as sensor array data, biomedical data, etc.) and the CP tensor model.

Description of drawings

Figure 1 is a basic block diagram of the present invention.

Fig. 2 is a flow chart of the present invention for analyzing multi-subject fMRI data.

Fig. 3 is a flowchart of task-related component cross-correlation calculation.

detailed description

A specific embodiment of the present invention will be described in detail below in conjunction with the technical scheme and accompanying drawings.

Assuming that the fMRI data of 16 subjects collected under the existing motion stimulation, the number of scans for each subject is 165. The present invention divides it into 4 subgroups, each subgroup contains 4 subjects, and then constitutes 4 subgroup tensors After tensor decomposition, the task-related components shared by each subgroup are respectively obtained, including time course components a _*(1) , a _*(2) , a _*(3) , a _*(4) and brain space activation area components s _*(1) ,s _*(2) ,s _*(3) ,s _*(4) . The specific steps are shown in Figure 2.

In the first step, the fMRI data of subject 1-subject 16 were subjected to PCA (principle component analysis) dimensionality reduction. According to the estimation results of the number of components of the motion stimulation fMRI data set in the existing literature, PCA was used to reduce the fMRI data of each subject from 165 dimensions to 20 dimensions, that is, it was assumed that the fMRI data contained 20 independent components.

In the second step, ICA is performed on the fMRI data of each subject after PCA dimensionality reduction. An information maximization algorithm with better performance can be used, and then the task-related components, including the time course component, can be selected from the 20 components separated by ICA and brain spatial activation area components

In the third step, the task-related time-course components and brain spatial activation area components 6 points) to perform cross-correlation calculations, the flow chart is shown in Figure 3:

(1) Obtain covariance matrix R _a and R _s , calculate the comprehensive matrix R=0.5|R _a |+0.5|R _s |.

(2) Let the grouping scheme pointer m=1.

(3) Let the subgroup pointer k=1, the initial subject p=m, repeat the search for two related subjects {p,q} for 3 times, and get {1,15}→{15,9}→{9,7}, Subgroup 1: {1,15,9,7} constructed.

(4) If the subgroup has not been constructed, that is, k<4, set k=k+1, and enter the next subgroup construction: let the initial subject p=the minimum number of ungrouped subjects, that is, set p=2, repeat two times Relevant subjects {p,q} searched for 3 subgroup construction operations to obtain subgroup 2: {2,11,10,12}. Similarly, construct subgroup 3: {3, 16, 6, 8} and subgroup 4: {4, 5, 13, 14}. So far, k=4, if k<4 is not satisfied, the mth grouping scheme has been constructed, including 4 subgroups: {1,15,9,7}, {2,11,10,12}, {3,16, 6,8}, {4,5,13,14}. Let m=m+1, enter the next grouping scheme construction.

(5) Determine whether m<16 is true. If so, repeat (3) and (4) to complete the next grouping scheme. For example, if the initial subject of the first subgroup (k=1) of the second grouping scheme is set as p=2, 4 new subgroups can be obtained: {2,15,9,7}, {1,4 ,11,10}, {3,16,6,8}, {5,12,13,14}. By analogy, until 16 grouping schemes are obtained.

(6) Calculate the sum of the average values of the cross-correlation coefficients among the 4 subgroups under the 16 grouping schemes is the average value of the cross-correlation coefficient between the task-related time course components and brain space activation area components of each subject in subgroup k, k=1,...,4, that is, to calculate the covariance matrix of 4 subgroups for each grouping scheme The sum of the average values of the triangular elements on R ₁ , R ₂ ,…,R _K α ₁ ,α ₂ ,…,α _m , and the grouping scheme with the largest sum value is the final grouping result. Here, the grouping scheme 7 with the largest sum value: {7,10,12,9}, {1,15,4,11}, {2,3,16,6}, {5,8,13,14} , as the final grouping result of the present invention.

The fourth step, according to the above grouping results, divide the fMRI data of 16 subjects after PCA dimensionality reduction into four subgroups {7, 10, 12, 9}, {1, 15, 4, 11}, {2, 3, 16 ,6}, {5,8,13,14}, and then construct 4 subgroup tensors

The fifth step is to use the tensor decomposition algorithm, such as the COMFAC (COMplex parallel FACtor analysis) algorithm, to decompose the 4 subgroup tensors to obtain 20 shared time process components and brain space activation area components, and select task-related components from them , including time course components a _*(1) ,a _*(2) ,a _*(3) ,a _*(4) and brain space activation area components s _*(1) ,s _*(2) ,s _{*(3 )} ,s _*(4) . If they are averaged separately, the task-related time-course components and brain-spatial activation area components common to all 16 subjects can be obtained.

Claims (1)

1. a kind of packet Tensor Method for how tested fmri data analysis, is characterized in that, by m tested fmri data The big tensor buildingIt is divided into k subgroup in units of tested Each subgroup comprises n_kIndividual tested, 2 ≤n_k<m；Tested cross-correlation coefficient mean value sum of k subgroupThis mean value sum M kind packet scheme is maximum,For tested task correlation time procedure component each in subgroup k and brain spatial activation Cross-correlation coefficient mean value between area's composition, k=1 ..., k；Each subgroup tensor is carried out respectively with cp decomposition:Obtain the shared time course composition a of tested of each subgroup_j(k)And Naokong Between active region ingredient s_j(k), the strength difference c of tested_j(k)And residual epsilon_(k)；J is the number of composition；

Cross-correlation coefficient calculate to as if m tested fmri data task correlation time procedure componentAnd active region composition between NaokongSeparated by ica method and obtain；Ask forWithCovariance matrix r_aAnd r_s, obtain different tested p, q, p, q ∈ { 1,2 ... m } Cross-correlation coefficient between the composition of active region between time course composition and NaokongWithCalculate synthetical matrix r= 0.5|r_a|+0.5|r_s|, carry out pairwise correlation tested { p, q } and find: started with initially tested p, by pth row element zero setting in r, Line number corresponding to maximizing in pth column element, is designated as q, then find the newly tested q mostly concerned with tested p；Repeat { p, q } finds n_k- 1 time, initially tested p=is not grouped tested middle numbering minimum to find order for the first time, finds afterwards and makes p=q, Then build and complete containing n_kAn individual tested subgroup；Repeat subgroup to build k time, build and complete k subgroup, then a kind of packet side Case builds and completes；

When m is tested, change in the first subgroup that { p, q } for the first time find initially tested be respectively 1,2 ..., m, Build altogether m kind packet scheme, then calculate every kind of packet scheme k subgroup covariance matrix r₁,r₂,…,r_kUpper triangle element Mean value sum α₁,α₂,…,α_m, the maximum packet scheme of value preset is final group result.