Open access

Phase Variations in fMRI Time Series Analysis: Friend or Foe?

Written By

Gisela Hagberg and Elisa Tuzzi

Published: 31 May 2014

DOI: 10.5772/58275

From the Edited Volume

Advanced Brain Neuroimaging Topics in Health and Disease - Methods and Applications

Edited by T. Dorina Papageorgiou, George I. Christopoulos and Stelios M. Smirnakis

Chapter metrics overview

2,473 Chapter Downloads

View Full Metrics

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, I, can generate an MR signal detectable by a scanner. Most MRI methods employ hydrogen atoms 1H, owing to a high nuclear gyromagnetic moment and the high water concentrations present in the human body (ca 55M).

When hydrogen atoms are dipped into an external field, B0, the spins tend to either align with the external field (lower energy state) or orient themselves in the opposite direction (upper energy state). At equilibrium, spins are in light excess in the former state, resulting in a non-zero net magnetic moment Mthat is aligned with the external field (Fig. 1). According to the rules of quantum mechanics, the nuclear magnetic moment cannot be exactly aligned with the external magnetic field, and are therefore oriented at an angle with respect to the field. This difference in orientation between the external field and the field of the spins causes a torque effect leading to a precession of the spins around the external field (see zoomed box, Fig. 1). The precession occurs at the Larmor frequency, ω , according to the fundamental equation of NMR, the Larmor equation:

ω = γ B 0 E1

where γ is the gyromagnetic ratio (in units of radian/T/s), and B0 the intensity of the magnetic field.

Figure 1.

Net magnetic moment Mof and ensemble of spins, each with an individual nuclear magnetic moment, I, placed in an external field B0. In the zoomed inset, the precession of I around the B0 axis at the Larmor frequency, ω is shown.

At equilibrium, no MR signal can be detected since the magnetization vector M is time independent. An MR signal can only be detected after the equilibrium condition has been perturbed, when the spin systems return to equilibrium and the spins flip from the higher to the lower energy level. The perturbation is made so that the spins can absorb the necessary energy to bridge the gap, ΔE, between the two spin states. A radiation field B1 oscillating exactly at the Larmor frequency is used to excite the spins and induce transitions from the lower to the higher energy state. The perturbing B1 field is applied orthogonally with respect to B0 (Fig. 2).

The induced transition to the upper energy level gives rise to a spin coherence, bringing the longitudinal equilibrium magnetization Mz into the transverse plane x-y, Mxy, in a fixed reference system (laboratory system) as shown in Fig. 2.

The transition from the longitudinal to the transverse direction is thus induced by a radio-frequency field B1, or a so-called 90° pulse. When the B1 90° pulse is interrupted,M starts to precess around the B0 z-axis at the Larmor frequency, whilst re-establishing its initial equilibrium value.This causes the gradual recovery of the longitudinal component Mz and the complete extinction of Mxy. The temporal variation of the transverse component, Mxy, originates the MR signal, also known as Free Induction Decay (FID),oscillates at the Larmor frequency and is detectable as an induced voltage in the radio-frequency coil, according to Faraday’s law of electromagnetic induction.

The evolution of the magnetization M of a spin system immersed in an external magnetic field, B0, perturbed by an oscillating field, B1, and the forces that rule the FID can be conjointly described by the phenomenological Bloch equations. The solutions to the Bloch equations govern the decay and recovery 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:

M x y = M 0 e t / T 2 E2
M z = M 0 ( 1 e t / T 1 ) E3

In the presence of local field homogeneities, that are not refocused by the imaging sequence the FID signaldecay is governed by the time constant T2*, or the effective transverse relaxation time:

M x y = M 0 e t / T * 2 E4

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 T2* arises due to a relative increase in the local blood concentration of oxyhemoglobin during activation. Oxyhemoglobin is less paramagnetic than deoxyhemoglobin and therefore generates a field shift which is diamagnetic relative to the baseline condition (see paragraph 4 and Fig 8). As will be discussed in greater detail further on, the presence of paramagnetic effects leads to an increase in the bulk magnetic field, with faster spin precession while the presence of diamagnetic effects causes a decreased local field, and a loss of angular precession velocity of the spins around Bo (see Eq. 1 and Fig. 2)

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, φ(t) depending on whether the field is uniform or varies across the voxel. Therefore a more general signal equation for a gradient echo sequences is:

