CN103202692A - Quantitative determination method for brain functional connectivity frequency range - Google Patents

Abstract

The invention provides a quantitative determination method for a brain functional connectivity frequency range. According to the method, accurate frequency ranges of functional connectivity can be determined automatically. The method comprises the following steps: S1, obtaining brain function imaging data; S2, preprocessing the obtained brain function imaging data; S3, generating a space template of a brain interest function system according to prior information; S4, selecting certain region in the brain interest function system as a seed point, and computing correlations between all voxels on each frequency point and the seed point through a signal synchronism spectral analysis method; S5, performing space weighting summation on correlation results between all the voxels and the seed point, which are computed frequency by frequency, to obtain a functional connectivity frequency curve; and S6, performing significance determination on the functional connectivity frequency curve in the S5. According to the quantitative determination method, not only the functional connectivity frequency curve can be computed based on data, and the accurate functional connectivity frequency range can be determined automatically based on the curve, but also interference of global noise in brain function imaging data with the functional connectivity frequency analysis can be overcome.

Description

Brain function connects the quantitative judgement method of frequency range

Technical field

The present invention relates to neuroimaging and calculate the field, be specifically related to a kind of quantitative judgement method of the frequency range that connects at brain function based on the cerebral function imaging data.

Background technology

It is research field that more and more receives publicity in the neuroscience that function connects (Functional connectivity).The synchronicity of the cerebral activity that it is separated from each other by the metric space position is studied the mutual transmission of brain block information and the interactive contact of function.Existing discovering, similar to brain dissection network by the brain function network that function connects and composes, and the specific object (mainly being intensity) that function connects can reflect the compromised brain function that disease causes, embody the Changes of Plasticity of brain function integration mode in growth, study and the rehabilitation process, predict the difference of tested behavior.This shows that it is the movable interactional embodiment of brain neuroblastoma unit that function connects, and can reflect the function framework of brain inherence, has very important functional meaning.

Except the intensity attribute, evidence suggests that in recent years the frequency attribute that function connects also may have certain functional meaning.Though the interacting activity frequency between the neuron very wide in range (0.01-500Hz), according to nerve-blood vessel coupled relation (neurovascular coupling), the function connection that is embodied by blood flow signal often concentrates on a certain specific low-frequency range.And this frequency range may be relevant with the function in affiliated brain district.For example, the function of inside, elementary sensorimotor area connects the frequency range that concentrates on below the 0.1Hz, and the connection of the function of limbic system (between lesion of bilateral amygdala) is distributed in more wide in range frequency range (0-0.25Hz) relatively.Research about function connection frequency attribute now also is in the starting stage, ten pieces of lazy weights, and its concrete functional meaning also needs further exploration.

Existing decision method about function connection frequency range mainly contains two kinds.The one, according to the Fourier transformation rule, the function bonding strength that time domain is calculated decomposes frequency domain, passes through the coarse frequency scope of its concentration of energy of naked eyes subjective judgment again; The 2nd, researcher is divided into some segments with frequency band earlier, and computing function bonding strength in each frequency segment is relatively judged the frequency general scope that its intensity is higher relatively at last again.These two kinds of methods all depend on the subjective judgment of researcher, and objectivity is poor; And frequency resolution is coarse as a result in differentiation, possibly can't distinguish the careful difference that the difference in functionality systemic-function connects frequency.The more important thing is that existing two kinds of methods all can't effectively be avoided global noise (as the noise of large tracts of land existence on the spaces such as physiological noises such as breathing, heart beating, a moving noise) interference problem in the cerebral function imaging data (as function NMR (Nuclear Magnetic Resonance)-imaging fMRI, function near infrared imaging fNIRS).Function near infrared imaging particularly, it is as a kind of emerging non-invasive brain imaging technique, receives the concern that function connects researcher rapidly in the period of nearly three, be successfully applied to healthy people, paralytic, and baby's function integrate research.But because it is based on the ABSORPTION AND SCATTERING characteristic imaging that returns the photon of scalp surface from scalp surface through cerebral tissue, include the blood flow signal of a large amount of scalp parts in the data.Existing discovering, this noise be enough to upset in addition the damage study person to the analysis of cerebral tissue active signal.In addition in the actual tests process, because the fixing undesirable imaging auroral poles that causes of auroral poles medicated cap and the relative slip once in a while of scalp also can make data a moving noise of large-area burr or step occur, thus the carrying out of interfering data analysis.Because the physiological noise in the scalp blood flow and the frequency distribution of a moving noise are all very wide in range and complicated, be connected with function and probably have frequency alias, so these noises residual in above two kinds of frequency method of discrimination will probably influence function and connect the accuracy that frequency range is judged.Also there is not the pretreated method of reliable and effective data to address this problem now.

Summary of the invention

In order to overcome the subjective judgment that above-mentioned prior art all depends on researcher, objectivity is poor; And differentiate the coarse defective of frequency resolution as a result; Simultaneously, the more important thing is, all can't effectively avoid global noise (as the noise of large tracts of land existence on the spaces such as physiological noises such as breathings, heart beating, a moving noise) interference problem in the cerebral function imaging data (as function NMR (Nuclear Magnetic Resonance)-imaging fMRI, function near infrared imaging fNIRS) in order to overcome existing two kinds of methods.The present invention proposes the method that a kind of automatic discrimination function connects the precise frequency scope, comprise the steps: S1, obtain the cerebral function imaging data; S2, the described cerebral function imaging data of obtaining are carried out pretreatment; S3, generate the space template of brain function interested system according to prior information; S4, choose the intrasystem a certain zone of described brain function interested as seed points, adopt signal synchronicity frequency spectrum analysis method to calculate the dependency of all voxels and seed points on each Frequency point; S5, by frequency ground the correlation results of described all voxels that calculate and seed points is carried out the spatial weighting summation, the function that obtains described brain function interested system connects frequency curve; S6, function among the step S5 is connected frequency curve carry out significance by frequency and judge.

