Sub-state Space Estimation Method Based on Functional Connectivity and Stochastic Process Model

Abstract

A sub-state space estimation method based on functional connectivity and stochastic process model, comprising the following steps: S1. Construct a stochastic process model, including an input module, a stochastic process module, and an output module; the input module receives functional connectivity data, the stochastic process module adopts an OU model, and the output module outputs a multi-dimensional OU sequence; S2. Calculate a functional connectivity matrix through BOLD data as the input of the input module; S3. Send the functional connectivity matrix into the stochastic process module, and output a stochastic process sequence through the output module; S4. Calculate dynamic functional connectivity through the stochastic process sequence; S5. Based on step S4, estimate the brain sub-state space to obtain sub-state probabilities and a state transition matrix. Compared with the existing calculation methods for sub-state space that are completely data-based, the present invention has stronger robustness and adaptability, helps to reveal the deep connection between functional connectivity and sub-state space, and provides more efficient technical support for the research of brain sub-states.

Description

Technical Field

[0001] The present invention belongs to the cross - technical field of biomedical engineering and brain science, and particularly relates to a sub - state space estimation method based on functional connectivity and stochastic process models. Background Art

[0002] Blood - oxygen - level dependent functional magnetic resonance imaging (BOLD - fMRI) is an imaging technique that indirectly detects neural activity by measuring changes in blood oxygen levels in the brain. It mainly includes two paradigms: the resting state and the task state. Among them, resting - state functional magnetic resonance imaging (rs - fMRI) highlights the rich structure of brain activity without tasks or stimuli. In the past few decades, the spontaneous fluctuations of brain functional activities in the resting state have attracted the attention of many researchers. The prevailing view in recent years is that the spontaneous fluctuation signals detected by rs - fMRI contain huge information on brain functional activities, have extremely important potential for data information mining, and can get closer to the essence of the organization and operation of the brain functional network through the analysis of spontaneous fluctuation information. When the brain completes a specific cognitive task, it is not responsible by a single brain region, but by the cooperation of multiple brain regions. In order to quantify the relationship between different brain regions, researchers have proposed the concept of functional connectivity (FC). It measures the synchrony of spontaneous activities in different brain regions of the brain in the resting state, usually obtained by calculating the Pearson correlation coefficient in a specific frequency band. In the past two decades, people have been committed to studying functional connectivity, that is, the functional interaction between different regions of the brain, which was considered static for a long time.

[0003] Previous studies on static functional connectivity technology assumed that the brain functional network was static and stable over a period of time. In fact, brain activities and networks are dynamically time-varying, and this dynamic nature is manifested in the resting state and different cognitive tasks. To address this issue, researchers proposed the dynamic functional connectivity method (dFC), which has been widely used to study the dynamic functional connectivity patterns between different brain regions. The dynamic nature can be characterized by the emergence and switching of probabilistic metastable substates (PMS). Past studies have used clustering analysis on dynamic functional connectivity or BOLD time series features, or constructed hidden Markov models, and found that these brain processes can be divided into several stable brain states, which correspond to the working modes of the brain under different cognitive tasks, emotional states, or physiological conditions. By calculating the occurrence probability of substates and the state transition matrix, the differences in substates among different populations can be understood, and the hierarchical structure of the brain state space in both spatial and temporal dimensions can be understood.

[0004] However, a series of past methods for calculating the substate space completely rely on the original BOLD sequence. In the case of data contamination or low quality in the collected BOLD sequence, especially when its temporal resolution is low or the number of time points is too short, it is often difficult to accurately calculate the substate space, which limits the application and comparison of substate space analysis across datasets and its promotion in clinical settings. At the same time, the relationship between resting-state functional connectivity and the substate space is not clear at present, which affects the joint analysis of data by researchers. Summary of the Invention

[0005] To solve the above technical problems, based on the basic mathematical properties of whole-brain functional connectivity and stochastic processes, and by using numerical simulation to solve the problem of the non-analyzable state of high-dimensional stochastic processes, the present invention designs an estimation method for the substate space based on functional connectivity and stochastic process models.

