The principle and method of fMRI data analysis ———— transferred from network

Principles and methods of fMRI data analysis and processing

Source: Tidy up the file when turned to, the source has not been found the embarrassed feeling written still good, paste here to save.

In recent years, the plasma oxygen level-dependent magnetic resonance brain function imaging (Blood oxygenation level-dependent functional Magnetic resonance imaging, BOLD-FMRI) technology has been very rapid development, In addition to the scanning hardware, scanning technology progress, more effective in the graphic image and other computer science as the core of the relevant disciplines support: image data post-processing technology becomes the key link in fMRI

I. The nature of functional image data

Functional magnetic resonance data includes anatomical (structural) imaging and functional imaging in two categories. The anatomy is based on a high-resolution T1, T2, and FSPGR three-dimensional imaging method. The processing of functional images is the key to fMRI data processing. Because the activity of the cerebral cortex changes rapidly, it requires fast imaging sequence to record the cortical activity caused by a stimulating task, and it is sensitive to the t2* shortening effect of the product of cerebral blood oxygen metabolism, the EPI (Echo planar Imaging) and FLASH ( Fast Low Angle Shot) and other sequences can meet these two conditions, most of which are now using EPI sequence acquisition fMRI function image.

EPI uses a series of inverse gradients on frequency coding, which generates all the signals of an Mr image with a single excitation, Gre-epi (Gradient echo-echo planar Imaging) technology based on small angle excitation, A series of images (several to dozens of) are obtained within a very short TR time. Each acquired image is composed of a brain volume (Volume), which is required to be designed in the fMRI experimental block (Epoch/block Paradigm), and each block must have an integer multiple of tr time. The actual flow of blood is correspondingly a slow process, after which the signal starts to rise after a small descent period, reaches a peak of 4-8 seconds and then slowly drops, recovering 11-14 seconds. In event-related design (event-related Paradigm), if you do not consider the interaction between the two (secondary) tasks, you need to ensure that the interval time is greater than the response time. However, studies have shown that short stimulation intervals do not have much effect on statistical results. (see Figure 1).

The EPI sequence produces hundreds of to thousands of images in a few minutes of experiment (Session) with very fast acquisition speed, and dozens of different time brain volumes become the time series of EPI images (time-series image). At the expense of the resolution of the image, the typical EPI image acquisition matrix is 64x64, which increases the sampling time and leads to more severe image geometric distortion. In addition, EPI sequence images are sensitive to the influence of external magnetic field environment, and the faint bold signal is accompanied by a large number of interfering components. The more prominent problems are:

1. The effect of head movement during the scanning process. Although a variety of physical methods can be used to limit, but the head of the motion is still difficult to completely eliminate, the side effect is far more than the function and structure like superposition fusion of mismatch. Small head of the movement to activate the change in the position of the voxel, resulting in a real functional signal change, the field strength of 1.5 Tesla, the bold signal itself only 0.5-2.0%, but usually two adjacent voxel signal difference is greater than 10%, the brain edge even up to 70%. And the head motion may be the activation related regular movement, will cause the activation area The complete error, seriously affects the experiment result.

2. Susceptibility artifacts. Due to the high-speed switching of the gradient magnetic field produced by the eddy current of the surface strength of the Mr Equipment, the difference in the magnetic susceptibility of the human head tissue, especially near the paranasal sinuses and other air-containing cavities, resulting in uneven local magnetic field, will make the reconstructed EPI image in the direction of phase coding to produce geometric deformation, so that the functional area is

3. The noise interference caused by the scanning equipment and physiological movement is high frequency noise. Physiological movements, including breathing, heartbeat, etc., especially when these movements are related to the task, have a greater effect on the detection of the bold signal. At the same time, because the bold effect is the blood flow regulation, the activation region signal change rate is limited, the physiological spontaneous activity causes the thermal noise and the high time frequency fluctuation, the scanning hardware instability may produce the low frequency drift.

4. N/2 Artifact or Ghost Shadow (N/2 Ghosts), due to inaccurate acquisition timing and uneven static magnetic field, K space alternating echo presents a certain phase difference, in the opposite frequency read out the gradient alternating Mr Signal odd, even echo epi sequence, When the signal is reconstructed by Fourier transform, it appears a pair of false shadows along the phase coding direction. Is the biggest cause of EPI image quality damage.

