Stability in magnitude images in terms of the physiological noise parameter
1. Introduction
Functional MRI studies (fMRI) are based on the blood-oxygenation-level-dependent effect (BOLD) that arises in brain areas where neuronal activity takes place (Ogawa et al., 1990, 1992; 1993). BOLD induces changes in the local magnetic susceptibility and these can be measured by Gradient Echo (GE) Echo-Planar-Imaging (EPI). The fMRI signal thus observed consists of a complex value, which is subdivided into a magnitude and a phase term. In most fMRI studies, the phase signal is discarded and only the magnitude changes are used to detect activated brain areas. In recent years, an increased interest in using both signal types has emerged, with the prospective of increasing both statistical power and spatial specificity.
This chapter describes the fMRI signal from a mathematical and biophysical view point, with emphasis on the phase signal. We will show which kind of information that the phase conveys and exemplify how it can be incorporated in the fMRI analysis pipeline. In section 2, we first describe the basics of an fMRI study, from the fundamental MR signal equations to how a whole brain image is generated in the scanner. We also describe the blood-oxygenation-level-dependent (BOLD) effect that follows neuronal activation. It alters the local magnetic susceptibility and the homogeneity of the local magnetic field and hereby alters the GE EPI signal.
Next, we discuss the phase signal in detail and we characterize the phase in comparison with the more familiar magnitude signal. We illustrate how the probability distributions of the magnitude and the phase values change as the signal-to-noise-ratio increases and how physical mechanisms, like paramagnetic and diamagnetic changes in the local field homogeneity, influence the magnitude and the phase signal in different ways. Regarding neuronal activity, we will review different biophysical models of BOLD with emphasis on the phase signal. Early models describe how the phase effect tend to cancel out for microvascular networks, while a detectable effect remains in larger vessels only (Yablonskiy and Haacke, 1994; Kiselev and Posse 1999; Menon, 2002). On the contrary, a more recent model, based on the Sphere of Lorentz effect, shows that a spatially extended, activated brain area produces a detectable change in the phase signal (Zhao et al., 2007; Feng et al., 2009).
An issue with the phase signal, addressed in the fourth section, is that the potentially valuable phase information is buried in strong phase variations that extend across the head, reflecting variations in the magnetic field homogeneity that vary both in space and in time. These noise factors are unrelated to the neuronal activity and affect the phase to a greater extent than the magnitude signal (Petridou et al., 2009; Hagberg et al., 2008). Methods that suppress or remove these unwanted effects represent a key factor fMRI analysis of phase data. We describe the sources of unwanted phase effects, how these can be modelled, and quantified. We also review several available post-processing methods that can be used to suppress or even remove these unwanted phase effects.
In the last section, we discuss emerging methods that utilize the information conveyed by the phase in the analysis of fMRI data. We will highlight a method that can be used to improve the analysis of the conventional magnitude signal. This method is derived by identifying physiological noise sources from the phase data, and by noise regression of the magnitude signal. The BOLD effect in phase images can also be usedwith the prospective of increasing statistical power and spatial specificity. For illustrative purposes, we show some phase and magnitude results obtained at the group-level in an fMRI study of a motor task performed at 3T.
2. Basics of fMRI physiology
Magnetic Resonance Imaging (MRI) is a non-invasive tomographic technique based on the principles of NMR. In virtue of being a tomographic technique it provides the reconstruction of sections of a 3D object. The capability to explore the structure of the human body non-invasively makes MRI widely employed in clinical imaging. Since the early 1990’s, MRI is used not only to study brain structure but also to study brain function. Theso-called functional Magnetic Resonance Imaging (fMRI), which is based on the BOLD contrast, discovered by Ogawa et al. (1990, 1992, 1993). The great advantage of fMRI resides in its non–invasiveness and high availability, since standard equipment already existing in hospitals can be used.
Neuronal activity causes a local increase in the available oxygen due to neurovascular coupling, either due to synaptic input from other brain areas, or due to local signal processing (Logothetis, 2008). The neurovascular coupling mechanism consists of several steps, leading from an increase in the local neuronal activity to changes in blood flow. Initially, the blood oxygen concentration decreases, causing the initial dip phenomenon of the BOLD effect. As a compensatory mechanism, vasoactive substances are released by glia cells into the extracellular space. When these substances reach nearby blood vessels, the blood flow increases, and so does the relative amount of oxygenated haemoglobin, as a result of changes in the diameter of the blood vessels, the density of the red blood cells, and the amount of O2 and CO2. The hemodynamic BOLD response thus apparently over-compensates the oxygen use in the activated neurons, giving rise to blood oxygen levels which are higher than necessary to re-establish local tissue oxygen availability. The onset of this BOLD response is typically delayed by 1-2 seconds and peaks 4-6 seconds after onset of the neural response. Spatially, the fMRI signal spreads out due to bigger veins that drain the capillaries close to the neurons at work.
2.1. The fundamental MRI signal equations
Any atom holding a net nuclear magnetic moment,
When hydrogen atoms are dipped into an external field,
where γ is the gyromagnetic ratio (in units of radian/T/s), and
At equilibrium, no MR signal can be detected since the magnetization vector
The induced transition to the upper energy level gives rise to a spin coherence, bringing the longitudinal equilibrium magnetization
The transition from the longitudinal to the transverse direction is thus induced by a radio-frequency field
The evolution of the magnetization
In the absence of local magnetic field inhomogeneities, the decay of the transverse component, Mxy, is described by the characteristic time constant, T2, also known as the spin-spin relaxation time, while the characteristic time constant T1, known as spin-lattice relaxation time, rules the recovery of the transverse component, Mz, according to:
In the presence of local field homogeneities, that are not refocused by the imaging sequence the FID signaldecay is governed by the time constant
This situation happens for gradient echo sequences used in fMRI studies. With this MR sequence, local field inhomogeneities are not refocused and in presence of a positive BOLD effect, a lengthening of
A diamagnetic or a paramagnetic shift in magnetic susceptibility causes variations in the local Larmor frequency and may lead to off-resonance and dephasing effects,
Implicitly, this equation takes into account the complex value of the MRI signal, which is described in more detail below. The decay term of [Eq. 5] is obtained from the magnitude component of the MR signal while the off-resonance dephasing effects,
2.2. Image generation and fMRI
The MRI signal is derived from the entire system of spins, which are excited by the RF pulse and precess at the Larmor frequency. In order to spatially distinguish between spins and to reconstruct an image, magnetic field gradients are added to the external magnetic field
The advantages of using gradient-echo EPI areitsrapid image generation with minimal energy deposition, sufficient spatial resolution and ease of acquisition, yielding whole brain coverage within a short imaging time. These advantages come however at the expense of a relative low signal-to-noise ratio, which is however combined with a great BOLD contrast-to-noise ratio.
To reach a sufficient BOLD contrast-to-noise ratio (CNR), several such GE-EPI images are acquired during presentation of selected stimuli, generating an fMRI time series.A typical stimulus paradigm for an fMRI experiment alternates “ON” states in which subjects are asked to perform a task, and “OFF” states (or “Control” state) in which they have to be at rest or perform a control task. Generally “ON” and “OFF” states are repeated every 20-40s, to maximize the BOLD signal, and several repetitions of such “ON”/”OFF” states are required. Typical durations of an fMRI time series is 5-10 min. An example of an fMRI experiment is shown in Fig. 4.
Besides the acquisition sequence and the stimulus paradigm, the success of an fMRI experiment depends uponthe statistical data analysis that is undertaken to individuate the cerebral regions where the variation of the magnitude MRI signal is stimulus correlated. The sensitivity of the statistical analysis depends on both the BOLD contrast and on the fMRI noise level. As we will show, the fMRI noise level can be reduced by using the information conveyed by the phase signal.
3. Characteristics of the MRI phase signal
As we will see, the MRI signal can be described by a complex number, consisting of a magnitude and a phase term. Most MRI methods only retain the magnitude term, but in some established applications, like phase-contrast angiography, Susceptibility Weighted Imaging or anatomical phase imaging methods at high field strengths, the phase is used as a sole ingredient or in combination with the magnitude to generate image contrast.
The phase and the magnitude signals are like two sides of a coin for the MRI signal. While the magnitude signal is related to the dispersion of the underlying local magnetic field, the phase is sensitive to coherent magnetic field shifts. The reason behind this may become clear if we consider the MRI signal as a vector sum of single spin contributions located at different spatial points inside the voxel, precessing around the B0 field according to the fundamental Larmor equation (e.g. Fig.1). In case of an ideal homogeneous magnetic field, without field dispersion or field shifts, all the vectors will precess ‘on-resonance’ at the same frequency generating a maximal magnitude signal at the echo time TE, and a zero phase. In the presence of a field inhomogeneity across the voxel, a field dispersion, each spin senses a slightly different field and precesses at a frequency which is slightly different from the ‘on-resonance’ condition. At TE, the vectors will therefore point into different directions and the vector sum generates a decreased magnitude signal, while the average phase is zero. In presence of a local deviation of the magnetic field which is uniform across the voxel, a field shift, all the vectors have the same frequency albeit different from the ‘on-resonance’ condition. At TE, there is a decreased magnitude signal, and an altered phase.
In the following we will discuss the different characteristics of the magnitude and phase of the MRI signal, with emphasis on some factors that are important for fMRI applications. Besides the signal and noise properties, we will describe three biophysical models for the effect of BOLD on the phase signal.
3.1. Definition of the magnitude and the phase of the fMRI signal
The MR detection systemis based on quadrature detection. The voltage induced by the oscillating net transversemagnetization
In the ideal case, the imaginary component of this complex vector is zero and all the necessary image information is contained in the real component. If this is not the case, it should in principle be possible to recover the real-valued component by ‘phase correcting’ the data. Thereby the complex vector would be rotated into the
In general, the phase evolves linearly with the local magnetic field
where
The phase evolution is defined as positive (clockwise) for increasing magnetic fields (paramagnetic shifts), and negative (counterclockwise) for decreasing magnetic fields (diamagnetic shifts). However this sign convention based on the physical torque effect on the magnetic vector is not compatible with the mathematical sign convention of complex numbers (see Fig. 2 and Hagberg, et al., 2010) and therefore often a negative sign is used in front of the expression in [Eq. 8].
3.2. Signal-to-noise ratio in magnitude and phase images
As we discussed above, the signal in the MRI detection system is acquired via two channels that correspond to the real and the imaginary components of a complex number. The electronic noise of each receiving channel has a Gaussian probability density function (PDF), with zero mean and variance determined by fluctuations around this mean. The PDF of the phase and magnitude noise are determined by the combined noise from both channels (Gudbjartsson and Patz 1996). The magnitude values are always greater than zero, and hence cannot have a zero mean. In presence of purely electronic noise, the noise in the magnitude image therefore follows the Rayleigh probability density function, while for SNR levels of 3 and above a Gaussian PDF is observed. The phase signal originating from pure electronic noise takes on any value in the ±π range with equal probability according to a uniform PDF (Eq. 9). As the SNR increases, the PDF of the phase signals evolve towards a Gaussian shape (Fig. 6). Since these the magnitude and the phase signals have a common origin in the real and imaginary components of the receiving channels, the SNR values are interrelated. For high SNR values, the PDF of the phase is obtained from the magnitude signal (Eq.9).
These considerations evidence how the noise properties of magnitude and phase signals are interrelated at high SNR levels. Although the theory by Gudbjartsson and Patz was developed for an MRI image acquired at a single time point, their observations are also fundamental for describing the magnitude and phase fluctuations that occur during an fMRI time series, described in section 4 (see also Hagberg et al., 2008). Keeping these fluctuations as small as possible is important in fMRI, since the statistical T-values depend on them. Less fluctuations lead to higher BOLD contrast-to-noise ratios and greater statistical T-values, yielding an enhancement of the fMRI detection sensitivity.
3.3. The blood-oxygenation level dependent effect and the phase signal
The contrast in a GE-EPIsequence is generated by the signal evolution between a 90º RF excitation pulse and signal read-out (that is the echo time, TE). While the magnitude signal follows an exponential decay, the phase is linearly accrued during TE (see Eqs. 5, 8). The BOLD effect causesachange in the magnetic susceptibility difference between veins and the surrounding tissue. During the activated state, the oxygenation fraction Y is greater than during baseline states, and therefore the susceptibility difference is diminished. The effect of BOLD on the magnitude signal isdriven by a change in the dispersion (variation) of the magnetic field inside the voxel, while a change in phase is driven by a coherent shift (a net change) of the magnetic field.
The magnitude/phase evolution during TE is influenced not only by the BOLD effect per se but also by the vessel geometry and the spatial orientations of the vessels with respect to the orientation of the static magnetic field of the scanner (see Eq. 10, Fig. 7). Ideally, only the BOLD effect arising in the capillary bed close to the site of neuronal activation should be detected, however larger draining veins also contribute to the activation maps.
In order to calculate the phase variation during BOLD, let’s consider the simple case of a single vessel crossinga voxel. The vessel can be modelled as an infinitely long cylinder with a different, blood-oxygenation and hematocrit dependent susceptibility than the surrounding brain tissue. In this case, the signal depends on the orientation of the vessel with respect to the static magnetic field of the scanner. Vessels oriented parallel to the field, do not ‘oppose’ the field as much as perpendicularly oriented vessels and hence produce a magnetic field shift that is localized inside the vessel, while in the perpendicular orientation, the vesselsproduce both extra- and intra-vascular effects. These are described by the following equations (Boxerman et al., 1995):
where Δχ0 is the susceptibility difference between the extra (
In Fig. 7, phase images of a cylindrical vial containing a solution with higher magnetic susceptibility than the solution contained in the surrounding, outer area are shown for a parallel and a perpendicular orientation, with respect to the static magnetic field.The greatest extra-vascular effects are obtained for the perpendicular orientation, with a spatial cos2Ψ modulation around the vessel, and a diamagnetic phase shift inside the vessel. The size of the phase shifts will depend on the oxygenation fraction and on the geometrical factor, determined by the orientation of the vessels, as expressed by [Eq.10]. The expected phase values that occur during the resting and active conditions of an fMRI study can be obtained by straightforward calculations based on [Eq. 10]. For instance, during rest (when the oxygenation fraction is ca
During an fMRI study the expected phase effect will be diminished due to the partial volume effects that arise when the vessel diameter is smaller than the measured voxel size (see Brainovich et al., 2010). In view of these phase shifts, Menon (2002) suggested that it should be possible to identify voxels dominated by draining vessels parallel to the field, due to the intravascular phase shift, while the contribution from perpendicular vessels should cancel out, due to signal modulation in the extravascular space that is averaged out as a consequence of partial volume effects.
The cortical region that encloses the activated neurons is characterized by a dense capillary network. Larger vessels that drain the activated region are more remotely located at the brain surface (e.g. Keller et al., 2011). The network consists of vessels with a range of different radii, and diffusion effects must also be considered to quantify the effect of the magnetic field distribution. This is not trivial to do and therefore calculations were initially done by Monte Carlo simulations (Yablonskiy and Haacke 1994) before an analytical expression was derived by Kiselev and Posse (1999, see Eqs. 5 and 25 in that paper). From these studies, it emerges that the strength of the BOLD- induced magnitude signal depends on the spatial vessel distribution, consistent with the fact that field dispersion dominates changes in the magnitude signal. On the other hand, if the single vessel model is used for the phase signal, positive and negative phase signals mutually cancel out for random vessel orientations, such as those found in the capillary network, while the phase signals add together for ordered orientations, as was recently explored in a simulation study (Chen and Calhoun, 2010).
These two models of the BOLD phase effect suggest that the interest of analyzing the phase signal resides in the possibility to exclude voxels dominated by draining veins, since predictions based on the summation of magnetic field changes in single vessels cancel out (Hoogenrad et al., 2001). A slightly different picture with regard to the BOLD effect and the variation of the phase signal emerges if the concept of the ‘sphere ofLorentz effect’ is taken into account (Durrant et al., 2003). Zhao et al., 2007 followed by Feng et al., 2009 showed by theory and by experiments that an extended activated brain area actually gives rise to a coherent phase shift. The sphere of Lorentz concept predicts the existence of two main effects of the BOLD response: local field changes caused by the red blood cells, and a bulk magnetic susceptibility shift which is further subdivided into a demagnetizing effect and a volume-averaged magnetization change (Durrant et al., 2003). The volume-averaged effect depends on the blood volume, the vessel geometry and the geometry of their ensuing magnetic field, as well as the size of the susceptibility effect in vessels located in the tissue surrounding the volume (Zhao et al., 2007). The spatial distribution of the phase effect that emerges from this model corresponds to a magnetic dipole and, similarly to the single vessel case, it depends on the geometry and orientation of the activated brain region with respect to the external field (Fig. 8).
Simulations of the BOLD effect that includes the sphere of Lorentz effect can be performed by a three-dimensional Fourier transform of the magnetic field distribution (Marques and Bowtell, 2005). Similar to the approach used by Feng et al. (2009), we performed such calculations for the BOLD response in a cortical region with an elliptic Gaussian field distribution described by:
where
With regard to the time evolution of the phase effect, studies have shown that it closely follows that of the magnitude signal (Zhao et al., 2007). This observation is in agreement with a common BOLD origin for both the phase and magnitude signals.
As we have seen, three differentbiophysical models for describing the phase signal changes caused by the BOLD effect are available in the literature. The first modelis based onthe signal behaviour in large draining veins, the second in capillary networks, while the third model describes the BOLD phase effect caused by a spatially extendedcortical region. The spatial pattern of the BOLD effect is different for the phase than for the magnitude images, but from the literature we know that their temporal evolution is similar (Zhao et al., 2007; Feng et al., 2009).
4. Unwanted phase effects and their removal
The phase value can be considered a fingerprint of the voxel averaged magnetic field. Therefore any effect that causes local or global dishomogeneities of the field leads to phase variations in space. Furthermore, during breathing and motion, the phase will vary from scan to scan in an EPI time series (Hagberg et al., 2008; Petridou et al., 2009). Since the expected BOLD-related phase changes are weak, correcting all unwanted phase fluctuations is of utmost importance.
MRI phase images are characterized by phase wrapping effects dueto an inherent limit of the MRI detection system: the range of phases that can be detected is limited to ±π and therefore phase wraps occur if this limit is exceeded. The true underlying phase value is thus equal to the observed phase ± an integer number of revolutions on the unit circle:
where k can take on any integer value. Such phase wraps occurin space and in time (Hagberg et al., 2008). For instance, close to air-tissue borders where the magnetic susceptibility difference is important, and therefore the phase value also varies greatly, especially if long TE times are used for the measurements, leading to phase jumps in the vicinity of such air-tissue borders (see [Eq.8], and Fig. 5). Temporal phase wraps may also occur in fMRI time series, as the magnetic field may vary with time due to physiological and scanner fluctuations.
4.1. Physiological and scanner noise in phase fMRI images
Unwanted field non-uniformities can be subdivided into static and dynamic magnetic susceptibility differences considering the duration of a 2D image acquisition and they are caused both by the subject and by the scanner (see Hagberg et al., 2008). Static effects are caused by factors that are considered to remain unchanged during the acquisition of a single image, while dynamic effects are caused by factors with a varying influence from slice to slice. Examples of static effects are imperfect shimming and slowly evolving field drifts, examples of dynamic effects are breathing, motion and system vibrations.
The magnetic susceptibility varies greatly across the brain, due to air-tissue bordersas well as due to more local susceptibility differences between the vessels,the sub-cortical structures, the grey and the white matter (Collins et al., 2002, see also the strong phase variations in brain regions just above the frontal sinuses in Fig. 5). When the brain is placed in an external magnetic field, these variations in magnetic susceptibility lead to field dishomogeneties. The scanner hardware compensates for some dishomogeneities by a procedure known as ‘shimming’, i.e. we try to make the field more uniform by activating the shim-coils. How good the shim can get depends on the hardware, and on some clinical scanners only linear correction terms are available, while other commercial systems are capable of performing corrections up to the 2nd order, therefore higher order field imperfections remain, despite the shimming procedure. Since the shim-values are generally set at the beginning of the fMRI experiment they do not compensate for field variations during scanning caused by respiration, that modifies the susceptibility difference between the lungs and the brain, see Ray et al., 2000, and by small micromovements of the brain during the respiratory and cardiac cycles (Le and Hu 1996; van Gelderen et al., 2007). These effects can be noticeable also in the fMRI magnitude signal (see Fig. 10), and must therefore be regressed out if more subtle effects are to be studied, like for instance in resting state studies. New hardware developments may in the future make such corrections unnecessary. One possibility is to change the shim-currents on-the-fly during the fMRI experiment, so that the effects caused by respiration are compensated (van Gelderen et al., 2007). Another possibility is to use magnetic field monitoring during image acquisition (Barmet et al., 2009). This information can either be used for image post-processing or for re-calibration of the shim-currents (Vannesjo et al., 2012).
Not only the subject, but also the scanner contributes to static effects, linked with the aforementioned imperfect shimming issue and slow magnetic field drifts. Time-dependent dynamic effects are caused by helium-pump related system vibrations, and concomitant field gradients arisingfrom Maxwell effects, eddy-current fields, and gradient non-linearities (Hagberg et al., 2008). In a typical fMRI time-series, the scan-to-scan variability of most of these factors can be neglected, except for system vibration, magnetic field drift and heating of the gradient system (Foerster et al., 2005).
4.2. Removal of unwanted phase effects
The unwanted phase effects that originate from different sources can be removed by different post-processing approaches that will be briefly described. We exemplify the effects of these methods on some data sets we acquired (Hagberg et al., 2008; 2012).
Phase wraps, defined in [Eq. 12], can be efficiently removed by performing phase unwrapping. Unwrapping is a straightforward process if it is applied across the time dimension, independently for each voxel in the fMRI time series (Hagberg et al., 2008). Alternatively, spatial phase unwrapping can be applied (e.g. Cusack and Papadakis, 2002; Jenkinson, 2003). Finally, the reference phase method has been proposed, that effectively removes phase wraps by referencing subsequent phase images to the first phase image acquired in the fMRI time series (Tomasi and Caparelli, 2007).
Logically, sources of unwanted phase effects must be unambiguously identified and characterized before removal. Recent studies have shown, that the dominating source of phase variations during fMRI time courses are induced by respiration and the cardiac beat (Petridou et al., 2009; Hagberg et al., 2008; 2012). The importance of respiratory and cardiac effects on the phase data becomes clear from Fourier analysis of the fMRI time series data (Fig.10). After spatio-temporal phase unwrapping, the physiological noise components are slightly less pronounced than after temporal unwrapping alone, while the reference phase method seems more efficient for removal of the respiratory component. The magnitude data is also affected by respiration and the heartbeat, and a very popular method for removal of physiological noise from magnitude data is RETROICOR (Glover et al, 2000). It uses the respiratory and cardiac traces to create regressors for noise removal. The benefit of this approach for phase fMRI time series acquired at 7T has been demonstrated (Petridou et al., 2009).
For long it has been known that information regarding respiration and heartbeat can be extracted from the phase value of the central k-space point (Hu et al., 1995; Le and Hu 1996). One possibility is thus toperform noise regression based on the phase of the central k-space point, instead or in combination with RETROICOR (Hagberg et al., 2012).
Other techniques specifically aim to compensate the scan-to-scan off-resonance effects in k-space. Available techniques are the DORK (acronym for dynamic off-resonance in k-space) and TOAST (temporal off-resonance alignment of single-echo time series technique) methods. In DORK, the global field-changes are estimated by the phase evolution of the central k-space point between subsequent TRs and then each read-out line is corrected in k-space assuming a linear phase accrual in time (Pfeuffer et al., 2002). In TOAST, phase rewinding is performed based on a voxel wise estimate of the B0 changes from the smoothed difference between the phase at each single time point and the average phase in the entire fMRI time series (Hahn et al., 2009). In the next paragraph, we will show how these methods compare in terms of eliminating physiological noise.
Finally, an alternative method that we have explored is based on spatial high-pass filtering. In the human brain at 3T, we observed that the unwanted phase effects were characterized by a low spatial frequency (Hagberg et al., 2012). This characteristic is present also for phantom measurements as can be seen in Fig. 7. We found that the phase fluctuations became similar to the ones observed in magnitude data in healthy human subjects,after application of Gaussian filters that removed all components extending 20mm or more (Hagberg et al., 2012). Fourier analysis also shows that the scanner Helium pump artefact at 2Hz and the respiratory artifact are completely removed by spatial filters, while some cardiac pulsation effects remained (Fig. 10, Hagberg et al., 2012).
The efficiency of each technique for removal of phase instabilities can be evaluated and compared with results from the magnitude data. Krueger and Glover (2001) proposed a model that quantifies the available temporal stability of magnitude data,
In accordance with the Gudbjartsson and Patz relationship between the noise in magnitude and phase images (Eq. 9), we can derive three models for the temporal fluctuations in phase images. In the first model, the noise components induce greater fluctuations in the phase than in the magnitude images, in the second model the contributions are similar in both image types and in the last model the noise is less pronounced in phase than in magnitude images, following the relationship (Hagberg et al., 2008):
where
The
5. On the use of phase information in fMRI studies
The phase signal possesses many interesting features that may be used to extract valuable information without the need for additional scanning time. The question that remains is howthe phase information best is used in order to improve currently available standard fMRI methods. At the time being, no clear-cut answer emerges in the literature. Ongoing work will probably give us an answer in the near future. What follows is a brief overview of some possible approaches that are currently being explored. In the last paragraph, we will also show an example of the outcome of an fMRI analysis performed for both the magnitude and the phase signals within the framework of the standard General Linear Model, using different phase post-processing methods.
We have seen how unwanted phase fluctuations due to respiration or heartbeat are well visible in the phase data. The phase can thus be used to identify and characterize the temporal evolution of different noise sources. This information can be used to remove contributions from common noise sources from the magnitude data sets.
We have seen how the BOLD response influences the phase signal as described by different biophysical models. The dipole of the ensuing magnetic field effect is manifest at different length-scales, depending on whether draining veins or microvascular networks drive the response. Indeed, a phase effect arises both at the level ofsingledraining vessels with diameters of 100-500μm and inextended cortical regionswith BOLD activation. Therefore, the phase signal may be used to exclude single voxels dominated by draining veins (Menon, 2002) or to potentiate the statistical fMRI analysis based on magnitude data (Rowe and Logan, 2004).
5.1. Improving the magnitude signal: noise reduction
Physiological noise is a known limiting factor in fMRI studies (Krüger and Glover, 2001). Improving temporal stability is desirable, as it directly influences the statistical T-values. As pointed out above, studies show that the phase evolution at the centre of k-space in fMRI time series closely reflects field variations due to respiration and the cardiac beat (Hu et al, 1995, Le et al., 1996). Therefore the phase signal can be used as a proxy for physiological (and scanner) fluctuations. The magnitude signal could then be de-noised by simply regressing out unwanted fluctuations (Hagberg et al., 2012). This kind of use of the phase data may thus improve the outcome of fMRI analysis based on magnitude data.
The improvement in temporal stability of the magnitude signal was evaluated for different post-processing strategies (Hagberg et al., 2012). The standard RETROICOR method is based on signal regression (Glover, 2000). The onset of respiration and heart beat is obtained from measurements with a respiration belt and a pulse oxymeter synchronized to the acquisition of the fMRI time series. A set of 8 or more regressors are then created and used for signal regression and de-noising of the magnitude signal. In our approach, that we termed phase based nuisance variable regression (NVR, Hagberg et al., 2012), we used the phase evolution in the centre of the k-space as the only noise regressor. We also evaluated the possibility to combine the RETROICOR noise regressors with the regressor based on the phase evolution in central k-space (Table 1 and Fig. 12). We found that the temporal stability of the magnitude fMRI time course improved by 5% after NVR, which was slightly less than RETROICOR (Hagberg et al., 2012). It should be noted that RETROICOR is based on the cardiac and respiratory traces that must be acquired in synchrony with the fMRI time course, while NVR can be applied without this information. In addition, we found that a combination of RETROICOR and phase-based magnitude correction could improve the overall signal stability more than RETROICOR alone.
Gray matter voxels | 0.0102±0.0042 | 0.0083±0.0041 | 0.0093±0.0041 | 0.0079±0.0039 |
White matter voxels | 0.0056±0.0022 | 0.0040±0.0024 | 0.0050±0.0022 | 0.0035±0.0026 |
CSF voxels | 0.0199±0.0076 | 0.0159±0.0066 | 0.0176±0.0065 | 0.0150±0.0058 |
5.2. Statistical analysis methods incorporating the phase effect
Basically, two methods that incorporate phase signals in the fMRI analysis have been suggested. The first method aims at excluding those voxels from the analysis where the phase changes are great. These voxels are most probably related to draining veins (Menon 2002). The second method is the constant phase model of Rowe and Logan (2004) that extends the General Linear Model GLM to complex values. This model was further refined by Rowe (2005) who combined the real and the imaginary terms of the MRI signal into a single matrix expression. The General Linear Model which is used to evaluate magnitude images is:
For complex numbers, GLM can be rewritten as (Rowe, 2005):
where
Examples of using this model indicate the difficulties that may arise when assessing fMRI results (Rowe 2005). Indeed at least three sets of maps showing activated voxels can be identified, dependent on which signals show significant activations: magnitude only (MO), phase only (PO), or activations in both magnitude and phase images (M&P). Considering the dipolar pattern that the BOLD response has (see Fig. 8), both positive and negative phase changes occur. Therefore we actually end up with as much as 5 sets of activation maps: MO; PO_positive; PO_negative, M&P_negative; M&P_positive (Rowe, 2005; Arja et al., 2010). At the time being, no method is available that brings the analysis of these maps one step further, so that a single activation map showing where in the brain activation takes place with a high sensitivity and specificity, is available. For instance, it is not clear what a PO activation which is visible in the phase but not in the magnitude images actually means. As we have seen in the earlier paragraphs, such an effect could be generated by a draining vessel, or by the sphere of Lorentz effect of the activated cortical patch, or perhaps by an unwanted phase effect that was not completely removed prior to analysis. It is desirable that methods that allow us to combine the activation maps are developedin the future. Such an approach could be based on the available biophysical models for the BOLD effect in phase images, perhaps even by performing quantitative susceptibility mapping (Balla et al., 2012), and may in the future bring us closer to the goal of correctly assigning the cerebral networks where the neuronal activity has its origin at the highest specificity and sensitivity.
5.3. On the use of phase images in an fMRI study
In the following, we will show an example of an fMRI study that incorporates both phase and magnitude images. The goal of the analysis was to investigate possible global phase effects according to the models that describe the phase signal changes caused by the BOLD effect. In order to
Five healthy right handed subjects (three females and two males, 24-27 years) volunteered to participate in the fMRI study (approved by the local IRB). They performed index finger tapping in synchrony with a visual stimulus flashing at 1.7 Hz, in a block design paradigm: left hand tapping for 18s, 9s rest; right hand tapping for 18s, 9s rest; repeated 4 times. Complex valued BOLD fMRI time series were collected (Gradient echo EPI, 75 volumes, TR=3s, TE=39 ms, slice thickness: 3 mm and a 50% gap between 24 slices, bandwidth: 3396 Hz/px, voxel size: 3 x 3 x 4 mm3). Phase wraps and unwanted phase effects related to spatial static magnetic field gradients were removed by three different methods: a) spatio-temporal phase unwrapping, b) pixel wise noise regression based on the phase of the central k-space point and c) spatial high-pass filtering. For the spatial filtering method, a homodyne k-space approach was used based on the following three window functions: Hanning window, Gaussian window (with three different window widths: 5, 10 and 20mm) and a spline window, that was adapted to the magnitude k-space data of each image slice and each timepoint in the fMRI time series (Hagberg et al., 2012). Normalization to the MNI brain space was based on the magnitude image and subsequently applied to the phase images, and conventional spatial smoothing of the pre-processed images was applied to both phase and magnitude images, prior to statistical analysis using SPM5. Statistical analysis was first performed at the single-subject level, then the contrast images were included in a second-level group analysis to reveal the common activation pattern. We thus did not specifically use the Rowe and Logan model in [Eq. 16], but could rely on standard statistical analysis within the framework of the General Linear Model, using standard statistical methods.
Activation patterns for magnitude and phase images achieved after post-processing and statistical analysis are shown in Fig. 13. No significant activations were found in the phase data after a) spatio-temporal unwrapping or b) phase-regression. Only after post-processing by c) spatial homodyne high-pass filtering, significant activation patterns were obtained. Voxels with significant magnitude changes were located in the primary motor cortex. Some of these voxels also showed significant positive changes of the phase, however most voxels with a significant phase effect were located outside M1. In general, the parameter estimates for the phase signal were smaller than for the activations revealed in the magnitude data, leading to lower statistical Z and T-values (Table 2). The T-values for the magnitude images reached levels of 15, while the T-values for the phase data depended on which spatial homodyne window function that had been used. Satisfactory values of the statistical T-value (≥ 8) were obtained with the Homodyne approach using a Hamming window and a Gaussian function with FWHM=20mm, which most probably removes all unwanted phase effects at this field strength and for this echo time (Hagberg et al., 2012). The statistical parametric maps (SPMs) achieved with three of the five filters could be thresholded with p<0.001 (uncorrected for multiple comparison). The other two sets of data were thresholded with p<0.1 uncorrected, in order to reveal the significant activations (Fig. 13).
|
||||||
Cont | # voxels in cluster | Tvalue | Zvalue | p-value | MNI coordinates | |
Hamming n=32 | Pos Neg |
Clmax 659 85 |
8.57 5.84 |
5.08 4.16 |
pu<0.001 pu<0.001 |
34, -26, 46 38, -6, 42 |
Gauss FWHM 20mm |
Pos Neg |
Clmax 792 92 |
8.04 5.74 |
4.93 4.11 |
pu<0.001 pu<0.001 |
32, -26, 48 38, -6, 40 |
Gauss FWHM 10mm |
Pos Neg |
Clmax 841 Clmax 238 |
7.80 5.38 |
4.86 3.96 |
pu<0.01 pu<0.01 |
42, -20 22 38, -6, 44 |
Gauss FWHM 5mm |
Pos Neg |
Clmax 1064 Clmax 381 |
6.61 5.02 |
4.46 3.79 |
pu<0.01 pu<0.01 |
34, -26, 44 -40, -30, 52 |
k-spline filter | Pos Neg |
Clmax 606 Clmax 213 |
8.62 5.28 |
5.10 3.90 |
pu<0.001 pu<0.001 |
42, -20, 18 -40, -30, 42 |
|
||||||
Raw data | Pos | Clmax 950 | 15.06 | 6.37 | FWE=0.5 | 52, -14, 54 |
After identification of voxels with significant signal changes in the second-level group analysis, the timeseries were extracted from these voxels for each subject. In Fig. 13, the extracted phase and magnitude timeseries, in brain regions with significant, stimulus correlated signal changes, are shown. In agreement with the work by Zhao et al. (2007), the time evolution of the phase and magnitude signals closely follow each other. The maximal change of the magnitude signal is 1-2% while the maximal phase change is 0.005-0.010 radians. The contrast-to-noise ratio was smaller for the phase than the magnitude data, both for the negative and for the positive phase changes. With respect to the size of the observed phase changes these were in agreement with those found for the simulated data (Figs. 8 and 9). The size of the activated clusters was 1000 voxels or less, corresponding to a region with an extension (FWHM) of 12mm. For this size, the observed phase changes of 0.005-0.01 radian would correspond to a diamagnetic field shift of ca 0.05-0.10 μT or less, which in turns corresponds to a change in the oxygenation fraction of ca 0.04-0.08 in a capillary network.
This preliminary study suggests how the phase information can be incorporated into the fMRI analysis. In a previous study, (Hagberg et al., 2008) we showed that it is possible to identify voxels with large phase changes in draining veins, by strong positive and negative phase changes close to the brain surface, where the largest vessels are found (Menon 2002). In the present study we used the same set of data but with a different post-processing approach, that included phase post-processing by Homodyne high pass filtering followed by normalization to a standard brain space and spatial smoothing of the final images. This approach allowed to remove contributions from strong positive and negative phase changes in draining vessels and enabled the detection of quite subtle phase changes. These results are consistent with a BOLD response in an extended brain area that generates a magnetic dipole effect (Durrant et al., 2003; Zhao et al., 2007).
6. Conclusion
Functional MRI methods based on the evaluation of magnitude image data may benefit from the inclusion of the information conveyed by the phase images. For instance, it is possible to directly improve the quality of the magnitude data by noise regression derived from the phase signal. In order to reveal subtle BOLD related activations in phase images, the many unwanted phase effects must first be accounted for. Spatial high-pass filtering that removes the spatially slowly evolving
Acknowledgments
The authors would like to thank Valentina Brainovich, who helped us with the fMRI study, and a grant from the Italian Ministry of Health (RF05.103).
References
- 1.
Arja S. K. Feng Z. Chen Z. Caprihan A. Kiehl K. A. Adali T. Calhoun V. D. 2010 Changes in fMRI magnitude data and phase data observed in block-design and event-related tasks 49 4 3149 3160 - 2.
Balla DZ, Sanchez-Panchuelo RM, Wharton SJ, Hagberg GE, Scheffler K, Francis ST, Bowtell R. (2014). Functional quantitative susceptibility mapping (fQSM). Neuroimage, doi: 10.1016/j.neuroimage.2014.06.011. [Epub ahead of print] - 3.
Barmet C. De Zanche N. Pruessmann K. P. 2008 Spatiotemporal magnetic field monitoring for MR Magn Reson Med.60 1 187 97 - 4.
Bernstein M. A. Thomasson D. M. Perman W. H. 1989 Improved detectability in low signal-to-noise ratio magnetic resonance images by means of a phase-corrected real reconstruction. Med Phys.16 5 813 817 - 5.
Boxerman J. L. Hamberg L. M. Rosen B. R. Weisskoff R. M. 1995 MR contrast due to intravascular magnetic susceptibility perturbations Magn Reson Med.;34 4 555 566 - 6.
Brainovich V. Sabatini U. Hagberg G. E. 2009 Advantages of using multiple-echo image combination and asymmetric triangular phase masking in magnetic resonance venography at 3 T. Magn Reson Imaging 27(1)23 EOF 37 EOF - 7.
Bretthorst G. L. 2008 Automatic phasing of MR images. Part II: Voxel-wise phase estim ationJournal of Magnetic Resonance191 193 - 8.
Chen Z. Calhoun V. D. Magnitude and Phase Behavior of Multiresolution BOLD Signal 2010 Conc Magn Reson B 37B(3)129 EOF 145 EOF - 9.
Cusack R. Papadikis N. 2002 New robust 3-D phase unwrapping algorithms: application to magnetic field mapping and undistorting echoplanar images. 16 754 764 - 10.
Durrant C. J. Hertzberg M. P. Kuchel P. W. 2003 Magnetic susceptibility: further insights into macroscopic and microscopic fields and the sphere of Lorentz ConcMagn. Reson. 18A(1): 72- 95. - 11.
Feng Z. Caprihan A. Blagoev K. B. Calhoun V. D. 2009 Biophysical modeling of phase changes in BOLD fMRI. 47 2 540 548 - 12.
Foerster B. U. Tomasi D. Caparelli E. C. 2005 Magnetic field shift due to mechanical vibration in functional magnetic resonance imaging. Magn Reson Med 54(5)1261 EOF 7 EOF - 13.
Glover G. H. 1999 Deconvolution of Impulse Response in Event-Related BOLD fMRI. 9 416 429 - 14.
Glover G. H. Li T. Q. Ress D. 2000 Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn Reson Med 44(1),162 EOF 7 EOF - 15.
Gudbjartsson H. Patz S. 1995 The Rician distribution of noisy MRI data. Magn Reson Med 34(6)910 EOF 4 EOF - 16.
Hagberg G. E. Bianciardi M. Brainovich V. Cassarà A. M. Maraviglia B. 2008 The effect of physiological noise in phase functional from blood oxygen level-dependent effects to direct detection of neuronal currents Magn Reson Imaging 26(7),1026 EOF 1040 EOF - 17.
Hagberg G. E. Welch E. B. Greiser A. 2010 The sign convention for phase values on different vendor systems: definition and implications for susceptibility-weighted imaging. Magn Reson Imaging.28 2 297 300 - 18.
Hagberg G. E. Bianciardi M. Brainovich V. Cassarà A. M. Maraviglia B. 2012 Phase stability in fMRI time series: effect of noise regression, off-resonance correction and spatial filtering techniques. 59 4 3748 3761 - 19.
Hahn A. D. Nencka A. S. Rowe D. B. 2009 Improving robustness and reliability of phase-sensitive fMRI analysis using temporal off-resonance alignment of single-echo timeseries (TOAST). 44 3 742 52 - 20.
Hoogenraad F. G. Pouwels P. J. Hofman M. B. Reichenbach J. R. Sprenger M. Haacke E. M. 2001 Quantitative differentiation between BOLD models in fMRI. Magn Reson Med.45 2 233 246 - 21.
Hu X. Le T. H. Parrish T. Erhard P. 1995 Retrospective estimation and correction of physiological fluctuation in functional MRI. Magn Reson Med 34(2)201 EOF 12 EOF - 22.
Jenkinson M. 2003 Fast, automated N-dimensional phase-unwrapping algorithm. Magn Reson Med. 49(1), 193-197. - 23.
Kiselev V. G. Posse S. 1999 Analytical model of susceptibility-induced MR signal dephasing: effect of diffusion in a microvascular network. Magn Reson Med.41 3 499 509 - 24.
Krüger G. Glover G. H. 2001 Physiological Noise in Oxygenation-Sensitive Magnetic Resonance Imaging. Magn Reson Med 46(4)631 EOF 7 EOF - 25.
Le T. H. Hu X. 1996 Retrospective estimation and correction of physiological artifacts in fMRI by direct extraction of physiological activity from MR data. Magn Reson Med 35(3)290 EOF 8 EOF - 26.
Logothetis N. K. 2008 What we can do and what we cannot do with fMRI. 453 7197 869 78 - 27.
Mansfield P. 1977 Multi-planar image formation using NMR spinechoes. 1. Phys. C. 10:L55 L58. - 28.
Marques J. P. Bowtell R. 2005 Application of a Fourier-based method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility Concepts Magn Reson Part B 25B65 78 - 29.
Menon R. 2002 Postacquisition Suppression of Large-Vessel BOLD Signals in High-Resolution fMRI Magn Reson Med 47(1), 1-9. - 30.
Ogawa S. Lee T. M. Nayak A. S. Glynn P. 1990 Oxygenation-sensitive contrast in magnetic resonance image of rodent brain at high magnetic fields Magn Reson Med.14 1 68 78 - 31.
Ogawa S. Tank D. W. Menon R. Ellermann J. M. Kim S. G. Merkle H. Ugurbil K. 1992 Intrinsic signal changes accompanying sensory stimulation: functional brain mapping with magnetic resonance imaging. Proc Natl Acad Sci U S A.89 13 5951 5955 - 32.
Ogawa S. Menon R. S. Tank D. W. Kim S. G. Merkle H. Ellermann J. M. Ugurbil K. 1993 Functional brain mapping by blood oxygenation level-dependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model. Biophys J64 3 803 812 - 33.
Pfeuffer J. 2002 Correction of Physiologically Induced Global Off-Resonance Effects in Dynamic Echo-Planar and Spiral Functional Imaging. Magn Reson Med.10 344 353 - 34.
Petridou N. Schäfer A. Gowland P. Bowtell R. 2009 Phase vs. magnitude information in functional time series: toward understanding the noise. Magn Reson Imaging 27(8)1046 EOF 57 EOF - 35.
Raj D. Paley D. P. Anderson A. W. Kennan R. P. Gore J. C. 2000 A model for susceptibility artefacts from respiration in functional echo-planar magnetic resonance imaging. Phys Med Biol 45(12)3809 EOF 20 EOF - 36.
Rowe D. B. Logan B. R. 2004 A complex way to compute fMRI activation 1078 EOF 1092 EOF - 37.
Rowe D. B. 2005 Modeling both the magnitude and phase of complex-valued fMRI data 25 1310 1324 - 38.
Schmitt F. Stehling M. K. Turner R. 1998 Echo-Planar Imaging: Theory, Technique and Application. Springer Verlag - 39.
Van Gelderen P. De Zwart J. A. Starewicz P. Hinks R. S. Duyn J. H. 2007 Real-time shimmingto compensate for respiration induced Bo fluctuations. Magn Reson Med57 362 368 - 40.
Tomasi D. G. Caparelli E. G. 2007 Macrovascular cobtribution in activation patterns of working memory. J Cereb Blood Flow Metab27 33 42 - 41.
Vannesjo SJ Haeberlin M Kasper L Pavan M Wilm BJ Barmet C Pruessmann KP 2012 . Gradient system characterization by impulse response measurements with a dynamic field camera. Magn Reson Med.doi: 10.1002/mrm.24263. - 42.
Yablonskiy D. A. Haacke E. M. 1994 Theory of NMR signal behaviour in magnetically inhomogeneous tissues: the static dephasing regime. Magn Reson Med32 6 749 763 - 43.
Zhao F. Jin T. Wang P. Hu X. Kim S. G. 2007 Sources of phase changes in BOLD and CBV-weighted fMRI. Magn Reson Med57 520 527