The present invention utilizes signal synchronicity frequency spectrum analysis methods such as coherence analysis, analyzes the dependency of all voxels and seed points on each Frequency point; And with in the analysis of spatial information pull-in frequency, eliminate global noise to the influence of synchronicity frequency analysis by the method for spatial weighting summation; Based on pursuing Frequency point to c ^SWCarry out statistical test, determine the precise frequency scope that function connects automatically.This method computational process is quick, can finish on common computer.The brain function that this method is mainly used in various metastable states (as quiescent condition) connects the frequency range exploratory study, is applicable to that also the brain function of other nonlinear states (as the influence of medicinal application to brain function) connects the analysis that frequency range changes.

Automatic discrimination function according to the present invention connects the method for precise frequency scope, it not only can go out function connection-curve of frequency distribution based on data computation, and determine accurate function automatically based on this curve and connect frequency range, can also overcome global noise in the cerebral function imaging data connects frequency analysis to function interference.

Description of drawings

Fig. 1 is the flow chart of implementation procedure of the present invention;

Fig. 2 is the sketch map of space-frequency connection matrix C;

Fig. 3 is the concrete auroral poles putting position of near infrared spectrum imaging sketch map in the instantiation;

Fig. 4 is the synchronicity spatial distribution map of different frequency in the instantiation;

Fig. 5 a is that the sensorimotor systemic-function of the method according to this invention gained connects frequency curve chart;

Fig. 5 b is that the sensorimotor systemic-function of traditional method gained connects frequency curve chart;

Fig. 6 a is spectrum distribution (left side) and the time graph example (right side) of the spontaneous fluctuation of simulation brain;

Fig. 6 b is the spatial distribution illustrated example of the spontaneous fluctuation of simulation brain;

Fig. 6 c is spectrum distribution (left side) and the time graph example (right side) of analog systems physiological noise;

Fig. 6 d be the analog systems physiological noise amplitude (on) and the spatial distribution illustrated example of phase place (descending);

Fig. 6 e is the time graph example of dummy head moving noise;

Fig. 6 f is the spatial distribution illustrated example of dummy head moving noise;

Fig. 7 a is under system's physiological noise disturbs, sequence example average time of region of interest 1 and region of interest 2;

Fig. 7 b is under system's physiological noise disturbs, the result that traditional coherence analysis method obtains;

Fig. 7 c is under system's physiological noise disturbs, space-frequency connection matrix example;

Fig. 7 d is under system's physiological noise disturbs, the result that method of the present invention obtains;

Fig. 7 e is under system's physiological noise disturbs, the ROC curve comparison diagram of traditional coherent analytical method and the inventive method;

Fig. 8 a is under a moving noise disturbs, the result that the traditional coherent method obtains;

Fig. 8 b is under a moving noise disturbs, the result that method of the present invention obtains;

Fig. 8 c is under a moving noise disturbs, the ROC curve comparison diagram of traditional coherent analytical method and the inventive method;

Fig. 9 a is that the inventive method is estimated excessive robust analysis result to the space template scope;

Fig. 9 b is that the inventive method is estimated too small robust analysis result to the space template scope;

Fig. 9 c is that the inventive method is to the robust analysis result of space template location estimation skew.

The specific embodiment

For making the purpose, technical solutions and advantages of the present invention clearer, below in conjunction with specific embodiment, and with reference to accompanying drawing, the present invention is described in further detail.

The theoretical premise of the method according to this invention is based on the existing research conclusion of following this area: movable similarity height between the voxel of (1) brain function internal system, and movable similarity is low between the inner voxel of function system and the outside voxel; (2) global noise can make the brain function between all voxels connect whole the rising; (3) on the main frequency that function connects, the global noise undercapacity in the actual tests is connected the difference of function system inside and outside to flood function fully.

The flow chart of Fig. 1, implementation procedure of the present invention.As shown in Figure 1, the present invention proposes the method that a kind of automatic discrimination function connects the precise frequency scope, it not only can go out function connection-curve of frequency distribution based on data computation, and determine accurate function automatically based on this curve and connect frequency range, can also overcome global noise in the cerebral function imaging data connects frequency analysis to function interference.

Referring to Fig. 1, after collecting the brain function data, choose the interior a certain voxel/zone of function system interested as seed points, calculate the synchronicity spatial distribution map of each voxel and seed points on each Frequency point, the space template of the function system interested of recycling priori, synchronicity spatial distribution map on each Frequency point is carried out the spatial weighting summation, obtain the synchronicity coefficient of spatial weighting, constitute function and connect frequency curve c ^SW, pursue Frequency point at last to c ^SWCarry out statistical test, the function that the significant frequency of intensity is brain function interested system connects frequency range.

The automatic discrimination function of describing the present invention's proposition with reference to the accompanying drawings in detail connects the method for precise frequency scope.Fig. 1 is the flow chart of implementation procedure of the present invention.As shown in Figure 1, the method according to this invention comprises the steps:

