CN1294875C - Time series analysis method of nuclear magnetic resonance for brain functions based on constrained optimization - Google Patents

Abstract

The present invention relates to a time sequence analysis method of nuclear magnetic resonance for brain functions based on constrained optimization, which belongs to the field of a nuclear magnetic resonance technology. The present invention comprises: 1. a functional magnetic resonance time sequence is obtained; 2. the hemodynamics function of each pixel is estimated, and inverse convolution-operation is carried out to the time sequence obtained in step 1; 3. different irritating hemodynamics functions are estimated, and an optimization method is used according to the hemodynamics functions of each pixel estimated in step 2; 4. statistical hypothesis tests are carried out, and are carried out to each pixel one by one. The present invention is used for preoperative brain function positioning in medicine clinic, the diagnosis of brain diseases and evaluation of the brain diseases after cure, brain functional area positioning in brain science researches, and the function connecting analysis of brain functional areas.

Description

Brain Functional MRI Time Series Analysis Method Based on Constrained Optimization

technical field

The present invention relates to the field of nuclear magnetic resonance technology, in particular to a time-series analysis method of brain function nuclear magnetic resonance based on constraint optimization, which is used in clinical medicine for brain function positioning before surgery, brain disease diagnosis and prognosis evaluation, and brain science research The location of the brain functional area and the functional connection analysis of the brain functional area belong to the intelligent information processing technology.

Background technique

Since the birth of functional magnetic resonance imaging fMRI technology, fMRI time series analysis has always been a hot research direction that fMRI researchers from all over the world pay attention to. Generally, fMRI time series analysis algorithms can be divided into two categories: model-driven and data-driven. Due to the physiological rationality and ease of use of the data-driven method, it is gradually favored by neuroscientists from all over the world. The representative ones in the Model-driven method are the general linear model and the deconvolution model. In short, the general linear model is to add the prior knowledge of hemodynamics into the model by artificially specifying the design matrix, and then perform multiple regression analysis, so that the fitness of the prior model and fMRI data can be obtained. Its disadvantage is that the specification of the design matrix is subjective. The deconvolution model first obtains the convolution kernel through the deconvolution operation of the time series and the stimulus sequence, and then performs multiple regression analysis, that is, the design matrix is estimated. In summary, the general linear model assumes that different subjects have the same hemodynamic changes in different brain regions. The deconvolution model assumes that different pixels have different hemodynamic changes. In this sense, the deconvolution model is more in line with the physiological characteristics of the human brain. Compared with the general linear model, although the deconvolution model improves the sensitivity to a certain extent, studies have shown that the hemodynamic changes between each stimulus (trial) of the human brain are different. The product model is powerless.

Contents of the invention

The purpose of the present invention is to propose a new brain function nuclear magnetic resonance time series analysis method for the deficiencies in the prior art, which takes into account the hemodynamic inconsistency between each stimulation (trial) of the human brain, and then improves Accuracy of detection of brain functional activation areas. Based on the deconvolution technology and the constraint optimization method, combined with the statistical hypothesis test, the present invention makes full use of the brain function nuclear magnetic resonance time series information, and proposes a novel brain function nuclear magnetic resonance time series analysis method. Due to the constrained optimization method, combined with the latest progress in the study of cerebral hemodynamic response, the model itself can be self-expanded by adding new constraints to the model.

The brain functional nuclear magnetic resonance time series analysis algorithm based on optimization proposed by the present invention includes three basic steps of estimating the hemodynamic function of a single pixel, estimating the hemodynamic function of different stimuli, and statistical hypothesis testing:

1. Estimate the hemodynamic function of a single pixel

For each pixel, there is an accompanying time series. In this method, we believe that this time series contains three components: 1) hemodynamic signals derived from external stimuli; 2) drift caused by physiological activities such as breathing and heartbeat and magnetic resonance systems; 3) noise . We assume that the hemodynamic process is a linear system, that is,

Time series = stimulus sequencehemodynamic function + drift + noise

Among them,  represents the convolution operation. Using the deconvolution technique and the least square method, we can estimate the hemodynamic function corresponding to each pixel.

2. Estimate the hemodynamic function of different stimuli

Based on the research results of hemodynamics and the constrained optimization method, we hypothesize that different stimuli cause different hemodynamic responses. The formula for calculating the hemodynamic function for different stimuli is as follows:

$\underset{{\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}}{\mathrm{<span\; class="notranslate">\; min</span>}}{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; Z</span>}}<span\; class="notranslate">\; -</span>\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; J</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; \&CircleTimes;</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

st H _j ∈ N(h, ε)

是去卷积后的时间序列，N(h，ε)代表h的邻域，h是步骤1中求得的单个象素的血液动力学函数。基于上式的基本框架，我们可以加入另一个约束条件Among them, _Hj is the hemodynamic function of the jth stimulus, J is the total number of stimuli,

is the time series after deconvolution, N(h, ε) represents the neighborhood of h, and h is the hemodynamic function of a single pixel obtained in step 1. Based on the basic framework of the above formula, we can add another constraint

