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

## 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 O_{2} and CO_{2}. 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

^{1}H, 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, _{0}, 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 * M*that 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,

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

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

B

_{1}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*.

B

_{1}field is applied orthogonally with respect to

B

_{0}(Fig. 2).

The induced transition to the upper energy level gives rise to a spin coherence, bringing the longitudinal equilibrium magnetization _{z}into the transverse plane x-y, _{xy}, 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 _{1}, or a so-called 90° pulse. When the _{1}90° pulse is interrupted,* M*starts to precess around the B

_{0}z-axis at the Larmor frequency, whilst re-establishing its initial equilibrium value.This causes the gradual recovery of the longitudinal component

M

_{z}and the complete extinction of M

_{xy}. The temporal variation of the transverse component,

M

_{xy}, 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,

B

_{0}, perturbed by an oscillating field,

B

_{1}, 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, M_{xy}, is described by the characteristic time constant, T_{2}, also known as the spin-spin relaxation time, while the characteristic time constant T_{1}, known as spin-lattice relaxation time, rules the recovery of the transverse component, M_{z}, 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 _{2}* *,*or the effective transverse relaxation time:

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 _{2}* **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

B

_{o}(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:

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.

### 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 _{0}. These gradients are applied in the three spatial directions (x,y,z), defining respectively the * slice-selection gradient*,

G

_{z}, the

phase-encode gradient, G

_{y}

*and the*,

*(also*frequency-encode

*)*readout

*,*gradient

G

_{x}. 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 mm

^{3}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 mm

^{3}(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*α

B

_{0}), 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

G

_{x}in Fig.3). Before each readout gradient, a brief phase-encode blip (

G

_{y}) 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.

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 B_{0} 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 _{xy}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 _{x}and the _{y}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 _{x}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):

In general, the phase evolves linearly with the local magnetic field _{loc}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:

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

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 (* ev*) - and intravascular (

*) space (in cgs units),*iv

*the oxygenation level of the blood,*Y

*is the haematocrit level,*Hct

B

_{0}is the static magnetic field strength,

*the vessel radius and*R

*the distance from the centre of the vessel cylinder, and*r

*are the angles defined in Fig 7.*θ, Ψ

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

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 _{cap}is the change in oxygenation in the capillaries (0.08) due to BOLD, * X*,

*and*Y

*are the Cartesian coordinates in 3D space,*Z

FWHM

_{i}expresses the spatial extent of the area in terms of the full width half maximum value for the coordinates in the

*-direction, where*i

*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.*i=X,Y,Z;

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 2^{nd} 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 B_{0} 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, * 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, SNR

_{0}, calculated from the average signal in each voxel across the GE-EPI timecourse,

*, divided by the Rayleigh corrected noise level,*S

*:*N

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 _{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 3x3x3mm

^{3}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.

## 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 _{t}is the complex valued fMRI time course, * x, X*are the design matrices for the complex vector and the magnitude signals, respectively;

β

_{t}

^{re},

β

_{t}

^{im},

*the parameter estimates obtained after model fitting; and*β

η

_{t}

^{re},

η

_{t}

^{im},

*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 ** m**ake 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).

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 | p_{u}<0.001p _{u}<0.001 | 34, -26, 46 38, -6, 42 |

Gauss FWHM 20mm | Pos Neg | Clmax 792 92 | 8.04 5.74 | 4.93 4.11 | p_{u}<0.001p _{u}<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 | p_{u}<0.01p _{u}<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 | p_{u}<0.01p _{u}<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 | p_{u}<0.001p _{u}<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 _{0} 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.

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