Introduction

Each activity a person performs is managed by a certain location of the brain. The location of a brain directly related to an activity can be visualized using images from functional magnetic resonance imaging (fMRI) instrument. These images are obtained using the changes between active and non-active state of location of a brain. The image contrast obtained this way is very small, and fMRI instrument is so sensitive that it picks unwanted signals or noise, that induce distortions of the actual experimental signals, which can be interpreted as false brain activities. The objective of this paper is to remove noise or distortions, which are unrelated to the experimental fMRI signals. The methods, both periodic discrete and translation invariant packet wavelets are suggested to remove noise from fMRI data due to their decorrelating effects in a temporal domain.

The main sources of noise are not fully understood. A number of possible sources have been suggested, for example, slow phase variations in the MR images due to respiration movements, cardiac and other physiological processes, patient movement, and local changes in the magnetic field due to scanner instabilities. During fMRI image acquisition process, high frequency components like heart rate (0.6 –1.2 Hz) and respiration (0.1 – 0.5 Hz) are under sampled with typical repeat times (TRs ) of 3 to 7s and can, according to Nyquist's theorem, be expressed as low-frequency (0.1 Hz) signal components or aliased higher frequency signals [₁].

Due to the fact that Rician distribution is used in the physics of magnetic resonance and the Gaussian distribution in functional neuroimaging, the noise, in the blood oxygenation dependent (BOLD) response of an fMRI data, is both Gaussian and Rician distributed. For high SNR fMRI data, Rician distributed noise is symmetric, thus, can be considered as Gaussian distributed. For low SNR fMRI data, there is a difference between Gaussian and Rician distributed noise, i.e, an image with low intensity and Rician distributed has probability density which is asymmetric. It was further shown that the difference between two Rician distributed images is symmetric and Gaussian distributed. Since an fMRI BOLD response is the difference between two BOLD responses, the active and passive state, thus, the noise in the resulting BOLD response is symmetrically distributed, which is characteristic of Gaussian distribution.

To remove noise from an fMRI signal, general wavelet method is presented. This method decorrelates and extracts noise from a BOLD response leaving a clean event related BOLD response of the fMRI data. To remove the noise, first, the noise is decorrelated and then the intensity of the noise is estimated statistically. Second, the estimated noise is removed from the fMRI data in a temporal domain. The wavelet methods, both periodic and packet based wavelets, a wiener filter is used with coefficients estimated from the noisy fMRI data adaptively. A wavelet transform approximately acts as a detrender, which facilitates the removal of noise from fMRI data in a temporal domain without affecting an event related BOLD response. Real and simulated fMRI data is used to confirm that the methods remove noise and preserve activities related BOLD response of an fMRI data.

After denoising, to reduce the false discovery rate, clustering is performed in spatial domain based on optimal minimal cluster size obtained from Monte Carlo simulation for a user-defined confidence level. In this work, it is shown that the suggested methods can recover a signal from a data with a noise that has intensity of any standard deviation away from the center of the data but an optimal performance is obtained when the noise level is less than seven standard deviations from the center.

Modeling of fMRI data

Overall, we need to model and remove deterministic components from the time series before proceeding with the statistical analysis. The fMRI signal y_m(t) at point M in the brain is given by:

y_m(t) = β_mx(t) + n(t) (Eq 3.1)

Where n(t) is stationary Gaussian noise, – β_m is a scalar that measures the strength of the response at the voxel M. When we process y_m(t) using wavelet filtering and the resulting equation will be:

Wy_m(t) = Wβ_mx(t) + Wn(t) (Eq 3.2)

where W is coefficient of wavelet. Note : y_m(t) can be extended to accommodate 2-D by deriving the locationM of y_m(t) from the row or column of the 2-D fMRI image and processing each row and column as separate time series and combining the final result. Equation 3-2 makes the noise uniform, i.e, N (0,σ² ) , throughout the fMRI time series. Once the noise is uniform, it is easy to remove it. During the noise removal process, the temporal behavior of the voxel time series is used to remove noise, i.e., the temporal domain decorrelates the noise and the noise level or standard deviation of the noise can be estimated and used in the noise normalization process to determines the coherent signal as described later.