The large number of EPI image data, low resolution and serious interference factors determine (1) must be removed with the bold signal-related interference signal, improve the signal-to-noise ratio. (2) A large number of time series of four-dimensional EPI images, to be transformed into three-dimensional form of expression. (3) The real bold signal is extracted by the appropriate algorithm. (4) Low-resolution functions like to be fused with high-res anatomical overlays, or to be represented in a known spatial anatomical structure. is the main task of fMRI data processing and analysis. It can be divided into data processing, analysis and presentation of results (see Figure 2).

Second, the function of image data processing

1. Correction (re-alignment).
The correction of head motion is an ideal single mode (Subject) single-mode (modality) registration, which is based on the rigid body motion model, which iteratively calculates the translation and rotation parameters to minimize the mismatch between the reference image (usually the first of the time series) and the subsequent sequence image. Achieve registration of all time series images. Three-dimensional space correction selects three-direction translation and three-axis rotation 6 parameters to describe the head rigid body model; three-dimensional timing also needs to take into account the effect of head motion in each brain volume (tr time), and calibrate each image with two-dimensional registration method respectively. Friston emphasizes the importance of the auto-regressive moving average model (autoregistration moving Average, ARMA) to eliminate the motion effects in the history of object spin excitation.

It is also important to note the slight difference in the acquisition time of each layer of the brain volume in the EPI multilayer acquisition process (dozens of milliseconds). During the block design experiment, because each task block has a long duration (seconds to dozens of seconds), these time differences may not be considered, but when the event-related design, the task excitation time requirements are high, it is necessary to each layer acquisition of different time differences to correct, Ensure that the dozens of-layer image that makes up each brain volume is completed in the same time period. The Sinc method is often used to interpolate the value.

Typically, hundreds of to thousands of images are collected per experiment, and a large amount of data makes the calibration process time-consuming, with some machines with commercial software that take this step to improve processing speed and achieve real-time results. The development of fast motion correction algorithms is very meaningful for real-time imaging (real-time imaging).

2. Registration (registration). The low-resolution EPI feature image often needs to be superimposed on high-resolution anatomical images for functional area identification, which is achieved through registration function activation mapping and anatomical images. Because of the effect of ghost and magnetic susceptibility, the deformation of EPI image is corrected by deconvolution (unwrapping). This is a single multi-modal registration. J. Asbnrner, such as the combined head size and morphology of the Bayesian maximum posterior estimator (Bayesian Maximum a posterior estimator, MAP) method, the use of intermediate images to achieve a variety of functional data and anatomical data accurate registration. But if the spatial normalization of the data is made, these transformations can also be resolved.

3. Normalization (Normalize). Accurate mapping of detected functional areas to high-resolution anatomical charts is the key to fMRI visualization, where functional activation maps do not contain any anatomical information and are not subject to anatomical mapping, but function maps and functional images can share the same coordinate system, so the functional images and anatomical images can be pre-registration, The resulting transformation is applied between the functional map and the anatomical image. The convolution parameters of the average image (Mean image) or the well-calibrated anatomical image of the spatial correction and the pre-designed standard anatomical space (template image) are applied to each tomographic image (Slice image). This ensures that the image data of different samples and different modes are evaluated in the same coordinate system. For single-sample analysis, it can not be normalized to the standard space, but to the template created on its own, for the brain occupying or infarction, such as the obvious damage of the brain structure of the sample image, must not be normalized to the Normal template, automatic algorithm of linear and non-linear conversion process will erase all damaged parts of the unique information, so that the failure For such samples, in addition to the template created by itself, a paid function covering (cost-function masking) technique can be used to cover the lesion site and then normalized to a standard space for comparison. The essence of normalization is a multi-body multi-modal registration.

The Talairach and Tournoux system is the most classic standard anatomical system, with data derived from solid anatomy, and the correspondence between Talairach and Tournoux systems and Brodmann ' s zoning is now well known, with abundant documentation. McGill University Montreal Neurological Institute established MNI system, using 305 normal human Mr Brain scan, mapped to Talairach and Tournoux obtained, such as the famous software SPM99, The standard template uses the MNI system. MNI system does not have the corresponding information with Brodmann ' s partition, MNI system is slightly larger than Talairach and Tournoux system, although some users use the same, but it is best to take a certain way to coordinate point interchange.

Because of the change of global cerebral blood flow and the instability of scanning hardware, the average signal intensity of the average image of time series images changes with the function activity independently, so that the response of each stimulus is not the same level, and the statistical detection function activates the information effect. Each image needs to be adjusted so that its average is equal to the global average, i.e. the normalization of the time series.

4. Smoothing (Smooth).