S1: obtain the cerebral function imaging data

The cerebral function imaging data can be the brain function nmr imaging datas according to the present invention, also can be brain function near infrared spectrum imaging datas.

Being captured on the magnetic resonance scanner that possesses plane echo-wave imaging (EPI) sequence of brain function nmr imaging data finished.Imaging parameters arranges as follows, and spatial resolution is generally several millimeters, as 3*3mm ², the repetition time, (repetition time TR) was generally the several seconds, and as 2s, the sampling time, the imaging scope enough covered functional areas interested and perimeter, partial function district preferably more than 5 minutes.

Being captured on the function near infrared spectrum imager of brain function near infrared spectrum imaging data finished.Imaging parameters arranges as follows:, temporal resolution is generally tens of or hundreds of milliseconds, and as 100ms, the auroral poles spacing is generally 3 ~ 5cm, sampling time, the imaging scope enough covered outside cortex part and the perimeter, partial function district of functional areas interested preferably more than 5 minutes.

Above gatherer process if be quiescent condition, requires experimenter's peace and quiet, and keeps head still as far as possible.

S2: pretreatment

After finishing, data acquisition needs to carry out conventional pretreatment.The pretreatment of function nuclear magnetic resonance data comprises: head moves and corrects, goes linear drift, space filtering (being also referred to as space smoothing), space criteriaization etc.The pretreatment of function near-infrared data comprises: resolve, go linear drift etc. by light intensity data to the blood oxygen concentration data.But above-mentioned these processes will determine according to concrete application target, except the resolving of near-infrared light intensity data, and nonessential process.

S3: the space template that generates brain function interested system according to prior information:

Determine the space voxel location of functional area interested according to prior information, generate the space template with the function system interested of brain function Data Matching.Prior information can be from the result of study of prior art, the also anatomic information that can provide from the structure picture scanning that the experimenter adds, or the analysis result of functional localization task scan.For example, function system interested is the sensorimotor system, and the experimenter is carried out the brain function scanning of finger motion task, and the brain district that this task is significantly activated is labeled as 1, and the zone marker that does not have significantly to activate is 0, can generate sensorimotor system space template m.According to concrete operating position, m can be carried out the space smoothing of suitable yardstick (as 8mm).

S4: choose seed points, computer memory-frequency connection matrix C generates the synchronicity spatial distribution map by frequency.

Choose the interior a certain voxel of function system interested or zone as seed points, be seed points time series s with the time series of this voxel or sequence definition average time that should the zone, calculate the time series x of each voxel in the image-generating unit and the synchronicity spectrum distribution c of s _Xs(λ), with c _Xs(λ) as column vector, be arranged in proper order together according to the locus of voxel, get final product the span-frequency connection matrix C.Fig. 2 is the sketch map of space-frequency connection matrix C, and as shown in Figure 2, the i of C is capable of Frequency point λ _iThe synchronicity size that goes up each voxel and seed points constitutes.According to this row vector of locus reconstruct, can obtain Frequency point λ _iSynchronicity spatial distribution map c _s(λ _i), shown in Fig. 2 c.

Calculate the synchronicity spectrum distribution and can adopt relevant (coherence) analytical method, also can adopt the method for Pearson correlation coefficient spectral factorization.

Wherein, relevant (coherence) analytical method is specific as follows:

$c_{\mathrm{xs}}\left(\mathrm{\λ}\right)=\frac{{\left|G_{\mathrm{xs}}\left(\mathrm{\λ}\right)\right|}^2}{G_{\mathrm{xx}}\left(\mathrm{\λ}\right)G_{\mathrm{ss}}\left(\mathrm{\λ}\right)}---\left(1\right)$

C in the formula _XsBe on Frequency point λ (λ), the coherence factor between time series x and the s, its value is 0 ~ 1,0 to be illustrated between this frequency time series and not have linear relationship, has desirable linear relationship between two time serieses of 1 expression.G _XsCrosspower spectrum between express time sequence x and the s, G _XxAnd G _SsThe auto-power spectrum of difference express time sequence x and s:

$G_{\mathrm{xs}}\left(\mathrm{\λ}\right)=F_x\left(\mathrm{\λ}\right)F_s^{*}\left(\mathrm{\λ}\right)$

$G_{\mathrm{xx}}\left(\mathrm{\λ}\right)=F_x\left(\mathrm{\λ}\right)F_x^{*}\left(\mathrm{\λ}\right)---\left(2\right)$

$G_{\mathrm{ss}}\left(\mathrm{\λ}\right)=F_s\left(\mathrm{\λ}\right)F_s^{*}\left(\mathrm{\λ}\right)$

In the formula, F _xAnd F _sThe frequency spectrum that obtains of express time sequence x and the analysis of s Fourier spectrum respectively, * represents conjugation.

The method of Pearson correlation coefficient spectral factorization is specific as follows:

Seasonal effect in time series counted is designated as for example N=300 of N(), Pearson correlation coefficient is calculated as follows:

$c_{\mathrm{xs}}=\frac{\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(t\right)s\left(t\right)}{\sqrt{\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}{x}^2\left(t\right)\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}{s}^2\left(t\right)}}---\left(3\right)$

According to discrete fourier inverse transformation rule:

$x\left(t\right)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{F}_x\left({\mathrm{\λ}}_i\right)e^{i2\mathrm{\π}{\mathrm{\λ}}_{i}t/N}---\left(4\right)$