M x y = M o e t / T * 2 e i φ ( t ) φ ( t ) = γ 0 t B ( τ ) d τ E5

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, φ(t) are obtained from the phase.

Figure 2.

At equilibrium, the net magnetization is said to be longitudinal, Mz and is oriented along the direction of the external magnetic field, B0. A radio frequency (RF) field B1oscillating at the Larmor frequency is applied orthogonally to B0. This so-called 90° pulse, rotates the magnetization vector into the transverse orientation Mxy. During the echo time, spin coherence is lost leading to a T2* signal decay, that depends on the local magnetic susceptibility. Spins in a diamagnetic environment precesses a little slower, while spins in a paramagnetic environment precesses a little faster than the on-resonance spins. The original equilibrium condition with a longitudinal net magnetization is recovered after a time that is about five times as long as the longitudinal T1 relaxation time.

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 B0. These gradients are applied in the three spatial directions (x,y,z), defining respectively the slice-selection gradient, Gz, the phase-encode gradient, Gy, and thefrequency-encode (also readout) gradient, Gx. The first gradient is used to select the magnetization located inside a slice through the object, while the latter two gradients are used to encode the spatial position of this magnetization in 2D space. This information is sampled in a reciprocal space, termed k-space, and the image of the selected slice is reconstructed by a 2D Fourier Transform that maps the original magnetization distribution of the object from k-space back into object space.In an fMRI study, the signal changes caused by the BOLD response is followed by whole brain measurements that are obtained every 1-3 seconds, with typical voxel sizes between 3x3x3 and 1x1x1 mm3 in human brain studies. Due to the sluggishness of the BOLD functional brain studies generally employ gradient-echo echo-planar imaging (GE-EPI) sequences for data acquisition: a 2D image (or a slice) is acquired following a single excitation pulse, and 20-30 such slices are imaged yielding whole brain coverage within a 2-3 s long repetition time, with typical voxel sizes between 1x1x1 and 3x3x3 mm3 (Schmitt et al., 1998). The GE-EPI sequence, in which the entire 2D k-space is filled using rapid gradient switching after a single excitation pulse, begins with a slice-selective excitation pulse with a flip angle α(i.e. an RF pulse that brings the Magnetization vector at an angle αwith respect to B0), and is followed by an EPI readout of gradient echoes (i.e. a train of gradient echoes made up of readout gradientswith alternating polarity, see Gx in Fig.3). Before each readout gradient, a brief phase-encode blip (Gy) is applied in order to quickly sample the whole k-space. Alternatively, a full 3D volume can be acquired by so-called echo-volume imaging (EVI, Mansfield, 1977). In this case no slice-selection is necessary, instead an entire volume is excited by a single RF-pulse and 3D k-space is covered by phase encoding along two spatial directions.

Figure 3.

The gradient echo (GE) echo-planar-imaging (EPI) sequence. An RF pulse with flip angle α (first row) is applied together with a magnetic field gradient in the slice selection direction (Gz. last row). Hereby all the spins located inside the selected slice have theirmagnetization vector rotated by an angle αwith respect to B0. k-space (the reciprocal of the image space) is sampled by a combination of brief phase-encoding blips (Gx second row) and readout gradientswith alternating polarity (Gy third row). Finally, a Fourier Transform is performed to recover the MR image of the selected slice.

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.

Figure 4.