[0006] The object of the present invention is achieved through the following technical solutions: An estimation method for the substate space based on functional connectivity and stochastic process models, comprising the following steps:

[0007] S1. Construct a stochastic process model, including an input module, a stochastic process module, and an output module; the input module receives functional connectivity data as the initial input of the model, representing the covariance matrix of a multi-dimensional OU process; the stochastic process module adopts the OU model and numerically simulates to solve the numerical solution by the Euler method; the output module outputs a multi-dimensional OU sequence;

[0008] S2. Calculate the functional connectivity matrix from BOLD data as the input to the input module; generate the functional connectivity matrix FC by calculating the Pearson correlation coefficient between the BOLD signals of different brain regions.

[0009] S3. Feed the functional connectivity matrix obtained in step S2 into the stochastic process module and output a stochastic process sequence through the output module.

[0010] The discretized form of the stochastic process model is represented by the following difference equation:

[0011]

[0012] where Z t is an N-dimensional standard normal random variable:

[0013]

[0014] X t represents the stochastic process signal at time t, I N represents the identity matrix, θ is the time regression coefficient; ∑ is the covariance matrix, which is the functional connectivity matrix input by the input module; Δt is the same as the sampling interval TR of the real BOLD data.

[0015] Numerically solve the stochastic process model using the Euler method to obtain a multi-dimensional OU process as the output of the output module.

[0016] S4. Calculate the dynamic functional connectivity from the stochastic process sequence: further calculate the dynamic functional connectivity based on the multi-dimensional OU process output by the output module using the phase difference method, sliding window method, and time derivative multiplication.

[0017] S5. Based on step S4, estimate the brain sub-state space to obtain the sub-state probability and state transition matrix; for the calculation of brain sub-states, first reduce the dimension of the dynamic functional connectivity at each moment using the LEiDA method, reducing the three initial dynamic functional connectivities of each subject to two-dimensional data of N*T. This method mainly realizes the reduction by extracting the maximum eigenvector:

[0018] dFC(t) = VΛV -1

[0019] where dFC(t) is the initial matrix, i.e., the dynamic functional connectivity matrix at time t, i.e., dFC_pd(t,i,j), dFC_sw(t,i,j) or dFC_td(t,i,j); V is the eigenvector matrix, Λ is the eigenvalue vector, and LEiDA is the eigenvector with the largest eigenvalue.

[0020] For the two-dimensional data of N*T, take the data at each moment as a sample, perform clustering on it, and then use the Dunn coefficient as the evaluation index for clustering, and select the number of clusters with the largest Dunn coefficient; after completing the clustering operation, according to each cluster centroid, summarize the data with N features at each moment into one state, and a state sequence of length T can be obtained from the original N*T data;

[0021] Or directly extract the state from the sequence features through a hidden Markov model to obtain a state sequence of length T;

[0022] After the state is estimated by the clustering algorithm or the hidden Markov algorithm, a state sequence of length T is obtained, and the sub-state probability and the state transition matrix are calculated through the state sequence.

[0023] The beneficial effects of the present invention: Based on the whole-brain functional connectivity and the basic mathematical properties of the stochastic process, and at the same time solving the drawback that the state of the high-dimensional stochastic process cannot be analyzed through numerical simulation, a method for estimating the sub-state space based on functional connectivity and stochastic process models is designed, overcoming the drawback that the original sequence-based calculation method overly depends on the time resolution and sampling length of the original data, which is helpful for the application and comparison of sub-state space analysis across datasets and the popularization in clinical practice. Compared with the existing sub-state space calculation methods completely based on data, it has stronger robustness and adaptability, which is helpful to reveal the deep connection between functional connectivity and sub-state space, and provides more efficient technical support for the research of brain sub-states. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] Figure 1 It is a flowchart of the method for estimating the sub-state space based on functional connectivity and stochastic process models of the present invention.