$s\left(t\right)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{F}_s\left({\mathrm{\λ}}_i\right)e^{i2\mathrm{\π}{\mathrm{\λ}}_{i}t/N}---\left(5\right)$

Can be with c _XsDecompose on the frequency domain:

$c_{\mathrm{xs}}\left(\mathrm{\λ}\right)=\frac{N[\mathrm{Re}\left(F_x\left(\mathrm{\λ}\right)\right)\mathrm{Re}\left(F_s\left(\mathrm{\λ}\right)\right)+\mathrm{Im}\left(F_x\left(\mathrm{\λ}\right)\right)\mathrm{Im}\left(F_s\left(\mathrm{\λ}\right)\right)]}{\sqrt{\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}{x}^2\left(t\right)\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}{s}^2\left(t\right)}}---\left(6\right)$

S5: the synchronicity spatial distribution map is carried out the spatial weighting summation by frequency ground:

Utilize the space template p of the brain function interested system that generates among the step S3, in conjunction with each Frequency point λ _iLast synchronicity spatial distribution map c _s(λ _i) characteristics, calculate the spatial weighting coefficient w of each voxel x correspondence on this Frequency point _x(λ _i)

$w_x\left({\mathrm{\λ}}_i\right)=\frac{m_x-\stackrel{\‾}{m}}{\sqrt{\underset{x=1}{\overset{M}{\mathrm{\Σ}}}{(m_x-\stackrel{\‾}{m})}^2\underset{x=1}{\overset{M}{\mathrm{\Σ}}}{[c_{\mathrm{xs}}\left({\mathrm{\λ}}_i\right)-{\stackrel{\‾}{c}}_s\left({\mathrm{\λ}}_i\right)]}^2}}---\left(7\right)$

M is the total number of voxel (for example 52) in the formula,

Be the space average number of m,

Be c _sThe space average number.

Utilize w _x, calculate synchronicity spatial distribution map c by frequency ground _sSpatial weighting and coefficient, that is:

$c^{\mathrm{SW}}\left({\mathrm{\λ}}_i\right)=\underset{x=1}{\overset{M}{\mathrm{\Σ}}}{w}_x\left({\mathrm{\λ}}_i\right)\×c_{\mathrm{xs}}\left({\mathrm{\λ}}_i\right)---\left(8\right)$

The c of gained ^SWThe function that is function system interested connects frequency curve.

S6: function is connected frequency curve carry out the significance judgement by frequency:

Error type I in the statistical test (" abandoning true ") adopts non-parametric multiple comparisons method here.For the null hypothesis that makes up frequency domain function bonding strength distributes, seed points seasonal effect in time series power spectrum is fixed, its phase frequency spectrum of change at random, obtain seed points time series behind the phase randomization by the Fourier inversion rule, the time series of other voxels is not done any variation, calculate seed points seasonal effect in time series space-frequency synchronization matrix behind all voxels and the phase randomization according to the method for top step S4, and corresponding function connects frequency curve.With randomisation process repeated several times (as 100 times), utilize the value in all functions connection frequency curves that obtain at random, the null hypothesis that generates frequency domain function bonding strength distributes.Pursue the significance level p that Frequency point ground calculates the function connection frequency curve value of function system interested, setting threshold p ₀, wherein, 0＜p ₀＜1, optional, as p ₀=0.05, at p＜p ₀The time, significant frequency is the main frequency that function connects.

Hereinafter provide and adopted the method according to this invention that function is connected the detailed process that the precise frequency scope is differentiated example automatically.

The function of sensorimotor system connects frequency range under this example planning studies quiescent condition.Choose 16 experimenters and participate in experiment, 19 ~ 26 years old, male 8 people, women 8 people.All experimenters are without any motion or nervous system disease.Through the test of Edinburg handedness scale, all experimenters are dextromanuality (standard is greater than 40).Two sections brain function near infrared spectrums that all experimenters are scheduled to scan.Be 11 minutes quiescent condition the last period, requires experimenter's peace and quiet to loosen and close one's eyes, and keeps head still as far as possible, systematically do not think deeply.The back is both hands finger sequence motion task for one section.In carrying out brain function near infrared spectrum scanning process, at first be 63 seconds rest (seeing "+"), be task block then, then be baseline block.Have 7 task block and 8 baseline block, the length of each block can be from 20 seconds to 28 seconds picked at random, but each task block is identical with length with its corresponding baseline block.At task block, " ~ * * * ~ " sequence of 0.5s at first appears on the screen, and this sequence disappears then, has four alphabetical sequences that " F " and " J " random order is formed, for example " FJFJ ".Require the experimenter see after this alphabetical sequence with the forefinger of left hand and the right hand not only fast but also accurate by the FJFJ on the lower keyboard, order and number of times all require identical with alphabetical sequence on the screen.1.5s back alphabetical sequence can disappear, occurring next group then stimulates.The setting sweep parameter is: temporal resolution is 10Hz, 17 emission auroral poles, and 16 receive auroral poles, form 52 altogether and lead, and cover other zones such as bilateral sensorimotor cortex and part superior temporal gyrus.Concrete auroral poles putting position as shown in Figure 3.