An fMRI experiment consisting of “OFF” and “ON” periods with a 20-40 seconds duration. Gradient-echo EPI images with whole brain coverage are acquired each 2-3 seconds, and 100 or more brain volumes make up the fMRI time series.

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 Mxy is split into two channels that function with a phase difference with respect to each other. Two output signals are thus detected, the first which is ‘in-phase’ and the second which is ‘in-quadrature’ (ie phase shifted) with respect to the first. These signals can be described in mathematical terms by the real and the imaginary components of a complex vector, corresponding to the Mx and the My components of the detected net magnetization of the FID.

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 Mx direction, and an absorption mode image of the magnetization vector would be available, as suggested bysome promising attempts (Bernstein et al., 1989, Bretthorst, 2008). However in general, due to phase imperfections that are difficult to correct, it is convenient to use the absolute (magnitude) value of the MRI signal, corresponding to the length of the complex vector. Such MRI images are produced by taking the square root of the sum of squares of the real and imaginary components of the measured magnetization vector. The phase of the complex vector, or the argument of the complex vector, is obtained by evaluating the arc-tangent of the ratio of the imaginary and the real components (see Fig. 5):

S = a b s ( M x y ) = M x 2 + M y 2 E6
φ = arg ( M x y ) = arctan ( M y M x ) E7

Figure 5.

The MRI detection system subdivides the image information into a real and an imaginary component, corresponding to the Mx and the My components of the transverse magnetization. These components define a complex vector and can be rearranged to form the length [Eq. 6] and the argument [Eq.7] of the vector, corresponding to the MRI magnitude and phase signals, respectively. The red arrow indicates the area just above the frontal sinuses with strong magnetic field variations, due to the difference in susceptibility between the air in the frontal sinuses and the brain tissue.

In general, the phase evolves linearly with the local magnetic field Bloc in a particular voxel during the time between the excitation pulse and signal read out, i.e. the echo time, TE. At the time of acquisition the phase is:

Δ φ = γ B l o c T E E8

where γ is the gyromagnetic ratio (in units of radian/T/s) for the observed nuclear spins (see Eq.2).

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).

Figure 6.

Histograms representing the signal-to-noise probability density functions (PDF) for the magnitude (left column) and phase (right column) signals. The following cases are shown. First row: with an SNR of zero (purely electronic noise), the magnitude values distribute according to a Rayleigh PDF, while the distribution of the phase values is uniform. In row 2-4, we show how the shape of the PDF changes as we go towards greater SNR. At low SNR levels of 1 (2nd row) a Rician PDF is observed for the magnitude data, while at an SNR of 3 (3rd row) and 10 (4th row) the shape is Gaussian. Above an SNR of 3, the phase signal also follows a Gaussian PDF, albeit with a more narrow shape, as the standard deviation of the phase signal is the inverse of the SNR observed for the magnitude signal.

S D φ = { π 3 if SNR = 0 1 S N R if SNR > > 1 E9

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):

Δ B v e s s e l = i v + e v i v = 2 π Δ χ o ( 1 Y ) H c t B o ( cos 2 θ 1 3 ) e v = 2 π Δ χ o ( 1 Y ) H c t B o ( R r ) 2 sin 2 θ cos ( 2 Ψ ) E10

where Δχ0 is the susceptibility difference between the extra (ev) - and intravascular (iv) space (in cgs units), Y the oxygenation level of the blood, Hct is the haematocrit level, B0 is the static magnetic field strength, R the vessel radius and r the distance from the centre of the vessel cylinder, and θ, Ψ are the angles defined in Fig 7.

Figure 7.