[0025] Figure 2 It is a schematic diagram of the result of estimating the state occurrence probability in the embodiment of the present invention.

[0026] Figure 3 It is a schematic diagram of the result of the state transition matrix in the embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0027] The technical solution of the present invention will be further described below with reference to the accompanying drawings.

[0028] As Figure 1 shown, a method for estimating the sub-state space based on functional connectivity and stochastic process models of the present invention includes the following steps:

[0029] S1. Construct a stochastic process model, including an input module, a stochastic process module, and an output module. The stochastic process model adopts a multi-dimensional OU process model, and the time series obtained by numerical simulation through the OU process is used as the output.

[0030] The input module receives functional connectivity data as the initial input of the model, representing the covariance matrix of the multi-dimensional OU process. The stochastic process module adopts the OU model and numerically solves the numerical solution through the Euler method. The output module outputs the multi-dimensional OU sequence.

[0031] S2. Calculate the functional connectivity matrix through BOLD data as the input of the input module. By calculating the Pearson correlation coefficient between the BOLD signals of different brain regions, generate the functional connectivity matrix FC. For the BOLD signal time series of brain regions i and j, the calculation expression of the Pearson correlation coefficient r ij is as follows:

[0032]

[0033] where, T represents the total time length, B i,t represents the BOLD signal value of brain region i at time t, represents the average value of the BOLD signal of brain region i;

[0034] S3. Send the functional connectivity matrix obtained in step S2 into the stochastic process module, and output the stochastic process sequence through the output module.

[0035] The discretized form of the stochastic process model is represented by the following difference equation:

[0036]

[0037] where Z t is an N-dimensional standard normal random variable:

[0038]

[0039] X t represents the stochastic process signal at time t, I N represents the identity matrix (indicating that each normal random variable is independently and identically distributed); θ is the time regression coefficient (reflecting the speed of regression to the mean); ∑ is the covariance matrix, that is, the functional connectivity matrix input by the input module; Δt is the same as the sampling interval TR of the real BOLD data;

[0040] By numerically solving the stochastic process model using the Euler method, a multi-dimensional OU process is obtained as the output of the output module.

[0041] S4. Calculate the dynamic functional connectivity through a sequence of random processes: Based on the multi-dimensional OU processes output by the output module, further calculate the dynamic functional connectivity through methods such as the phase difference method, the sliding window method, and the time derivative multiplication method;

[0042] The dynamic functional connectivity under the phase difference definition is expressed as:

[0043] dFC_pd(t,i,j) = cos(angle o (t) - angle j (t))

[0044] where dFC_pd(t,i,j) represents the dynamic functional connectivity strength between brain regions i and j at time t. The closer it is to 1, the more synchronous the two brain regions are at that time; the closer it is to -1, the more anti-synchronous the two brain regions are at that time; and 0 indicates asynchronous;

[0045] angle i (t) represents the phase of the brain region at that time:

[0046]

[0047] where x(t) represents the sequence of random processes output by S3; is the Hilbert transform of x(t), defined as:

[0048]

[0049] The sliding window method divides a time series into T moments by translating the observation window. The dynamic functional connectivity at a certain moment is expressed as the Pearson correlation coefficient between the sequences under a certain window:

[0050]

[0051] where dFC_sw(t,i,j) represents the dynamic functional connectivity strength between brain regions i and j at time t, N represents the window size, B i,t1 represents the BOLD signal value of brain region i at time t1 within the observation window, represents the average value of the BOLD signals of brain region i within the observation window;

[0052] The time derivative multiplication method calculates the dynamic functional connectivity by calculating the product of the time derivatives at a certain moment and then performing a certain smoothing on the result:

[0053]

[0054] where dFC_td(t,i,j) represents the dynamic functional connectivity strength between brain regions i and j at time t, w represents the smoothing coefficient, and MTD(t,i,j) represents the product of time derivatives:

[0055]

[0056] where dt it is the derivative of the time series, and σ i is the standard deviation of dt in brain region i;

[0057] dt it = ts it - ts it-1