Background

Statistical Parametric Mapping (SPM) tool and many research works use Gaussian filtering to remove nosie from an fMRI data. Gaussian filtering is low-pass filtering, that is being used traditionally to process fMRI data, but it can remove relevant detail information. Gaussian filtering requires a window or a kernel size to be derived from the fMRI data itself to avoid processing errors . The other most widely used methods are Fast Fourier Transform (FFT) and SPline, but FFT unable to identify sharp transient events that are similar to fMRI data,. SPline is FFT and wavelet based, its sharp frequency characterization makes a good fit to process an fMRI or time related signal[₂] but it needs optimization [₃].

Most of the existing literatures suggest performing both the analysis and filtering in the time domain or temporal domain, and after the analysis phase in the time domain, an inverse transform is applied to reconstruct an activation map from the coefficients that are designated as significant. While this reconstructed map is very useful for visualization purposes, it does not have a direct statistical interpretation, that is, the statistical parameters, such as, t or z values are computed in the time domain and there is no straightforward way to map the statistics to the spatial domain.

In this paper both statistical analysis and filtering are performed separately, this approach helps optimize the detection of false positive active points of the voxel time series.

Methodology – Estimation of Noise and its removal from fMRI

For a continuous variate x, the median Ḿ defined, using the (cumulative) distribution function F(x), by the condition:

F(m) = ½ ( 5 1)

This means that one half of the observations are below and the other half above the median. For a sample of ordered variables x₁,x₂, …x_n the sample median, denoted as ḿ = median {x_i} is given by (with integer k ); Note : these ordered variables are fMRI time series high passed using wavelet, and decimated.

Figure 1

MAD = med{|x_i - ḿ|} for i = 1,2,…,n (Eq. 5.3)

For a normal distribution, the probability density function is

Figure 2

As for any symmetrical distribution, mean and median coincide, thus µ = ḿ. In addition, we assumed fMRI mean µ is centered at zero (0); thus MAD = median{|x|}. According to ( 5 1) we then have to evaluate the limits, -α and +α, for which ∫^α _-αp(x) = ½ OR ∫^α ₀p(x) = ¼ which gives the numerical solution α = 0.6745, and thus MAD can be linked, for a normal distribution, with the average standard deviation σ(x) of a single observation x, by MAD = ασ(x) which is MAD/0.6745 = σ(x); Through out this work this relationship is used to estimate noise for 1-D and 2-D fMRI signal at each voxel location. It is shown that the noise estimated this way is almost the same as known noise level, i.e., when adding a known noise level to a signal and filtering using wiener filter as specified in this paper has the same effect as estimating the noise and filtering it the same way.

The DICOM data with activation points is collected form patients and transformed into AFNI [₄] format. The signal from AFNI is processed as specified below and ported to AFNI for visualization and wavelet processed to determine the activation points in the brain and to compare and contrast the processed signal with othe fMRI processing methods. The wavelet filter used is symmlet with size four(4). The procedures employed to separate the coherent signal form the noise are as follows:

1. For processing 1-D fMRI data :Divide the signal into 1-D representations or into 1-D time-series; if the 1-D signal is not power of two, padding with zeros several hundred positions is necessary

2. For Processing 2-D fMRI data: Divide the 3D fMRI data into slices and if the size of a slice is not power of two, padding with zeros several hundred positions is necessary

3. Wavelet decompose the fMRI signal, this step makes the noise distribution N(0, σ2)

4. 1. Decompose the fMRI data into discrete wavelet coefficients (wi)

   2. Decompose the fMRI data into packet wavelet coefficients (wi)

5. Estimation of the noise level σ from its distribution N(0, σ2): σ = MAD(wi) / 0.6745 for both periodic discrete and packet wavelets, separately