Biophysical cylinder model for a single vessel placed with its length axis oriented parallel (A,D, symbol: =) or perpendicular (C, F, symbol: ) to the external magnetic field, B0. B: The angles θ, Ψ and the distance from the centre of the vessel, r, define the position for a point in space where a magnetic field is generated due to the blood-tissues susceptibility difference, as described in [Eq. 10]. A,C: Simulated magnetic field map for a cylinder placed in the centre of a sphere (first row: Δχ0=0.18ppm, hematocrit=0.4, oxygenation fraction Y=0.54). For a cylinder oriented parallel to the field, no phase modulation occur in the extravascular space (the field inside the cylinder is visible as a blue spot in the centre). For a cylinder oriented perpendicular to the external field, a cos2Ψ modulations occurs around the vessel. D,F: Measured MRI phase images (Magnetic field strength: 3T, GE-EPI, TE=20ms) for a cylinder with higher susceptibility than the surrounding water. The pattern corresponding to the dipole field is evident for the perpendicular orientation. Unwanted background variations of the field are present in the measured but not in the simulated data.E: Depending on the vessel orientation (= or ) with respect to the external field, diamagnetic or paramagnetic field shifts can be observed.

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 Y=0.54) the maximal field difference between the extra- and intravascular space amounts to ±0.62μT at 3T and appears in the tissue surrounding a perpendicular vessel (extravascular space). This value corresponds to a phase shift as high as ±4.98 radian if a GE-EPI sequence with a TE of 30 ms is used. In presence of a BOLD response during the activated state (Y=0.68) the susceptibility difference between the blood and the tissue diminishes and therefore the maximal field difference becomes ±0.43μT (±3.45 radian). In the intravascular space of a parallel vessel, the field undergoes a diamagnetic shift and is therefore reduced from 0.42 μT (3.37 radian) at rest to 0.29 μT (2.33 radian) during activation. Consequently, the expected BOLD induced field change is maximal in the extravascular space for perpendicular vessels (±0.19μT, 1.52 radian), while for the parallel case it is greatest in the intravascular space (-0.13μT, 1.04 radian). These observations are also in line with previous work (Hoogenrad et al., 2001).

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).

Figure 8.

Simulation of the spatial pattern of the phase effect generated by the BOLD response in a brain area with an elliptic Gaussian shape (contours of the region are overlaid on the phase images). The external, static, magnetic field, B0 is applied along the Z-direction, and the phase effect is shown in three planes: XZ, YZ, and YX. Similarly to the case of a single vessel, the phase effect depends on the orientation of the cortical regions with respect to B0. A. When the longest axis of the region is oriented parallel to B0, a paramagnetic shift that increases the local field occurs. B. When the longest axis of the region is oriented perpendicular to B0, a diamagnetic effect that tends to decrease the field is generated.

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:

Δ B = 2 π Δ χ o Δ Y c a p H c t B o ( exp(- X 2 2 σ x 2 - Y 2 2 σ y 2 Z 2 2 σ z 2 ) σ i = F W H M i 2 2 ln ( 2 ) E11

where ΔYcap is the change in oxygenation in the capillaries (0.08) due to BOLD, X, Y and Z are the Cartesian coordinates in 3D space, FWHMi expresses the spatial extent of the area in terms of the full width half maximum value for the coordinates in the i-direction, where i=X,Y,Z; and the other parameters are the same as in [Eq.10] and in Fig. 7. We varied the size of the activated region between 5-40mm (FWHM of the shortest length axis) and the strength of the BOLD response. Values for the BOLD induced magnetic field change were close to those observed above for the single vessel model (± 0.05-0.25 μT). The results from these simulations are shown in Figs. 8 and 9. We found that both the spatial extent of the activated area and the strength of the BOLD response influence the size of phase change. For instance, magnetic field shifts of ±0.1 μT lead to phase shifts in the range of 0.01-0.04 radian, depending on the spatial extent of the region. As shown in Fig. 9, the greatest phase effect was found in small areas with a strong BOLD response. The size of these expected phase changes are thus about two orders of magnitude lower than those found for fully sampled large draining vessels, despite similar size of the BOLD induced magnetic field change. This fact highlights the difference between biophysical models.

Figure 9.

Simulations of the phase effect with positive (red) and negative (blue) sign generated by the BOLD response in a brain area, with an ellipsoid Gaussian shape. As in Fig. 9, the longest axis is oriented along B0. A: Variation of the phase effect with the spatial extent of the activated area (indicated as the FWHM of the shortest length axes, see Eq. 11 and Fig.8). The BOLD response is fixed (±0.10 μT). The phase effect is greatest when the size of the region is small, but diminishes as the size of the activated region increases. B: Variation of the phase effect with the strength of the BOLD responses. The spatial extent of the activated brain area is fixed (FWHM=10mm). The positive (negative) phase increases (decreases) with increasing BOLD field changes. The arrow indicates equivalent points in Fig. 9 A and 9B where the FWHM is 10mm and the BOLD response has a field change of ±0.10 μT.

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:

φ true = φ obs ± 2 π k E12

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).