[0058] where ts it is the signal activity intensity of brain region i at time t;

[0059] S5. Based on step S4, estimate the brain sub-state space to obtain the sub-state probability and the state transition matrix; for the calculation of the brain sub-state, first reduce the dimension of the dynamic functional connectivity at each moment through the LEiDA method, and reduce the three initial dynamic functional connectivity strengths dFC (all three-dimensional data of N*N*T) of each subject to two-dimensional data of N*T. This method mainly realizes the reduction by extracting the maximum eigenvector:

[0060] dFC(t) = VΛV -1

[0061] where dFC(t) is the initial matrix, that is, the dynamic functional connectivity matrix at time t, namely dFC_pd(t,i,j), dFC_sw(t,i,j) or dFC_td(t,i,j); V is the eigenvector matrix, Λ is the eigenvalue vector, and LEiDA is the eigenvector with the largest eigenvalue;

[0062] For the two-dimensional data of N*T, take the data at each moment as a sample, cluster it through methods such as kmean++ or deep learning, and then use the dunn coefficient as the evaluation index for clustering. Select the number of clusters with the largest dunn coefficient. The calculation method of the dunn coefficient is:

[0063]

[0064] where C is the set of clusters, containing k clusters; dist(c i , c j ) is the minimum distance between clusters c i and c j ; Δ(C i ) is the maximum distance metric within the cluster C i ;

[0065] After completing the clustering operation, according to the centroid of each cluster, generalize the data with N features at each moment to a state, then a state sequence with a length of T is obtained from the original N*T data.

[0066] Alternatively, states are directly extracted from sequence features through a Hidden Markov Model to obtain a state sequence of length T; the Hidden Markov Model is represented as follows:

[0067] λ = (A, B, II)

[0068] where A represents the state transition probability matrix, B represents the observation state transition probability, and II represents the initial state probability;

[0069] After the states are estimated by the clustering algorithm or the Hidden Markov algorithm, a state sequence of length T is obtained, and the sub-state probability and the state transition matrix are calculated through the state sequence.

[0070] In this embodiment, the specific method for evaluating the estimation accuracy is as follows:

[0071] Estimation accuracy evaluation: Calculate the difference between the PMS space calculated based on this method and the PMS space completely based on real data to evaluate the fitting accuracy, that is, use the KL distance as a quantization index:

[0072]

[0073] where X is the domain of the probability distribution (i.e., the possible event space), and P(x) is the probability of event x under distribution P.

[0074] Figure 2 For the result display of the PMS space output based on this method and the comparison with the calculation result of real data.

[0075] Calculate the difference between the state transition matrix calculated based on this method and the state transition matrix completely based on real data to evaluate the fitting accuracy, that is, use the Markov entropy distance as a quantization index:

[0076]

[0077] where D MED (P1, P2) represents the Markov entropy distance between two state transition matrices, A represents the steady-state distribution probability, and P ij represents the probability from state i to state j.

[0078] Figure 3 For the result display of the state transition matrix output based on this method and the comparison with the calculation result of real data.

[0079] In summary, the method of the present invention aims to solve the problem that the sub-state space highly depends on the quality of real data. By only using functional connectivity and stochastic process models, the sub-state space can be effectively estimated. Based on the basic mathematical properties of whole-brain functional connectivity and stochastic processes, and at the same time, by means of numerical simulation, the disadvantage that the states of high-dimensional stochastic processes are not resolvable is solved. It overcomes the disadvantage that the original sequence-based calculation method overly depends on the time resolution and sampling length of the original data, which helps the application and comparison of sub-state space analysis across datasets, as well as its promotion in clinical practice. Compared with the existing calculation methods of sub-state space completely based on data, it has stronger robustness and adaptability, helps to reveal the deep connection between functional connectivity and sub-state space, and provides more efficient technical support for the research of brain sub-states.