$\underset{{\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}}{\mathrm{<span\; class="notranslate">\; min</span>}}{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; Z</span>}}<span\; class="notranslate">\; -</span>\underset{\mathrm{<span\; class="notranslate">\; j</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; J</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; \&CircleTimes;</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}_{\mathrm{<span\; class="notranslate">\; 2</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

st H _j ∈ N(h, ε)

FWHM(H _i ): FWHM(H _j )=RT _i :RT _j , i≠j

Wherein, FWHM(H ₁ ) is the full width at half maximum of the hemodynamic function H ₁ , and RT _i is the response time of the i-th stimulus.

3. Statistical hypothesis testing

To determine whether a certain pixel is active, we perform a statistical hypothesis test.

Null hypothesis:

Alternative Hypothesis:

The statistic F is

$\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; =</span>\frac{\frac{\mathrm{<span\; class="notranslate">\; SB</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; SF</span>}}{{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; B</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; f</span>}}}}{\frac{\mathrm{<span\; class="notranslate">\; SF</span>}}{{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; f</span>}}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

Among them, H _min is the optimal solution of constrained optimization, $\mathrm{SB}={\left|\right|\hat{Z}\left|\right|}_2^2,\mathrm{SF}={\left|\right|\hat{Z}-\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}{f}_j\left(h_{\min }\right)\&CircleTimes;f_j\left|\right|}_2^2,$ d _B =NP-2, d _F =NP-2-(P+1). Under the null hypothesis, the statistic F obeys the F(d _B -d _F , d _F ) distribution, and a larger F indicates a greater possibility of activation of the corresponding pixel.

The present invention adopts the constraint optimization method, which can take into account the inconsistency of hemodynamic responses between different stimuli, and by adding new constraints, our method can be flexibly expanded, which is a simple and effective brain functional MRI Time Series Analysis Methods. The invention can be used for the location of brain function before operation in medical clinic, the diagnosis and prognosis evaluation of brain disease, the location of brain function area and the analysis of function connection of brain function area in brain science research.

Description of drawings

Fig. 1 is the schematic diagram of the brain function nuclear magnetic resonance time series analysis method based on constraint optimization of the present invention;

Fig. 2 and Fig. 3 are time series diagrams selected by the time series analysis method of brain function nuclear magnetic resonance based on constraint optimization in the present invention.

Detailed ways

In order to better understand the technical solution of the present invention, further description will be made below in conjunction with the accompanying drawings and specific embodiments.

The principle of the present invention based on the optimized brain function nuclear magnetic resonance time series analysis method is shown in FIG. 1 .

Step 1: Acquire fMRI time series. Acquisition of brain functional MRI time was completed on a magnetic resonance scanner equipped with an echo planar imaging (EPI) sequence. There are no special requirements for the specific parameters of the imaging, but generally there are no less than 3 layers, the sampling time points are generally dozens or more, and the spatial resolution is generally several millimeters, such as 3×3mm2.

Step 2: Estimate the hemodynamic function of a single pixel. Perform deconvolution operation on the time series obtained in step 1, and the result of deconvolution is the hemodynamic function of a single pixel.

Step 3: Estimate the hemodynamic function for different stimuli. The hemodynamic function of a single stimulus can be estimated by using the optimization method (formula (2)) according to the hemodynamic function of a single pixel estimated in step 2.

Step 4: Statistical hypothesis testing. A statistical hypothesis test (formula (3)) is performed on each pixel one by one to detect activated pixels.

In Figure 2, the selected time series is shown in Figure 2. There are 13 stimuli and 91 time points.

In Figure 3, the dotted line represents the original time series, the solid line represents the hemodynamic function of 13 stimuli, and the sharp dotted line below represents the time of stimulus presentation.

Example

1. Estimate the hemodynamic function of a single pixel

The selected time series is shown in Figure 2. There are 13 stimuli and 91 time points.

We first estimate the hemodynamic function of the pixel, the result is: [3.26 5.38 0.50 -3.92 -3.96 -4.46 -2.57]

2. Estimate the hemodynamic function of different stimuli

Using the hemodynamic function of the pixel estimated in the first step and constrained optimization (see equation (2)), we can obtain the hemodynamic function of each stimulus (see Fig. 3).

3. Statistical hypothesis testing

Using the formula (3), the value of the calculated F statistic is 18.53, obeying the F(7,195) distribution, and the corresponding probability value is 2.5618e-018. In general, if the p-value is taken as 0.01. Then this pixel is an active pixel. In addition, compared with the traditional deconvolution method, the following table is obtained: method F-statistic Traditional method (deconvolution) 16.94 The method of the invention 18.53

By comparison, the F statistic of the method of the present invention is 18.53, while the statistic of the traditional deconvolution method is 16.94. It can be seen that this method is superior to the traditional deconvolution method.

Claims (1)

1. brain function nuclear magnetic resonance, NMR Time series analysis method based on constrained optimization is characterized in that comprising the hemodynamics mathematic(al) function of estimating single pixel, hemodynamics mathematic(al) function and three basic steps of assumed statistical inspection of estimating different stimulated:

1) estimates the hemodynamics mathematic(al) function of single pixel: be to utilize contrary convolution technique under the hypothesis of a linear system in the hemodynamics variation process,, estimate the hemodynamics mathematic(al) function of each pixel correspondence by least square method;

2) the hemodynamics mathematic(al) function of estimation different stimulated: based on the basic framework of hematodinamics and constrained optimization method

$\underset{H_j}{\min }{\left|\right|\hat{Z}-\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}{H}_j\&CircleTimes;f_j\left|\right|}_2^2$

$s.t.H_j\&Element;N(h,\mathrm{\&epsiv;})$

(1)

Wherein, H _jBe j hemodynamics mathematic(al) function that stimulates, J is the total number that stimulates,

Be the time series after deconvoluting, (h ε) represents the neighborhood of h to N, and h is the hemodynamics mathematic(al) function of the single pixel of trying to achieve in the step 1;

3) assumed statistical inspection: for whether definite certain pixel activates, carry out assumed statistical inspection,

Null hypothesis:

Alternative hypothesis: Under null hypothesis, statistic F obeys F (d _B-d _F, d _F) distribute, and bigger F represents that the activated probability of corresponding pixel is big more.

Images