Figure 10.

Fourier analysis of the heart beat and the respiratory cycle (A) and their relation to the MRI signal in magnitude (B) and phase (C-F) images. A: Respiration occurs between 0.12 and 0.25 Hz, and the first and second harmonics of the heart-beat around 1.2 and 2.4 Hz, respectively. B: In magnitude images, respiration induced fluctuationscause the typical 1/f decrement at low frequencies, but instabilities are also prominent around the first cardiac harmonic. Note the relative absence of disturbances around 2Hz, caused by the scanner Helium pump and visible in most of the phase time-series. Results are shown for gray matter (solid line), white matter (dotted line), and CSF (dashed line). C: temporal unwrap and detrending of the phase time-series (Hagberg et al., 2008); D: spatial unwrapping followed by temporal unwrapping and detrending, (Hagberg et al., 2008); E: Reference phase method (Tomasi and Caparelli, 2007); F: k-space filtering (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, tSNR based on the factor λM that expresses the proportionality between the magnitude signal and the physiological noise. Basically this model implies that the available temporal signal-to-noise level reaches a maximum limiting value for high signal-to-noise-ratios at a single time point, SNR0, calculated from the average signal in each voxel across the GE-EPI timecourse, S, divided by the Rayleigh corrected noise level, N:

S N R o = S 1.53 N t S N R = S N R o 1 + λ M 2 S N R o 2 E13

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):

t S D φ = 1 t S N R = 1 + λ φ 2 S N R 0 2 S N R 0 E14

where λM< λφ; λM= λφ; or λM> λφ, respectively for the three noise models.

The tSNR obtained for the magnitude signal and 1/tSDφ images obtained after different phase processing methods are shown in Fig. 11 and illustrates that the phase noise is generally greater than the magnitude noise, but can be efficiently reduced by adequate post-processing methods. (Fig. 11, see also Hagberg et al., 2012). Spatial high-pass filtering currently represents the best approach for removal of unwanted phase fluctuations. If adequately adapted, similar or even reduced λφ values are obtained for the phase compared to magnitude data. In general, the amount of filtering must be adapted to the field strength of the scanner, the echo time and voxel size of the GE-EPI sequence. For instance, at 3T, a TE of 30ms, and a voxel size of 3x3x3mm3 we found that a Gaussian filter with a 20mm FWHM suffices to reach the a similar stability in both phase and magnitude time series (Hagberg et al., 2012). Filters with a FWHM close to the voxel size, reduces the phase noise to less than half the magnitude noise. For most other post-processing methods, including spatio-temporal unwrapping, regression based on the phase evolution central k-space point (NVR), RETROICOR, and TOAST, the phase noise was always greater than the magnitude noise.

Figure 11.

Temporal stability in magnitude (A) and phase (B-F) images. The phase images were post-processed by different methods prior to calculation of the inverse of the standard deviation of the phase across the time series (see Eq. 14). (B) RETROICOR; (C) TOAST; (D) noise regression based on the phase of the central k-space point; (E) combined approach based on both thecentral k-space point and RETROICOR; (F) : k-space High pass filtering of the phase data by a Homodyne filter based on a Gaussian kernel in k-space (Hagberg et al., 2012).


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.

Tissue Raw magnitude data RETROICOR NVR RETROICOR and NVR combined
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

Table 1

Stability in magnitude images in terms of the physiological noise parameter λ M (see Eq. 13) for raw images and after noise-regression.

Figure 12.