6. Select the representative coefficients - wj form the wavelet packet coefficients - wi of step 2b; to select the representative coefficients, Shannon entropy is used . The selection process involves the following steps :

7. 1. Start from the bottom of the decomposition tree and mark every thing

   2. Determine the entropy of each element of the tree, and then if the entropy of a parent node is less than or equal to the sum of entropy of its the two children, them unmark the children and mark the parent, otherwise leave the children marked, continue this way until you reach the top. Shanon entropy is determined using the formula E(X) = -Σnpnlog2pn), where pn = (xn)2/ ||x||2 , where xn is the nth component of x, ||x|| is the signal length and E(X) is the entropy of the fMRI signal x.

8. Normalize the noise in the temporal domain to make the noise distribution N(0,1) : wi = wi/σ for discrete wavelet and for packet wavelet wj = wj / σ, and σ is a result from discrete or packet wavelets used in step 4

9. Determine the coefficient of wiener filtering of the wavelet coefficient: Φ = |wi 2| / (|wi 2| + |σ2|) for discrete wavelet and for packet wavelet Φ = |wj 2| / (|wj 2| + |σ2|)

10. Get the wavelet coefficients from the noisy fMRI time series and or slice from step 3, and then apply wiener filter from step 7 to determine the clean wavelet coefficients: wi clean = wi⋅Φ

11. Reconstruct the clean wavelet coefficient (wi clean ) into an fMRI time series or slice

12. Go to step 3 and repeat the steps until the entire time series or fMRI slices are processed.

Results

The brain processes different information in different way, for example, when the subject was instructed to process items according to their meanings ( is the word hot or cold ?) or when the subject is instructed to perform bilateral finger tapping at a set of interval of time. The former task was associated with activations in a set of brain regions including left lateral prefrontal cortex (PFC) and left medial temporal cortex. The latter showed relatively greater activation in right and left PFC.

The mapping of functions to brain regions presupposes, among other things, a thorough understanding of the cognitive functions that are to be mapped onto the brain regions. This understanding, however, while developing, is still rudimentary. Someone with schizophrenia, when confronted with a psychological task, might tackle it in a very different way, in terms of the cognitive strategies used, and from a healthy person confronted with the same task. The observation that brain activity differs across the two individuals would only be interpretable insofar as one thoroughly understood the processes that each individual invoked in response to the task demands.

Recent studies have found that the spatial extent of fMRI activation in healthy older adults is approximately half that of younger subjects, for both visual and task related activities. Also, it is very important to note this fact: the difference, in activation, between healthy older and young subject, is not associated with reduced hemodynamic response (HDR) amplitude in older subjects, because young and elderly adults have similar HDR amplitudes [₅], and have similar distributions of HDR amplitudes across voxels. Furthermore, greater head motion in the elderly does not cause this difference, although head motion differences can affect spatial extent of activation [₅]. Instead, spatial extent differences are associated with higher voxelwise noise levels in elderly adults, perhaps due to changes in cardiac or respiratory effects upon the fMRI signal .

In this experiment no comparison across subject is performed and an individual subject data is analyzed independently, and the subject was instructed to perform bilateral finger tapping as brain activation mapping of this action can be seen on both sides of lateral prefrontal cortex, and this finding, Figure ,is consistent with theoretical basis of neural functions for a healthy individual, the complete result of the processing an fMRI obtained from a healthy person is shown in Figure.

Figure 3

Figure 1: Cluster of Activation Points when the subject is finger tapping

Figure 4

Figure 2: Cluster of Activation Points when the subject is finger tapping generated using euclidean distance and minimum cluster size.

The other most important thing to notice is, during the experiment, sometimes the well instructed subjects do all kinds of things in addition to the task they were asked to perform and this showed up in the fMRI images, as seen in Figure, and interpretation of the results should take into consideration this kind of events. Furthermore, the best strategy is to mask the data to include only the areas or regions that corresponding to the activities in the brain.

Acknowledgements