The function near infrared spectrum data that collects is carried out following pretreatment, comprising: light intensity data removes linear drift to the conversion of brain oxygen density data.The unstable data of 30s before then both hands finger sequence motion task data being removed, high-pass filtering (0.017Hz), by leading generalized linear model (GLM) analysis of carrying out, carry out single-sample t-test in the population of subjects level at last, significance level p＜0.05(proofreaies and correct through false discovery rate multiple comparisons) to lead be leading of significantly activating of task.

The beacon that all tasks are significantly activated is designated as 1, and the beacon that significantly activates is not designated as 0, and utilizes the matrix of 3*3: [0.25,0.5,0.25; 0.5,1,0.5; 0.25,0.5,0.25] and for smoothing kernel carries out space smoothing, generate the space template m of sensorimotor system.

Location tasks is activated the most significant the 24th to be led and is defined as seed points.Extract the 24th and lead time series as the seed points time series, utilize the coherence analysis method to calculate itself and other 51 synchronicity spectrum distribution of leading under the quiescent condition.Each synchronicity spectral distribution curve of leading as column vector, is arranged in according to the number order of leading and forms space-frequency connection matrix C together.Extract the delegation of space-frequency connection matrix, according to locus reconstruct, namely obtain the synchronicity spatial distribution map of this frequency.Fig. 4 has shown the synchronicity spatial distribution map on 4 Frequency points.

Utilize space template m and space-frequency connection matrix C, calculate each spatial weighting coefficient of leading on each frequency according to formula (7), by the spatial weighting of frequency computation part synchronicity spatial distribution map (be space-frequency connection matrix row vector) and, systematic function connects frequency curve.At last function is connected frequency curve and carry out significance judgement (p＜0.05) by frequency.Shown in Fig. 5 a, function connection frequency range is 0.01-0.0732Hz under the quiescent condition of sensorimotor system.

In order to contrast, also presented the result that traditional arbitration functions connects the methods analyst of frequency range here.Based on identical tranquillization attitude brain function near infrared spectrum imaging data, all time serieses of leading of average left and right sides sensorimotor system realm (as shown in Figure 3) are as the time series of sensorimotor area, the left and right sides respectively, and utilize coherence analysis to calculate synchronicity spectrum distribution between two time serieses, with synchronicity significantly high frequency connect frequency range as function.Shown in Fig. 5 b, based on traditional method, the frequency range that function connects has been determined about 0-0.1025Hz and 1.3Hz.

From above with traditional method and the inventive method result's relatively discovery, the present invention suppressed the many places that exist in the traditional method by noise cause synchronously false, for example＜intrasonic of 0.01Hz (may be drifted about relevant with data synchronously, also might be relevant with physiological intrasonic noise), 0.1Hz about low frequency synchronously (may with 0.1Hz about blood pressure and the variation of heart rate volatility, and Mayer ' swave is relevant), and the high frequency about 1.3Hz synchronous (may fluctuate relevant with heart beating).This explanation the present invention both can be automatically accurately discrimination function connect frequency range, can overcome noise of overall importance in the data connects frequency analysis to function influence again.

Based on method of the present invention, carried out a series of simulation experiment below.Wherein separating method part and result partly set forth.

One, method part:

Analog data D is made of three major parts: the spontaneous fluctuation D of brain ^SF, global noise D ^GN, and machine noise D ^IN, i.e. D=D ^SF+ D ^GN+ D ^INLength of time series is 15 minutes, and sample rate is 10Hz.Observing matrix is 10 * 25.

The spontaneous fluctuation D of brain ^SFAccording to generating with drag:

D ^SF＝T ^SF×S ^SF （9）

Wherein, as shown in Fig. 6 a, T ^SFBe sequence basic time of the spontaneous fluctuation of low frequency, it is made of the phase place sinusoidal wave linear superposition at random of 0.01-0.1Hz, and in the 0.01-0.1Hz section, spectrum distribution is obeyed 0.69*1/f.S ^SFThe spatial distribution of representing spontaneous fluctuation.Shown in Fig. 6 b, it comprises two squared region, represents two region of interest, constitutes function system interested jointly.S ^SFBeing the binaryzation matrix, namely is 1 in the region of interest, and region of interest is outward 0.The outer Base computing of symbol " * " expression vector in the formula (6).

Machine noise D ^INBe that 1 random Gaussian white noise generates by signal to noise ratio, each observation is separate between leading.

Global noise D ^GNComprise two classes, i.e. system's physiological noise D ^GN-1With a moving noise D ^GN-2They produce according to characteristics simulation between the Space Time of each self noise respectively.Concrete grammar is as follows:

1. first kind global noise---the physiological noise D of system ^GN-1

System's physiological noise mainly comprises by systemic physiological activity, the blood oxygen activity that causes as vasomotion of heart beating, breathing and some low frequencies etc.May connect the influence that frequency range produces to measuring ability in order to assess these complicated system's physiological noises, we generate D at each locus k according to following modeling ^GN-1:

D _k ^GN-1=T ^GN-1(t+P _k ^GN-1)×S _k ^GN-1 (10)