Increase in temporal stability of the magnitude signal for different magnitude post-processing methods. The data shown is the percent increase in temporal signal-to-noise ratio with respect to no post-processing for (A)after NVR; noise regression based on the evolution at the centre of k-space; (B) after RETROICOR; and (C) after a combined approach based on both NVR and RETROICOR methods. Data were obtained during a 5min acquisition of resting-state fMRI data. A sagittal plane showing the results for all acquired transverse slices is shown.

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:

a b s ( y t ) = X β + η E15

For complex numbers, GLM can be rewritten as (Rowe, 2005):

[ Re ( y t ) Im ( y t ) ] = [ cos φ t 0 0 sin φ t ] [ x 0 0 x ] [ β r e β i m ] + [ η r e t η i m t ] E16

where yt is the complex valued fMRI time course, φ t is the phase at each time point, x, X are the design matrices for the complex vector and the magnitude signals, respectively; βtre, βtim, β the parameter estimates obtained after model fitting; and ηtre, ηtim, η the noise of the real, imaginary and magnitude signals.

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 make sure that we adequately handle unwanted phase effects, different approaches for removal of these were explored and their impact on the final activation maps were investigated.

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).

Figure 13.

Activation patterns for a motor task obtained from a) magnitude (M) and b-f) phase (P) images. Unwanted phase effects were removed by five different homodyne (spatial high-pass) filters: b) Hamming window; c) Gaussian 20mm; d) Gaussian 10mm; e) Gaussian 5mm and f) k-space spline window. The statistical maps have been thresholded at a statistical significance level of p<0.001(uncorrected for multiple comparisons). Color bars indicate statistical T-value ranges.

Cont # voxels in cluster Tvalue Zvalue p-value MNI coordinates
Hamming n=32 Pos
Clmax 659
34, -26, 46
38, -6, 42
Gauss FWHM
Clmax 792
32, -26, 48
38, -6, 40
Gauss FWHM
Clmax 841
Clmax 238
42, -20 22
38, -6, 44
Gauss FWHM
Clmax 1064
Clmax 381
34, -26, 44
-40, -30, 52
k-spline filter Pos
Clmax 606
Clmax 213
42, -20, 18
-40, -30, 42
Raw data Pos Clmax 950 15.06 6.37 FWE=0.5 52, -14, 54

Table 2

Results obtained in a second level analysis. T and Z statistical values and size of the activated clusters (agglomerate of activated voxels) for magnitude and phase images are reported. Corresponding p-values and Montreal Neurological Institute (MNI) coordinates of their centre are shown for positive (Pos) and negative (Neg) phase changes, and for the positive magnitude activation.

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.

Figure 14.