[0080] Those of ordinary skill in the art will realize that the embodiments described herein are for helping the reader understand the principles of the present invention, and it should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those of ordinary skill in the art can make various other specific deformations and combinations without departing from the essence of the present invention according to these technical revelations disclosed in the present invention, and these deformations and combinations are still within the protection scope of the present invention.

Claims

1. A sub-state space estimation method based on functional connectivity and stochastic process models, characterized by: The following steps are involved: S1. Construct a random process model, including input module, random process module, and output module; The input module receives functional connectivity data as the initial input of the model, which represents the covariance matrix of the multidimensional OU process. The random process module uses the OU model to perform numerical simulation and solve the numerical solution through the Euler method. The output module outputs a multi-dimensional OU sequence; S2. Calculate the functional connectivity matrix using BOLD data as the input to the input module; generate the functional connectivity matrix FC by calculating the Pearson correlation coefficient between the BOLD signals of different brain regions. For the BOLD signal time series of brain regions i and j, the Pearson correlation coefficient r ij The calculation expression is as follows: Among them, T represents the overall time length, B i,t represents the BOLD signal value of brain area i at time t, represents the average value of BOLD signal in brain area i; S3, sending the functional connectivity matrix obtained in step S2 to the random process module, and outputting the random process sequence through the output module; The discretized form of the stochastic process model is expressed by the following difference equation: where Z t is an N-dimensional standard normal random variable: X t represents the random process signal at time t, I N represents the identity matrix, θ is the temporal regression coefficient; ∑ is the covariance matrix, i.e., the functional connectivity matrix of the input module; Δt is the same as the sampling interval TR of the real BOLD data; By numerically solving the random process model using the Euler method, a multidimensional OU process is obtained as the output module output; S4. Calculate dynamic functional connectivity through random process sequence: Based on the multidimensional OU process output by the output module, dynamic functional connectivity is further calculated through phase difference method, sliding window method, and time derivative multiplication; The dynamic functional connectivity under the phase difference definition is expressed as: dFC_pd(t,i,j)=cos(angle i (t)-angle j (t)) angle i (t) represents the phase of the brain region at that moment: where x i (t) represents the random process sequence of brain region i output by S3, is x i Hilbert transform of (t); The sliding window method divides a time series into T moments by shifting the observation window. The dynamic functional connectivity at a certain moment is expressed as the Pearson correlation coefficient between the sequences under a certain window: N represents the window size, B i,t1 represents the BOLD signal value of brain area i at time t1 within the observation window, represents the average value of the BOLD signal of brain area i within the observation window; The time derivative multiplication method calculates the dynamic functional connectivity by calculating the product of the time derivatives at a certain time and then smoothing the result: w represents the smoothing coefficient, MTD(t,i,j) represents the time derivative product; S5. Based on step S4, estimate the brain sub-state space and obtain the sub-state probability and state transition matrix. For the calculation of brain sub-states, first reduce the dimension of the dynamic functional connectivity at each moment, and reduce the three initial dynamic functional connectivity of each subject into N*T two-dimensional data. This method is mainly achieved by extracting the largest eigenvector: dFC(t)=VΛV -1 Where dFC(t) is the initial matrix, that is, the dynamic functional connectivity matrix at time t, namely dFC_pd(t,i,j), dFC_sw(t,i,j) or dFC_td(t,i,j); V is the eigenvector matrix, Λ is the eigenvalue vector, and the matrix after dimensionality reduction is the eigenvector matrix with the largest eigenvalue; For N*T two-dimensional data, the data at each moment is taken as a sample and clustered. Then, the Dunn coefficient is used as the clustering evaluation indicator, and the number of clusters with the largest Dunn coefficient is selected. After the clustering operation is completed, the data with N features at each moment is summarized into one state according to the centroid of each cluster. In this way, a state sequence of length T is obtained from the original N*T data. Or a hidden Markov model is used to extract states from sequence features to obtain a state sequence of length T; After the clustering algorithm or the hidden Markov algorithm estimates the state, a state sequence of length T is obtained, and the sub-state probability and the state transition matrix are calculated through the state sequence.

Images