Wherein, T ^GN-1Sequence basic time of expression system physiological noise is shown in Fig. 6 c.There are some researches show that the frequency range span of system's physiological noise is very wide, wherein extremely low frequency fluctuation (~ 0.04Hz), the variation of arteriotony and heart rate fluctuation (~ 0.1Hz), and breathe relevant fluctuation (~ 0.25Hz) may with the spontaneous fluctuation D of brain ^SFFrequency have overlapping.And the relevant fluctuation of heart beating (~ 1Hz) relative high frequency, substantially not can with D ^SFFrequency overlap is arranged.Therefore, in order to consider the overlapping and not overlapping two kinds of situations of frequency of signal frequency comprehensively, the frequency range that we simulate physiological noise is 0-0.15Hz.This frequency range not only covers the frequency range (0.01-0.1Hz) of the spontaneous fluctuation of brain, and is also wide than it, and not lap of frequency is arranged.T ^GN-1By amplitude, phase place all at random the sinusoidal wave linear superposition of 0-0.15Hz generate.For the spatial distribution of system's physiological noise, existing studies show that, although there is very strong concordance in observation between leading, there is certain difference in different observations between leading on phase place and amplitude.Therefore, we have simulated phase space scattergram and amplitude spatial distribution map.They lead in each observation different values.Particularly, they are by being generated by the white Gaussian noise matrix behind two-dimentional gaussian kernel (full width at half maximum the is 10) space smoothing.P ^GN-1The amplitude of matrix is limited between 0 to 2.S ^GN-1The amplitude of matrix is limited between 1/3 to 1/2.In order to simulate the time space characteristic of the complexity of system's physiological noise in the actual experiment all sidedly as far as possible, the stochastic simulation experiment has repeated 100 times.

2. the second class global noise---a moving noise D ^GN-2

For the dummy head moving noise may connect the influence that frequency range is differentiated to function, we generate a moving noise D according to following modeling ^GN-2:

D ^GN-2=T ^GN-2×S ^GN-2 (11)

Wherein, T ^GN-2Be sequence basic time of a moving noise, shown in Fig. 6 e.It is made of 1-12 the pulse signal that occurs at random on the time.Approach in the frequency of occurrences of this stature moving noise and the actual experiment.S ^GN-2It is a spatial distribution map of moving noise amplitude.Shown in Fig. 6 f, between difference observation is led, moving noise amplitude is slightly different, this is owing to a moving noise in the actual experiment is because the relative displacement at random of auroral poles medicated cap and scalp causes, and the curvature of diverse location auroral poles medicated cap is different, may cause the amplitude of a moving noise different at diverse location thus.This spatial distribution map is by being generated by the white Gaussian noise matrix behind two-dimentional gaussian kernel (full width at half maximum the is 10) space smoothing.The amplitude of matrix is limited between 3 to 10.In order as far as possible comprehensively to simulate the characteristic that a moving noise in the actual experiment produces on the time at random, and different tested difference, a moving noise is by stochastic simulation 100 times.

3. space masterplate

In the present invention, masterplate m in space has very crucial effect to removing the global noise interference.Therefore, the present invention is extremely important to the robustness of m.In order to simulate the estimation difference of the function system size of m, we enlarge the function system scope of simulating respectively on the original basis or dwindle 0%, 10%, and 20%.In addition, in order to simulate among the m estimation difference to the function system position, we move horizontally 0%, 10% with the function system position of simulation respectively on original position, and 20%.For every kind of situation, estimate that all corresponding function connects frequency range, and the error of tolerance m is to result's influence.

In order to estimate that quantitatively the gained function connects the accuracy of the testing result of frequency range, we have calculated specificity (specificity) and the sensitivity (sensitivity) of testing result according to following formula in simulation experiment:

$\mathrm{specificity}=\frac{F{C}_{\mathrm{df}}^{\mathrm{sw}}\∩F{C}_{\mathrm{df}}^{\mathrm{ideal}}}{F{C}_{\mathrm{df}}^{\mathrm{ideal}}}$

(12)

$\mathrm{sensitivity}=\frac{I-F{C}_{\mathrm{df}}^{\mathrm{sw}}\∪{\mathrm{FC}}_{\mathrm{df}}^{\mathrm{ideal}}}{I-{\mathrm{FC}}_{\mathrm{df}}^{\mathrm{ideal}}}$

Wherein

Expression detects the function that obtains and connects frequency range,

Represent that real function connects frequency range, here we are modeled as 0.01-0.1Hz.I represents the frequency band that specificity and sensitivity are calculated, and we are set at 0-0.25Hz here.In addition, consider that the threshold value selection may connect the influence of frequency range testing result to function, we have also adopted experimenter's performance curve (receiver operator characteristic (ROC)) further to analyze specificity and the sensitivity of this method.

Two, part as a result:

1, first kind global noise---the physiological noise D of system ^GN-1

In the simulation experiment of first, data are according to D=D ^SF+ D ^GN-1+ D ^INGenerate, be used for the gauging system physiological noise function is connected the interference that frequency range detects.

In order to compare, also provided the correlated results that utilizes traditional coherent approach measuring ability to connect frequency range with method of the present invention.Traditional coherent approach is directly with the degree of coherence c of the average time in two representative areas (i.e. Mo Ni two region of interest) in the function system between the sequence ^ROI1-ROI2Function as this system connects frequency characteristic.The average c of 100 stochastic simulation experiments ^ROI1-ROI2Curve presents (variance is thin black line among the figure) in Fig. 7 b.Shown in below black horizontal stripe among Fig. 7 b, it is from 0.0001Hz to 0.1510Hz that the average tranquillization function that traditional method obtains connects frequency range, and the bound standard deviation is respectively 0.0008Hz and 0.0032Hz.This scope is not only located but also has been located all to substantially exceed the goldstandard (0.01-0.1Hz, the middle part Dark grey part among the figure) that function that experiment sets connects frequency range at high frequency (the light grey part in the right side among the figure) in intrasonic (among the figure left side light grey part).The function that detection obtains connects frequency range substantially and the frequency range (0-0.15Hz) of system's physiological noise of experimental simulation is coincide.The specificity of testing result is 0.591 ± 0.020(means standard deviation), sensitivity is 1.This explanation utilizes traditional coherent approach, and system's physiological noise can connect the frequency range detection by obvious interference function.