Timeseries extracted from brain regions with significant, stimulus correlated signal changes in three different subjects (SS, BM, and OF). Positive (left) and negative (right) phase signals are shown, overlaid on the magnitude signal (Amp, green). The phase signal clearly has a lower contrast-to-noise ratio than the magnitude signal. The temporal evolution of the phase and magnitude timeseries are similar.

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 B0 signal components currently represent the best approach for achieving this goal. We show that phase changes consistent with a dipolar activation pattern caused by an extended cortical region can be detected. Future work must address the question of how the BOLD related phase activation patterns can be incorporated in the statistical analysis pipeline, in order to improve the specificity and the sensitivity of fMRI.



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).


  1. 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 Neuroimage 49 4 3149 3160
  2. 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. 3. Barmet C. De Zanche N. Pruessmann K. P. 2008 Spatiotemporal magnetic field monitoring for MRMagn Reson Med. 60 1 187 97
  4. 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. 5. Boxerman J. L. Hamberg L. M. Rosen B. R. Weisskoff R. M. 1995MR contrast due to intravascular magnetic susceptibility perturbations Magn Reson Med.;34 4 555 566
  6. 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. 7. Bretthorst G. L. 2008 Automatic phasing of MR images. Part II: Voxel-wise phase estimationJournal of Magnetic Resonance 191 193
  8. 8. Chen Z. Calhoun V. D. Magnitude and Phase Behavior of Multiresolution BOLD Signal 2010Conc Magn Reson B 37B(3) 129 EOF 145 EOF
  9. 9. Cusack R. Papadikis N. 2002 New robust 3-D phase unwrapping algorithms: application to magnetic field mapping and undistorting echoplanar images. Neuroimage 16 754 764
  10. 10. Durrant C. J. Hertzberg M. P. Kuchel P. W. 2003 Magnetic susceptibility: further insights into macroscopic and microscopic fields and the sphere of LorentzConcMagn. Reson. 18A(1): 72- 95.
  11. 11. Feng Z. Caprihan A. Blagoev K. B. Calhoun V. D. 2009 Biophysical modeling of phase changes in BOLD fMRI. Neuroimage 47 2 540 548
  12. 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. 13. Glover G. H. 1999 Deconvolution of Impulse Response in Event-Related BOLD fMRI. NeuroImage 9 416 429
  14. 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. 15. Gudbjartsson H. Patz S. 1995 The Rician distribution of noisy MRI data.Magn Reson Med 34(6) 910 EOF 4 EOF
  16. 16. Hagberg G. E. Bianciardi M. Brainovich V. Cassarà A. M. Maraviglia B. 2008 The effect of physiological noise in phase functional magnetic resonance imagingfrom blood oxygen level-dependent effects to direct detection of neuronal currentsMagn Reson Imaging 26(7), 1026 EOF 1040 EOF
  17. 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. 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. Neuroimage 59 4 3748 3761
  19. 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). Neuroimage 44 3 742 52
  20. 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. 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. 22. Jenkinson M. 2003Fast, automated N-dimensional phase-unwrapping algorithm. Magn Reson Med. 49(1), 193-197.
  23. 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. 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. 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. 26. Logothetis N. K. 2008 What we can do and what we cannot do with fMRI. Nature 453 7197 869 78
  27. 27. Mansfield P. 1977Multi-planar image formation using NMR spinechoes. 1. Phys. C. 10:L55L58.
  28. 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 susceptibilityConcepts Magn Reson Part B 25B 65 78
  29. 29. Menon R. 2002Postacquisition Suppression of Large-Vessel BOLD Signals in High-Resolution fMRI Magn Reson Med 47(1), 1-9.
  30. 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 fieldsMagn Reson Med. 14 1 68 78
  31. 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. 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 J 64 3 803 812
  33. 33. Pfeuffer J. Van de Moortele P.-F, Ugurbil K, Hu X, Glover GH, 2002 Correction of Physiologically Induced Global Off-Resonance Effects in Dynamic Echo-Planar and Spiral Functional Imaging.Magn Reson Med. 10 344 353
  34. 34. Petridou N. Schäfer A. Gowland P. Bowtell R. 2009 Phase vs. magnitude information in functional magnetic resonance imagingtime series: toward understanding the noise.Magn Reson Imaging 27(8) 1046 EOF 57 EOF
  35. 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. 36. Rowe D. B. Logan B. R. 2004 A complex way to compute fMRI activation Neuroimage 1078 EOF 1092 EOF
  37. 37. Rowe D. B. 2005 Modeling both the magnitude and phase of complex-valued fMRI data Neuroimage 25 1310 1324
  38. 38. Schmitt F. Stehling M. K. Turner R. 1998 Echo-Planar Imaging: Theory, Technique and Application.Springer Verlag
  39. 39. Van Gelderen P. De Zwart J. A. Starewicz P. Hinks R. S. Duyn J. H. 2007Real-time shimmingto compensate for respiration induced Bo fluctuations. Magn Reson Med 57 362 368
  40. 40. Tomasi D. G. Caparelli E. G. 2007Macrovascular cobtribution in activation patterns of working memory. J Cereb Blood Flow Metab 27 33 42
  41. 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. 42. Yablonskiy D. A. Haacke E. M. 1994Theory of NMR signal behaviour in magnetically inhomogeneous tissues: the static dephasing regime. Magn Reson Med 32 6 749 763
  43. 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 Med 57 520 527

Written By

Gisela Hagberg and Elisa Tuzzi

Published: 31 May 2014