For the disturbance signal produced by the hardware instability and physiological motion, it can be smoothed out: the spatial smoothing reduces the random noise of Mr Image, improves the signal-to-noise ratio and the detection ability of the function activation data. By integrating the fMRI data with a three-dimensional Gaussian function convolution integral to form a filter, the smoothing range of the filter can be represented by the full width half height (FWHM) of the Gaussian nucleus (Gaussian kernel). Theoretically the Gaussian nucleus should be the same as the scale of the reaction zone, but to ensure that the Gaussian nucleus must be greater than a voxel scale, otherwise it will result in data re-sampling, so that the intrinsic resolution decreases. When the signal-to-noise ratio is low, a wide filter is used, and the active area detected is covered with a larger range. In the case of multi-sample comparisons, the FWHM is also larger (8mm), so that each sample data can be projected onto a common functional anatomy to reduce differences between samples. Although the filter can effectively filter out the noise of a particular frequency, it also sacrifices a real bold signal that is quite a bit of frequency.

For the low frequency drift of time series signals, a filter (fwhm=2.8mm) similar to the bold signal waveform can be used to smooth the time series of each voxel, and if a short tr is used to capture the function image, the frequency band suppression or the minimum mean square Adaptive filter is applied to remove the physiological noises associated with the respiratory heartbeat. Improve the signal-to-noise ratio of the reaction time process, increase the ability of statistical detection signals.

In addition, although the true bold signal is mainly derived from the activation of brain tissue capillary blood oxygen metabolism contribution, but due to the flow of large vessels in the air inflow effect, in the inactive area also has a large number of deoxidation hemoglobin inflow, resulting in increased signal, called "inflow of Pseudo-shadow", appears in more drainage vein cortex area. The pseudo-signal of the low-field-strength machine is more serious, the increase of field strength can reduce the large vessel effect, the SE-EPI sequence can also reduce the inflow effect, and for single-layer EPI imaging, the blood flow sensitivity can be limited by increasing the action time of the RF flip pulse. However, the multi-layered EPI is unable to accommodate a long enough rollover pulse time per layer. Some scholars have reduced their effects by weighting the statistical parameter graph T contrast of various organizations.

Third, the analysis of functional data

After preprocessing the data, we use the appropriate algorithm to extract the real representative pixels, that is, the analysis of the function data. In the first fMRI study, a simple method of image subtraction was used to demonstrate task-dependent brain regions. This is based on the pure insertion hypothesis of the psychology of cognitive subtraction (cognition subtraction) principle, with the task state of the image minus the control state of the image, the difference image high gray value response is the task generated effective active area. This method is too hasty for judging the threshold value of active and inactive voxel, and is particularly sensitive to motion-related effects and other unknown cause disturbances. A more reliable brain graph is a method of using parametric and nonparametric tests.
Commonly used there are 0 hypothesis T test, based on each voxel calculation, weighted average signal difference, t value is greater than the set threshold (such as p=0.05) the voxel is considered to be active, often in the form of pseudo-color expression. Correlation coefficient method in brain graph, each voxel is related to the linear cross correlation coefficient R value, which represents the correlation between the voxel signal strength and the reference function of the time series, and the relationship between the gray scale of the voxel and the desired oxygen metabolic reaction in the process of time series, and the voxel with the correlation greater than the threshold value is considered to be active. There are also F-tests, Z-tests, etc. From a statistical standpoint, these parameter tests can be considered to be a special case of generalized linear models (general Linear model, GLM). GLM is a standard statistical tool developed by K.J Friston and its colleagues for pet data processing, and can be used for analysis of fMRI data by including all interested and non-interested factors in the design matrix, if the time-space autocorrelation between time series can be fully considered.

Hypothesis test, first constructs the statistical parameter mapping of a statistic, calculates the linear relationship between the time process of each voxel reaction and the reference function, determines a threshold value according to the significance level of the test, tests the 0 hypothesis, and identifies the activation and non-activation by the threshold statistical parameter mapping. When constructing the statistical mapping parameter graph, it is important that the gray level time process of each voxel is similar to that of the expected blood flow function. It is very important to set up the corresponding model of hemodynamics accurately by using the corresponding pulse function of the hemodynamic and the convolution integral of an ideal on-off function as the reference function. Several methods have been proposed to establish the corresponding model of hemodynamics. such as Bandettini Fourier spectrum analysis Technology, Bullmore positive and cosine wave linear combination fitting experimental data.