Fig. 7 c is for utilizing the present invention, and seed points is selected the voxel of the middle of left side region of interest, and the space-frequency connection matrix that obtains and low frequency part thereof are amplified demonstration.Based on this connection matrix, and the space masterplate m as showing among Fig. 6 b, obtain function and connect frequency curve c ^SWAs shown in Fig. 7 d, on average from 0.0090Hz to 0.1033Hz, the bound standard deviation is 0.0026Hz and 0.0041Hz to the function connection frequency range that obtains (below black horizontal stripe among the figure).This result is very approaching with the goldstandard (0.01-0.1Hz, the middle part Dark grey part among the figure) that the function that experiment is set is connected frequency range.The specificity of testing result is 0.931 ± 0.026(means standard deviation), be higher than traditional coherent approach far away.The sensitivity of testing result and traditional coherent approach are similar, are 0.999 ± 0.011.This explanation the present invention can successfully suppress the harmful effect of overall physiological noise.In addition, shown in Fig. 7 e, ROC area under curve of the present invention is 0.975, is significantly higher than the ROC area under curve 0.862 of traditional coherent method.This result has further proved the present invention and traditional method ratio, connects the frequency range detection for function and has higher specificity and sensitivity.

2. the second class global noise---a moving noise D ^GN-2

In the simulation experiment of second portion, data are according to D=D ^SF+ D ^GN-2+ D ^INGenerate, be used for the interference of a research moving noise.

The c that traditional coherent approach obtains ^ROI1-ROI2Curve presents (demonstration of index coordinate) in Fig. 8 a.It has obviously covered very wide frequency range, and on average between 0.0002Hz and 0.754Hz, the bound standard deviation is 0.0014Hz and 0.173Hz to the function that obtains connection frequency range (below black horizontal stripe among the figure).The function that this scope is extensively set in experiment far away connects the goldstandard (0.01-0.1Hz, the middle part Dark grey part among the figure) of frequency range.The specificity of testing result is 0.037 ± 0.178(means standard deviation), sensitivity is 1.The frequency spectrum analysis of correct moving noise shows, the frequency range that the function that the traditional coherent method obtains connects and the frequency range basically identical (result is not presented on the figure) of a moving noise.This explanation is when utilizing the traditional coherent method, and a moving noise can connect the frequency range detection by serious interference function.

The c that the present invention obtains ^SWCurve presents in Fig. 8 b.It concentrates on low frequency part, and the function that obtains connects frequency range and is positioned between 0.0077 ± 0.0039Hz and the 0.113 ± 0.032Hz.Specificity is 0.866 ± 0.197, is significantly higher than traditional method.Sensitivity is 1, and is the same with traditional method.In addition, shown in Fig. 8 c, ROC area under curve of the present invention is 0.976, is higher than 0.850 of traditional method.Further proof the present invention with respect to traditional method to overcoming the effectiveness of an overall situation moving noise.

3. space masterplate

The present invention is presented among Fig. 9 the robust analysis result of space masterplate m.Shown in Fig. 9 a and 9b, when to the excessive or too small estimation of function system size, although c ^SWIntensity can slightly reduce, still be consistent with ideal situation substantially but final function connects the frequency range result.Similarly, shown in Fig. 9 c, when the function system position deviation was estimated, it still was consistent with ideal situation substantially that final function connects the frequency range result.Therefore, experimental result shows, the present invention has certain robustness to the space masterplate to the estimation of function system size and locus.

Above-described specific embodiment; purpose of the present invention, technical scheme and beneficial effect are further described; be understood that; the above only is specific embodiments of the invention; be not limited to the present invention; within the spirit and principles in the present invention all, any modification of making, be equal to replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (8)

1. the method for an automatic discrimination function connection precise frequency scope comprises the steps:

S1, obtain the cerebral function imaging data;

S2, the described cerebral function imaging data of obtaining are carried out pretreatment;

S3, generate the space template of brain function interested system according to prior information;

S4, choose the intrasystem a certain zone of described brain function interested as seed points, adopt signal synchronicity frequency spectrum analysis method to calculate the dependency of all voxels and seed points on each Frequency point;

S5, by frequency ground the correlation results of described all voxels that calculate and seed points is carried out the spatial weighting summation, the function that obtains described brain function interested system connects frequency curve;

S6, function among the step S5 is connected frequency curve carry out significance by frequency and judge.

2. method according to claim 1, it is characterized in that: calculate in the described step 4 that the dependency of all voxels and seed points comprises on each Frequency point: the synchronicity spectral distribution curve that calculates all voxels and seed points in the described zone, form space-frequency connection matrix C, generate the synchronicity spatial distribution map by frequency.

3. method according to claim 1 and 2, it is characterized in that: the signal synchronicity frequency spectrum analysis method in the described step 4 is coherence coherence analysis method.

4. method according to claim 1 and 2, it is characterized in that: the signal synchronicity frequency spectrum analysis method in the described step 4 is the method for Pearson correlation coefficient spectral factorization.

5. method according to claim 3, it is characterized in that: described coherence coherence analysis method is specially:

$c_{\mathrm{xs}}\left(\mathrm{\λ}\right)=\frac{{\left|G_{\mathrm{xs}}\left(\mathrm{\λ}\right)\right|}^2}{G_{\mathrm{xx}}\left(\mathrm{\λ}\right)G_{\mathrm{ss}}\left(\mathrm{\λ}\right)}$

C wherein _XsBe on Frequency point λ (λ), the coherence factor between time series x and the s, its value is 0 ~ 1,0 to be illustrated between this frequency time series and not have linear relationship, has desirable linear relationship between two time serieses of 1 expression.G _XsCrosspower spectrum between express time sequence x and the s, G _XxAnd G _SsRespectively express time sequence x and s from power:

$G_{\mathrm{xs}}\left(\mathrm{\λ}\right)=F_x\left(\mathrm{\λ}\right)F_s^{*}\left(\mathrm{\λ}\right)$

$G_{\mathrm{xx}}\left(\mathrm{\λ}\right)=F_x\left(\mathrm{\λ}\right)F_x^{*}\left(\mathrm{\λ}\right)$

$G_{\mathrm{ss}}\left(\mathrm{\λ}\right)=F_s\left(\mathrm{\λ}\right)F_s^{*}\left(\mathrm{\λ}\right)$

In the formula, F _xAnd F _sThe frequency spectrum that obtains of express time sequence x and the analysis of s Fourier spectrum respectively, * represents conjugation.

6. method according to claim 4, it is characterized in that: the method for described Pearson correlation coefficient spectral factorization comprises:

Seasonal effect in time series counted is designated as N, and Pearson correlation coefficient is calculated as follows:

$c_{\mathrm{xs}}=\frac{\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(t\right)s\left(t\right)}{\sqrt{\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}{x}^2\left(t\right)\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}{s}^2\left(t\right)}}$

Wherein, according to discrete fourier inverse transformation rule:

$x\left(t\right)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{F}_x\left({\mathrm{\λ}}_i\right)e^{i2\mathrm{\π}{\mathrm{\λ}}_{i}t/N}$

$s\left(t\right)=\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{F}_s\left({\mathrm{\λ}}_i\right)e^{i2\mathrm{\π}{\mathrm{\λ}}_{i}t/N}$

Subsequently, with c _XsDecompose on the frequency domain and obtain:

$c_{\mathrm{xs}}\left(\mathrm{\λ}\right)=\frac{N[\mathrm{Re}\left(F_x\left(\mathrm{\λ}\right)\right)\mathrm{Re}\left(F_s\left(\mathrm{\λ}\right)\right)+\mathrm{Im}\left(F_x\left(\mathrm{\λ}\right)\right)\mathrm{Im}\left(F_s\left(\mathrm{\λ}\right)\right)]}{\sqrt{\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}{x}^2\left(t\right)\underset{t=0}{\overset{N-1}{\mathrm{\Σ}}}{s}^2\left(t\right)}}.$

7. according to claim 1,2,5 or 6 one of any described methods, it is characterized in that: by frequency ground the correlation results of described all voxels that calculate and seed points is carried out the spatial weighting summation among the described step S5 and further comprise:

S51, utilize the space template m of the brain function interested system that generates among the described step S3 and each Frequency point λ _iLast synchronicity spatial distribution map c _s(λ _i), calculate the spatial weighting coefficient w of each voxel x correspondence on this Frequency point _x(λ _i)

$w_x\left({\mathrm{\λ}}_i\right)=\frac{m_x-\stackrel{\‾}{m}}{\sqrt{\underset{x=1}{\overset{M}{\mathrm{\Σ}}}{(m_x-\stackrel{\‾}{m})}^2\underset{x=1}{\overset{M}{\mathrm{\Σ}}}{[c_{\mathrm{xs}}\left({\mathrm{\λ}}_i\right)-{\stackrel{\‾}{c}}_s\left({\mathrm{\λ}}_i\right)]}^2}}$

Wherein, M is the total number of voxel, Be the space average number of m, Be c _sThe space average number;

S52, utilize w _x, calculate synchronicity spatial distribution map c by frequency ground _sSpatial weighting and coefficient, that is:

$c^{\mathrm{SW}}\left({\mathrm{\λ}}_i\right)=\underset{x=1}{\overset{M}{\mathrm{\Σ}}}{w}_x\left({\mathrm{\λ}}_i\right)\×c_{\mathrm{xs}}\left({\mathrm{\λ}}_i\right)$

The c of gained ^SWThe function that is function system interested connects frequency curve.

8. according to claim 2,5 or 6 one of any described methods, it is characterized in that: described step S6 comprises:

S61, seed points seasonal effect in time series power spectrum is fixed, its phase frequency spectrum of change at random obtains seed points time series behind the phase randomization by Fourier inversion, and the time series of other voxels is not done any variation;

S62, calculate seed points seasonal effect in time series space-frequency synchronization matrix behind all voxels and the phase randomization according to the method for described step S4, and corresponding function connects frequency curve;

S63, with the randomisation process repeated several times, utilize functions that all obtain at random to connect value in frequency curves, the null hypothesis that generates frequency domain function bonding strength distributes;

S64, calculate the significance level p of the function connection frequency curve value of function system interested, setting threshold p by Frequency point ground ₀, wherein, 0＜p ₀＜1, at p＜p ₀The time, significant frequency is the main frequency that function connects.

Images