The main problem of the above method is how to select the threshold separating statistic parameter mapping to determine the activation and the inactive voxel. It is very difficult to determine the appropriate uncorrected single-pixel significant threshold value, which can increase the sensitivity of activation detection, but increases the likelihood of inactive zone as the active zone and increases the false positive rate of the test results. And the GLM method hypothesis tests each voxel into a multi-hypothesis (multiple comparisons) test, and the overall voxel test results in more error rates. In order to control false positive, the common Boferroni method is corrected, but it is too conservative, which leads to the decrease of detection rate. Gaussian stochastic theory (Gaussian random Field) can be used to guarantee the occurrence of false positives above the level of voxel, and the same Gaussian check image should be used to smooth the data into Gaussian distribution. For such data, K.J. Friston The grading theory of test can be divided into set (set) level, cluster (Cluster) level and voxel (Voxel) level, although the statistical ability is enhanced, but the spatial distribution characteristic is reduced. It may also be considered that the probability of exceeding the threshold of the real active voxel adjacent cluster is also greater, and the statistical parameter mapping method is separated by the combination intensity threshold and the cluster size threshold, which can reduce the false positive occurrence without reducing the statistical ability. Monte Carlo simulation technology does not require many assumptions, but it is time-consuming.

The previously mentioned voxel-dependent method is only applicable to the experimental data analysis of task design which has been clearly known by time parameters, and it will not be applied to data analysis of spontaneous physiological activities such as sleep, epileptic discharge and so on. The data for such experiments can be analyzed using multivariate analysis of PCA (Principal Component Analysis, PCA) and independent component analysis (independent Component, ICA). The fMRI data is decomposed into orthogonal spatial components or independent components with different time processes, extracting the function information contained in the time series image, without requiring any time process data of the hemodynamic response and a priori hypothesis of the cortical amplitude, and the experimental design is not dependent on any experimental model (such as block or event correlation). Therefore, the single-variable method of voxel-dependent is called mode-driven (paradigm driven), and the corresponding multivariable analysis is called Data Driven analysis mode. PCA determines the distribution characteristics of the functional systems associated with the reaction by detecting the temporal form of the spatial feature pattern with the beginning of the variation of the experimental conditions, and focuses on describing the distribution of functional systems rather than positioning, and is used to explore the interlinkages between the functional areas. ICA, by extracting a series of spatially independent spatial models, focuses more on spatial positioning than PCA, and is best suited to explore the occurrence of a new hypothesis model rather than the test of a known hypothesis. As fMRI studies the central mechanisms of drug action, sleep, and hunger, there have been recent studies of time-clustering (temporal clustering analysis) for epileptic foci with no EEG combined. The disadvantage of PCA and ICA is that it is difficult to give a physiological explanation for the data correlation of most of the different components.

Four, functional magnetic resonance data visualization method

fMRI data is processed and analyzed in an intuitive form to facilitate the observation and reference of results. In addition to the anatomical and mapping parameters, the reconstruction of the cerebral cortex can be used to provide anatomic and geometrical characteristics of the cerebral cortex, and to localize the functional area of the response accordingly. The anatomic image of the standard T1 was divided into gray and cerebrospinal fluid components (Segment), and the anatomical reconstruction of the cortex was performed. In the Retina Brain map (retinotopic map) technique of functional MRI, the cerebral sulcus gyrus structure of occipital lobe was unfolded and evaluated on the plane image. There are many reconstruction methods, such as Voxel-based, 2D contour reconstruction methods, and so on, it is important to maintain the anatomical topological structure of the cortex.

Brodmann's division is now widely used in the localization of the brain functional area, due to the specificity of the cortical structure, in addition to the primary motor and sensory cortex region is more constant, other functional areas and anatomical relationship between the variation is common, unless the nerve anatomy and Talairach and Tournoux systems are familiar with, In general, it is difficult to locate the reaction zone under the naked eye. As mentioned earlier, the individual brain map normalized to the standard brain structure, it is convenient to the reaction zone coordinate point by Brodmann's division to confirm, there are professional software automatic processing.

This paper briefly introduces the principle and method of fMRI data processing and analysis. The implementation of these steps is done by software based on different algorithms. There are a variety of professional software, but the methods and procedures are basically the same. The more general international functional imaging software has integrated processing analysis software, such as the SPM (statistical parametric Mapping) series software and MCW Afni (Medical Wellcome) of the Department of Neurological Imaging Science, University of London, UK College of Wisconsin Analysis of Functional neuroimages) software. As well as special function processing software, like browsing, format conversion of MRIcro software

(http://www.psychology.nottingham.ac.uk/staff/cr1/mricro.html), motion-corrected air software (http://bishopw.loni.ucla.edu/AIR3/), etc. Many are open free software that can be downloaded on the relevant website.

The principle and method of fMRI data analysis ———— transferred from network